shoot.cpp

Visual C++ 常用数值算法集 源代码

void shoot(int nvar, double v[], double delv[], int n2, double x1,
		   double x2, double eps, double h1, double hmin, 
		   double f[], double dv[])
{
    double y[21], dfdv[401],xp[201], yp[2001];
	int ypn=100;
	int indx[21]; 
    load(x1, v, y);
	int i,iv,nok,nbad;
    odeint(y, nvar, x1, x2, eps, h1, hmin, nok, nbad, xp, yp,ypn);
    score(x2, y, f);
	double sav;
    for (iv = 1; iv<=n2; iv++)
	{
        sav = v[iv];
        v[iv] = v[iv] + delv[iv];
        load(x1, v, y);
        odeint(y,nvar,x1,x2,eps,h1,hmin, nok, nbad, xp, yp,ypn);
        score(x2, y, dv);
        for (i = 1; i<=n2; i++)
		{
            dfdv[(i-1)*20+iv] = (dv[i] - f[i]) / delv[iv];
        }
        v[iv] = sav;
    }
    for (iv = 1; iv<=n2; iv++)
	{
        dv[iv] = -f[iv];
    }
	int det;
    ludcmp(dfdv, n2, 20,indx, det);
    lubksb(dfdv, n2, 20,indx, dv);
    for (iv = 1; iv<=n2; iv++)
	{
        v[iv] = v[iv] + dv[iv];
    }
}