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

/***************************************************************************

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

****************************************************************************/

#include

#include

#define steplength 0.1 //步?长¡è可¨¦根¨´据Y需¨¨要°a调Ì¡Â整?;

#define FuncNumber 2 //FuncNumber为a微¡é分¤?方¤?程¨¬的Ì?数ºy目?;

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)//当Ì¡À右®¨°方¤?程¨¬改?变À?时º¡À，ê?需¨¨要°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;jAuxiliary[j]=*(Yn[j]+i);

rightfunctions(x[i],Auxiliary,Func);

for(j=0;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;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;j{

K[j][2]=Func[j];

Auxiliary[j]=*(Yn[j]+i)+steplength*K[j][2];

}

rightfunctions(x[i],Auxiliary,Func);

for(j=0;jK[j][3]=Func[j];

for(j=0;jYn[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;

}

}