打靶法.cpp

打靶法c程序

#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>
#include <string>
#include <process.h>
# include  <malloc.h>
# include  <stdio.h>
# include  <float.h>
double x;
double hh,h6,xh;//,dydx[1],yout[1];
int i,j,n;
double d,h=0.0;
double y[1],z[1],dydx[1],dzdx[1],yout[1],zout[1],yding[1],zding[1];
void derivs(double x,double y[],double z[],double dydx[],double dzdx[])
{
	dydx[0]=z[0];//y'=z[0]  
	dzdx[0]=1.5*y[0]*y[0];//z'=f(x,y,z)=1.5*y*y
}
//  子过程RK4

void rk4(int n,double h,double x,double y[],double z[],double dydx[],double dzdx[],double yout[],double zout[])
{
	
	double zt[2],dzt[2],dzm[2],dzf[2],yt[2];
 	hh=h*0.5;
	h6=h/6;
	xh=x+hh;
	derivs(x,y,z,dydx,dzdx);
	for(i=0;i<n;i++)
	{
		
		zt[i]=z[i]+hh*dzdx[i];//z(n)+h/2*L1
		yt[i]=y[i]+hh*z[i];
	}
	derivs(xh,yt,zt,dydx,dzt);//x+h/2,y(n)+h/2*zn,zn+h/2*L1,L2
	for(i=0;i<n;i++)
	{
		yt[i]=y[i]+hh*z[i]+h*h/4*dzdx[i]; 
		zt[i]=z[i]+hh*dzt[i];//z(n)+h/2*L2
	}
	derivs(xh,yt,zt,dydx,dzm);//dzm=L3
	for(i=0;i<n;i++)
	{
		yt[i]=y[i]+h*z[i]+h*h/2*dzt[i];//dzt=L2
		zt[i]=z[i]+h*dzm[i];//z(n)+h/2*L3
		//dym[i]=dyt[i]+dym[i];///K2+k3
	}
	derivs(x+h,yt,zt,dydx,dzf);/////x+h*******dzf=L4
	for(i=0;i<n;i++)
	{
		yout[i]=y[i]+h*z[i]+h*h/6*(dzdx[i]+dzt[i]+dzm[i]);
	 	y[i]=yout[i];
		zout[i]=z[i]+h6*(dzdx[i]+2.0*dzt[i]+2.0*dzm[i]+dzf[i]); //h/6=h6 
	//	cout<<zout[i]<<endl;//yout从这里输出,切记！！！//	 cout<<h6;
	 	z[i]=zout[i];

	}
	 for(i=0;i<n;i++)
	{
		
		dzm[i]=0.0;
		dzt[i]=0.0;
		yt[i]=0.0;
		zt[i]=0.0;
		dzf[i]=0.0;
	 } 
}
void main()
{
	double F[1],Fs[1],zjia[1],dlts=0.00001;
	double Fpie[1],dltF;
	n=1;
	x=0;
	y[0]=4;
	z[0]=-1.0000;
	//dydx[0]=z[0];//y'=z[0]  
	//dzdx[0]=1.5*y[0]*y[0];//z'=f(x,y,z)=1.5*y*y
do
	{		
	for(j=0;j<5;j++)
		{ 
			if(j<1)
			{
			zding[0]=z[0];//zding[0]=z[0.0]
			yding[0]=y[0];//yding[0]=y[0.0]
	   		}
		h=0.2+0.2*j;
		rk4(n,h,x,y,z,dydx,dzdx,yout,zout);
//////		cout<<"yout["<<h<<"]="<<yout[0]<<endl;//如果输出的yout[0]，也可以!!!!
		}
	F[0]=yout[0]-1.000000;//输出F(s);
	zjia[0]=zding[0]+dlts;//zjia[0]=s(i)+dlts(i)	
	dydx[0]=zjia[0];//y'=z[0]  
	dzdx[0]=1.5*yding[0]*yding[0];//z'=f(x,y,z)=1.5*y*y
	for(j=0;j<5;j++)
		{            
		h=0.2+0.2*j;
		rk4(n,h,x,yding,zjia,dydx,dzdx,yout,zout);
	 //	cout<<"yout["<<h<<"]="<<yout[0]<<endl;
		}
	Fpie[0]=yout[0]-1.000000;
	dltF=(Fpie[0]-F[0])/dlts;
	z[0]=zding[0]-F[0]/dltF;//输出新的z[0] 也就是s(i+1)
	Fs[0]=F[0];
 	 //cout<<Fpie[0]<<endl;
    //cout<<yout[0]<<endl;
	//	F[0]=Fpie[0];
	y[0]=4;
	cout<<"F[0]="<<F[0];
	cout<<endl;
	}while(F[0]>0.01);//	&&F[0]<-0.00001);
	cout<<"ok!"<<endl;
	cout<<"z[0]="<<z[0]+Fs[0]/dltF<<endl;
	cout<<"F[0]="<<F[0]<<endl;
	cout<<"yout[0]="<<yout[0];
}