1. 基本算法练习

double f(double x, double y) {

    return  .... ;

}

//evolve df one step size

double Euler_predict_correct(double y, double x, double h){

    double k1 = ....;

    double k2 = ....;

    return y + h*(k1+k2)/2;

}

int main() {

    double h=0.1,x=0, y=1;

    printf("------Euler prediction correction method------\n");

    for(int i=0; i<10; i++){

        y = Euler_predict_correct(y, x, h);     //调用函数，演化一步

        x += h;

        printf("%3d  %.2f   %f  exact=%f\n", i, x, y, sqrt(2*x+1));

    }

    return 0;  

}

龙格库塔方法类似。

2. 磁场中带电粒子的运动

const int N=6;

double omegac=2.0;

//根据t,Y[]计算F[]， 注意Y，F都是数组

void Feval(double t, double Y[], double F[]){

    F[0] = Y[3] ;

    F[1] = Y[4] ;

    F[2] = Y[5] ;

    F[3] = omegac*Y[4] ;

    F[4] =-1*omegac*Y[3] ;

    F[5] = 0;

}

void Runge_Kutta_vector_solve(double t, double dt, double Y[]){

    double k1[N], k2[N], k3[N], k4[N];

    double y2[N], y3[N], y4[N];

    Feval(t, Y, k1);    //计算k1数组

    for(int i=0; i<N; i++){

        y2[i] = Y[i] + k1[i]*dt/2 ;

    }

    Feval(t+dt/2.0, y2, k2);  //计算k2数组

    //k3、k4数组的计算类似k2

    ....

    for(int i=0; i<N; i++){

        Y[i] = Y[i] + dt*(k1[i] + 2*k2[i] + 2*k3[i] + k4[i])/6.0;

    }

}

int main(){

    double Y[N]={1,1,0,2,1,8};    //初始状态

    double dt=0.1, t0=0;

    FILE *fp = fopen("xyz.dat","w");      //打开文件

    for(int i=0; i<100; i++){

        double t =t0 + i*dt;

        Runge_Kutta_vector_solve(t,dt,Y);    //调用程序

        //将（x,y,z)数据写入文件。

        ...

    }

    fclose(fp);

    return 0;

}