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

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

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微?é分¤?方¤?程¨?的ì?数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;j

Auxiliary[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],Auxil