difeq.cpp

实用的算法

void difeq(int k, int k1, int k2, int jsf, int is1, int isf,
		   int indexv[], int ne, double s[], int nsi, int nsj,
		   double y[], int nyj, int nyk)
{
    int m = 41;
    if (k == k1)
	{
        if (n + mm % 2 == 1)
		{
            s[2*nsj+3 + indexv[1]] = 1.0;  //方程(15-37)
            s[2*nsj+3 + indexv[2]] = 0.0;
            s[2*nsj+3 + indexv[3]] = 0.0;
            s[2*nsj+jsf] = y[1];  //方程(15-31)
		}
        else
		{
            s[2*nsj+3 + indexv[1]] = 0.0;  //方程(15-37)
            s[2*nsj+3 + indexv[2]] = 1.0;
            s[2*nsj+3 + indexv[3]] = 0.0;
            s[2*nsj+jsf] = y[nyk+1];  //方程(15-31)
        }
	}
    else
	{
		if (k > k2)
		{
			s[3 + indexv[1]] = -(y[2*nyk+m] - c2) / (2.0 *
				         double(mm + 1.0));    //方程(15-38)
			s[3 + indexv[2]] = 1.0;
			s[3 + indexv[3]] = -y[m] / (2.0 * double(mm + 1.0));
			s[jsf] = y[nyk+m] - (y[2*nyk+m] - c2) * y[m] / (2.0 *
				         double(mm + 1.0));   //方程(15-32)
			s[nsj + 3 + indexv[1]] = 1.0;     //方程(15-39)
			s[nsj + 3 + indexv[2]] = 0.0;
			s[nsj + 3 + indexv[3]] = 0.0;
			s[nsj + jsf] = y[m] - anorm;      //方程(15-33)
		}
		else
		{
			s[indexv[1]] = -1.0;              //方程(15-34)
			s[indexv[2]] = -0.5 * h;
			s[indexv[3]] = 0.0;
			s[3 + indexv[1]] = 1.0;
			s[3 + indexv[2]] = -0.5 * h;
			s[3 + indexv[3]] = 0.0;
			double temp = h/(1.0 - pow((x[k]+x[k - 1]),2)* 0.25);
			double temp2= 0.5 * (y[2*nyk+k] + y[2*nyk+k - 1]) - c2
				          * 0.25 * pow((x[k] + x[k - 1]) , 2);
			s[nsj+indexv[1]] = temp * temp2 * 0.5;  //方程(15-35)
			s[nsj+indexv[2]] = -1.0 - 0.5 * temp * double(mm + 1.0)
				               * (x[k] + x[k - 1]);
			s[nsj+indexv[3]] = 0.25 * temp * (y[k] + y[k - 1]);
			s[nsj+3 + indexv[1]] = s[nsj+indexv[1]];
			s[nsj+3 + indexv[2]] = 2.0 + s[nsj+indexv[2]];
			s[nsj+3 + indexv[3]] = s[nsj+indexv[3]];
			s[2*nsj+indexv[1]] = 0.0;  //方程(15-36)
			s[2*nsj+indexv[2]] = 0.0;
			s[2*nsj+indexv[3]] = -1.0;
			s[2*nsj+3 + indexv[1]] = 0.0;
			s[2*nsj+3 + indexv[2]] = 0.0;
			s[2*nsj+3 + indexv[3]] = 1.0;
			s[jsf] = y[k] - y[k - 1] - 0.5 * h * (y[nyk+k] +
				       y[nyk+k - 1]);  //方程(15-26)
			double dum = (x[k] + x[k - 1]) * 0.5 * double(mm + 1.0)
				        * (y[nyk+k] + y[nyk+k - 1]);    //方程(15-27)
			dum = dum - temp2 * 0.5 * (y[k] + y[k - 1]);
			s[nsj+jsf] = y[nyk+k] - y[nyk+k - 1] - temp * dum;
			s[2*nsj+jsf] = y[2*nyk+k] - y[2*nyk+k - 1]; //方程(15-30)
		}
	}
}