shoot.c

Numerical Recipes 是国际公认的高水平的、关于数值计算的书

void shoot(nvar,v,delv,n2,x1,x2,eps,h1,hmin,f,dv)
int nvar,n2;
float v[],delv[],x1,x2,eps,h1,hmin,f[],dv[];
{
	int nok,nbad,iv,i,*indx,*ivector();
	float sav,det,*y,**dfdv,**matrix(),*vector();
	void odeint(),ludcmp(),lubksb(),derivs(),rkqc(),
		free_matrix(),free_vector(),free_ivector();
	void load();	/* ANSI: void load(float,float *,float *); */
	void score();	/* ANSI: void score(float,float *,float *); */

	y=vector(1,nvar);
	indx=ivector(1,nvar);
	dfdv=matrix(1,n2,1,n2);
	load(x1,v,y);
	odeint(y,nvar,x1,x2,eps,h1,hmin,&nok,&nbad,derivs,rkqc);
	score(x2,y,f);
	for (iv=1;iv<=n2;iv++) {
		sav=v[iv];
		v[iv] += delv[iv];
		load(x1,v,y);
		odeint(y,nvar,x1,x2,eps,h1,hmin,&nok,&nbad,derivs,rkqc);
		score(x2,y,dv);
		for (i=1;i<=n2;i++)
			dfdv[i][iv]=(dv[i]-f[i])/delv[iv];
		v[iv]=sav;
	}
	for (iv=1;iv<=n2;iv++) dv[iv] = -f[iv];
	ludcmp(dfdv,n2,indx,&det);
	lubksb(dfdv,n2,indx,dv);
	for (iv=1;iv<=n2;iv++) v[iv] += dv[iv];
	free_matrix(dfdv,1,n2,1,n2);
	free_ivector(indx,1,nvar);
	free_vector(y,1,nvar);
}