crank_nicolson.h

四种方法求积分：runge-kutta法

vector<double> crank_nicolson(double a,			//a
							  vector<double> u0,//初始条件
							  vector<double> us,//边界条件 
							  vector<double> um,//边界条件
							  double x0,		//x0
							  double t0,		//t0
							  double r,			//a*tao/h/h
							  double tao,		//tao	
							  int total)		//(tn-t0)/tao
{
	vector<double> u, utemp;

	int n = u0.end() - u0.begin();
	
	vector<double> aa, b, c, d;

	for(int i = 0; i < n - 3; i++)
	{
		aa.push_back(- r / 2);
		c.push_back(- r / 2);
	}

	for(i = 0; i < n - 2; i++)
	{
		b.push_back(1 + r);
		d.push_back(0);
	}
	

	for(i = 0; i < n; i++)
		u.push_back(u0[i]);


	for(int k = 0; k < total - 1; k++)
	{
		u.push_back(us[k]);


		d[0] = tao * fn_chafen(sqrt(tao * a / r) + x0, t0 + (k + 0.5) * tao)
			+ r / 2 * (arpha(t0 + (k) * tao) + arpha(t0 + (k + 1) * tao))
			+ u[k * n + 1] * (1 - r) + u[k * n + 2] * r / 2;
		for(i = 1; i < n - 3; i++)
		{
			d[i] = tao * fn_chafen(sqrt(tao * a / r) * (i) + x0, t0 + (k + 0.5) * tao)
				+ u[k * n + i + 1] * (1 - r) 
				+ (u[k * n + i] + u[k * n + i + 2]) * r / 2;
		}
		d[n - 3] = tao * fn_chafen(sqrt(tao * a / r) * (n - 3) + x0, t0 + (k + 0.5) * tao)
			+ r / 2 * (beta(t0 + (k) * tao) + beta(t0 + (k + 1) * tao))
			+ u[k * n + n - 2] * (1 - r)
			+ u[k * n + n - 3] * r / 2;

		utemp = zhuigang(aa, b, c, d);

		
		for(i = 0; i < n - 2; i++)
			u.push_back(utemp[n - 3 - i]);

		u.push_back(um[k]);
	}

	return u;
}