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

/***************************************************************************Thisprogramistosolvetheinitialvalueproblemoffollowingsystemofdifferentialequations:dx/dt=x+2*y,x(0)=0,dy/dt=2*x+y,y(0)=2,x(0.2)andy(0.2)aretobecalculated****************************************************************************/#includestdio.h#includemath.h#definesteplength0.1//步?长¡è可¨¦根¨´据Y需¨¨要°a调Ì¡Â整?;#defineFuncNumber2//FuncNumber为a微¡é分¤?方¤?程¨¬的Ì?数ºy目?;voidmain(){floatx[200],Yn[20][200],reachpoint;inti;x[0]=0;Yn[0][0]=0;Yn[1][0]=2;//初?值¦Ì条¬?件t;reachpoint=0.2;//所¨´求¨®点Ì?可¨¦根¨´据Y需¨¨要°a调Ì¡Â整?;voidrightfunctions(floatx,float*Auxiliary,float*Func);voidRunge_Kutta(float*x,floatreachpoint,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();}voidrightfunctions(floatx,float*Auxiliary,float*Func)//当Ì¡À右®¨°方¤?程¨¬改?变À?时º¡À，ê?需¨¨要°a改?变À?;{Func[0]=Auxiliary[0]+2*Auxiliary[1];Func[1]=2*Auxiliary[0]+Auxiliary[1];}voidRunge_Kutta(float*x,floatreachpoint,float(*Yn)[200]){inti,j;floatFunc[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;jFuncNumber;j++)K[j][3]=Func[j];for(j=0;jFuncNumber;j++)Yn[j][i+1]=Yn[j][i]+(K[j][0]+2*K[j][1]+2*K[j][2]+K[j][3])*steplength/6.0;x[i+1]=x[i]+steplength;}}