龙格库塔算法解微分方程组c语言.docxVIP

/*************************************************************************** This program is to solve the initial value problem of following system of differential equations: dx/dt=x+2*y,x(0)=0, dy/dt=2*x+y,y(0)=2, x(0.2) and y(0.2) are to be calculated ****************************************************************************/ #includestdio.h #includemath.h #define steplength 0.1 //步?长?è可¨|根¨′据Y需¨¨要°a调ì??整?; #define FuncNumber 2 //FuncNumber为a微?é分¤?方¤?程¨?的ì?数oy目?; void main() { float x[200],Yn[20][200],reachpoint;int i; x[0]=0;Yn[0][0]=0;Yn[1][0]=2; //初?值|ì条??件t; reachpoint=0.2; //所¨′求¨?点ì?可¨|根¨′据Y需¨¨要°a调ì??整?; void rightfunctions(float x ,float *Auxiliary,float *Func); void Runge_Kutta(float *x,float reachpoint, float(*Yn)[200]); Runge_Kutta(x ,reachpoint, Yn); printf(x ); for(i=0;i=(reachpoint-x[0])/steplength;i++) printf(%f ,x[i]); printf(\nY1 ); for(i=0;i=(reachpoint-x[0])/steplength;i++) printf(%f ,Yn[0][i]); printf(\nY2 ); for(i=0;i=(reachpoint-x[0])/steplength;i++) printf(%f ,Yn[1][i]); getchar(); } void rightfunctions(float x ,float *Auxiliary,float *Func)//当ì?à右?¨°方¤?程¨?改?变à?时o?à，ê?需¨¨要°a改?变à?; { Func[0]=Auxiliary[0]+2*Auxiliary[1]; Func[1]=2*Auxiliary[0]+Auxiliary[1]; } void Runge_Kutta(float *x,float reachpoint, float(*Yn)[200]) { int i,j; float Func[FuncNumber],K[FuncNumber][4],Auxiliary[FuncNumber]; for(i=0;i=(reachpoint-x[0])/steplength;i++) { for(j=0;jFuncNumber;j++) Auxiliary[j]=*(Yn[j]+i); rightfunctions(x[i],Auxiliary,Func); for(j=0;jFuncNumber;j++) { K[j][0]=Func[j]; Auxiliary[j]=*(Yn[j]+i)+0.5*steplength*K[j][0]; } rightfunctions(x[i],Auxiliary,Func); for(j=0;jFuncNumber;j++) { K[j][1]=Func[j]; Auxiliary[j]=*(Yn[j]+i)+0.5*steplength*K[j][1]; } rightfunctions(x[i],Auxiliary,Func); for(j=0;jFuncNumber;j++) { K[j][2]=Func[j]; Auxiliary[j]=*(Yn[j]+i)+steplength*K[j][2]; } rightfunctions(x[i],Auxiliary,Func); for(j=0;j