rkqc.cpp

此程序为VC++常用数值算法这本书附赠的光盘中包含了本书中全部的源代码

void rkqc(double y[], double dydx[], int n, double &x, double htry, double eps, double yscal[], double &hdid, double &hnext)
{
	int i,flag;
	double one,safety,errcon,fcor,ytemp[11],ysav[11];
	double dysav[11],xsav,h,hh,pgrow,pshrnk,errmax;
    one = 1.0;
    safety = 0.9;
    errcon = 0.0006;
	fcor = 0.066667;
    pgrow = -0.2;
    pshrnk = -0.25;
    xsav = x;
    for( i = 1; i<=n; i++)
	{
        ysav[i] = y[i];
        dysav[i] = dydx[i];
    }
    h = htry;
    do 
	{
        hh = 0.5 * h;
        rk4(ysav, dysav, n, xsav, hh, ytemp);
        x = xsav + hh;
        derivs(x, ytemp, dydx);
        rk4(ytemp, dydx, n, x, hh, y);
        x = xsav + h;
        if( x == xsav)
		{
            cout<<" stepsize not significant in rkqc."<<endl;
            _c_exit();
        }
        rk4(ysav, dysav, n, xsav, h, ytemp);
        errmax = 0.0;
        for( i = 1; i<=n; i++)
		{
            ytemp[i] = y[i] - ytemp[i];
            if (errmax < fabs(ytemp[i] / yscal[i]))
			{
                errmax = fabs(ytemp[i] / yscal[i]);
			}
        }
        errmax = errmax / eps;
        if (errmax > one)
		{
            h = safety * h * pow(errmax , pshrnk);
            flag = 1;
		}
        else
		{
            hdid = h;
            if( errmax > errcon )
			{
                hnext = safety * h * pow(errmax , pgrow);
			}
            else
			{
                hnext = 4.0 * h;
			}
            flag = 0;
		}
	}
	while( flag == 1);
    for (i=1; i<=n; i++)
	{
        y[i] = y[i] + ytemp[i] * fcor;
    }
	for(i=1; i<=10; i++)
	{
		dysav[i]=0.0;
		ysav[i]=0.0;
		ytemp[i]=0.0;
	}
}