shootf.cpp

实用的算法

void shootf(int nvar, double v1[], double v2[], double delv1[],
			double delv2[], int n1, int n2, double x1, double x2,
			double xf, double eps, double h1, double hmin,
			double f[], double dv1[], double dv2[])
{
    double y[21], f1[21], f2[21], dfdv[401], xp[201], yp[2001];
	int indx[21];
	int nok,nbad,ypn=200;
    load1(x1, v1, y);
    odeint(y, nvar, x1, xf, eps, h1, hmin, nok, nbad, xp, yp,ypn);
    score(xf, y, f1);
    load2(x2, v2, y);
    odeint(y, nvar, x2, xf, eps, h1, hmin, nok, nbad, xp, yp,ypn);
    score(xf, y, f2);
    int iv,i,j = 0;
	double sav;
    for (iv = 1; iv<=n2; iv++)
	{
        j = j + 1;
        sav = v1[iv];
        v1[iv] = v1[iv] + delv1[iv];
        load1(x1, v1, y);
        odeint(y,nvar,x1,xf,eps, h1, hmin, nok, nbad, xp, yp,ypn);
        score(xf, y, f);
        for (i = 1; i<=nvar; i++)
		{
            dfdv[(i-1)*20+j] = (f[i] - f1[i]) / delv1[iv];
        }
        v1[iv] = sav;
    }
    for (iv = 1; iv<=n1; iv++)
	{
        j = j + 1;
        sav = v2[iv];
        v2[iv] = v2[iv] + delv2[iv];
        load2(x2, v2, y);
        odeint(y,nvar,x2,xf,eps,h1, hmin, nok, nbad, xp, yp,ypn);
        score(xf, y, f);
        for (i = 1; i<=nvar; i++)
		{
            dfdv[(i-1)*20+j] = (f2[i] - f[i]) / delv2[iv];
        }
        v2[iv] = sav;
    }
    for (i = 1; i<=nvar; i++)
	{
        f[i] = f1[i] - f2[i];
        f1[i] = -f[i];
    }
	int det;
    ludcmp(dfdv, nvar, 20, indx, det);
    lubksb(dfdv, nvar, 20, indx, f1);
    j = 0;
    for (iv = 1; iv<=n2; iv++)
	{
        j = j + 1;
        v1[iv] = v1[iv] + f1[j];
        dv1[iv] = f1[j];
    }
    for (iv = 1; iv<=n1; iv++)
	{
        j = j + 1;
        v2[iv] = v2[iv] + f1[j];
        dv2[iv] = f1[j];
    }
}