/***************************************************************************
This program is to solve the initial value problem of following system
of differential equati ons:
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<stdi o.h>
#include<math.h>
#define steplength 0.1 //步?长¡è可¨¦根¨´据Y需¨¨要°a调Ì¡Â整?;
#define FuncNumber 2 //FuncNumber为a微¡é分¤?方¤?程¨¬的Ì?数ºy目?;
void main()
{
float x[200],Yn[20][200],reachpoint;i nt 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 R unge_Kutta(float *x,float reachpoint, float(*Yn)[200]);
Runge_Kutta(x ,reachpoi nt, 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 R unge_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<FuncNumber;j++)
Auxiliary[j]=*(Yn[j]+i);
rightfunctions(x[i],Auxiliary,Func);
for(j=0;j<FuncNumber;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<FuncNumber;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<FuncNumber;j++)
{
K[j][2]=Func[j];
Auxiliary[j]=*(Yn[j]+i)+steplength*K[j][2];
}
rightfunctions(x[i],Auxiliary,Func);
for(j=0;j<FuncNumber;j++)
K[j][3]=Func[j];
for(j=0;j<FuncNumber;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;
}
}。