odeint.cpp

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

void odeint(double ystart[5], int &nvar, double &x1, double &x2, double &eps, double &h1, double &hmin, double &nok, double &nbad, double xp[201], double yp[11][201])
{
	int i,nstp,maxstp;
	double yscal[11], y[11], dydx[11],two,zero,tiny,x,h,xsav,hdid,hnext;
    maxstp = 10000;
    two = 2.0;
    zero = 0.0;
    tiny = 1e-30;   
    x = x1;
    h = fabs(h1) * sgn(x2 - x1);
    nok = 0;
    nbad = 0;
    kount = 0;
	hdid=0.0;
	hnext=0.0;
    for (i = 1; i<=nvar; i++)
	{
        y[i] = ystart[i];
    }
    if (kmax > 0)
		xsav = x - dxsav * two;
    for( nstp = 1; nstp<=maxstp; nstp++)
	{
        derivs(x, y, dydx);
        for( i = 1; i<=nvar; i++)
		{
            yscal[i] = fabs(y[i]) + fabs(h * dydx[i]) + tiny;
        }
        if( kmax > 0 )
		{
            if( fabs(x - xsav) > fabs(dxsav))
			{
                if( kount < kmax - 1)
				{
                    kount = kount + 1;
                    xp[kount] = x;
                    for( i = 1; i<=nvar; i++)
					{
                        yp[i][kount] = y[i];
                    }
                    xsav = x;
                }
            }
        }
        if( (x + h - x2) * (x + h - x1) > zero )
		{
			h = x2 - x;
		}
        rkqc(y, dydx, nvar, x, h, eps, yscal, hdid, hnext);
        if( hdid == h)
		{
            nok = nok + 1;
		}
        else
		{
            nbad = nbad + 1;
		}
        if( (x - x2) * (x2 - x1) >= zero )
		{
            for (i = 1;i<=nvar;i++)
			{
                ystart[i] = y[i];
            }
            if( kmax != 0 )
			{
                kount = kount + 1;
                xp[kount] = x;
                for( i = 1; i<=nvar; i++)
				{
                    yp[i][kount] = y[i];
                }
            }
			for(i=1; i<=10; i++)
			{
				dydx[i]=0.0;
				y[i]=0.0;
				yscal[i]=0.0;
			}
            goto over;
        }
        if (fabs(hnext) < hmin)
		{
            cout<<" stepsize smaller than minimum."<<endl;
            _c_exit();
        }
        h = hnext;
    }
    cout<<" too many steps."<<endl;
over:;
}