龙格库塔算法解微分方程组 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;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],Auxiliary,Func);

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

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

    }

}