【C/C++】龙格库塔+亚当姆斯求解数值微分初值问题

/*
    解数值微分初值问题： 
    龙格-库塔法求前k个初值 + 亚当姆斯法 
*/
#include<bits/stdc++.h>
using namespace std;

double f(double x,double y){
    //y(0) = 1 
    return (y - 2*x/y);
}
void getRungeResult(double *Runge_k,double x0,double y0,double h,int N){
    //求解N个初值，保存在Runge_k[1 to N]中
    double K1,K2,K3,K4;
    double x1,y1;
    for(int i = 1;i<=N;i++){
        x1 = x0+h;
        K1 = f(x0,y0);
        K2 = f(x0+h/2,y0+h/2*K1);
        K3 = f(x0+h/2,y0+h/2*K2);
        K4 = f(x1,y0+K4);
        y1 = y0 + h/6*(K1+2*K2+2*K3+K4);
        Runge_k[i] = y1;
        x0 = x1;
        y0 = y1;
    }
    return;
}

//亚当姆斯多步法 
void Adams(double *Runge_k,double *predict,double x0,double y0,double h,int N){
    Runge_k[0] = y0;
    //（0）龙格库塔法求前4个初值 
    getRungeResult(Runge_k,x0,y0,h,3);
    double y1,y2,y3,dy0,dy1,dy2,dy3;
    y1 = Runge_k[1];
    y2 = Runge_k[2];
    y3 = Runge_k[3];
    dy0 = f(x0,y0);
    dy1 = f(x0+h,y1);
    dy2 = f(x0+2*h,y2);
    dy3 = f(x0+3*h,y3);
    double x3 = x0+3*h;
    double x4,y4,yp,dyp,dy4;
    for(int i = 4;i<=N;i++){
        x4 = x3+h;
        //(1)预测 
        yp = y3 + h/24*(55*dy3-59*dy2+37*dy1-9*dy0);
        predict[i] = yp;//保存预测值 
        //预测要用dyp 
        dyp = f(x4,yp);
        //(2)校正 
        y4 = y3 + h/24*(9*dyp + 19*dy3 -5*dy2+dy1);
        //存起来
        Runge_k[i] = y4;
        //求下一次需要用到导
        dy4 = f(x4,y4);
        //为下一次循环做准备 
        x3 = x4;
        y3 = y4;
        dy0 = dy1;
        dy1 = dy2;
        dy2 = dy3;
        dy3 = dy4;
    }
    return;
} 


/*假设这里保证四阶精度*/
int main(){
    /*说明：x0,y0是初值，h是小区间长度，N是要求的个数*/
    double x0,y0,h;
    int N;
    cout<<"输入初值x0,y0,小区间h,需要的初值个数N：";
    cin>>x0>>y0>>h>>N;
    //保存Runge求的4个初始值,龙格法求3个就可以;之后也用这个保存最终的Adams结果 
    double Runge_k[100];
    //保存预测值，方便以后比较 
    double predict[100];
    memset(predict,0,sizeof(predict));
    memset(Runge_k,0,sizeof(Runge_k));
    Adams(Runge_k,predict,x0,y0,h,N);
    cout<<endl;
    printf("预测值："); 
    for(int i = 0;i<=N;i++){
        if(i<4){
            printf("%.6lf ",0);
        }else{
            printf("%.6lf ",predict[i]);
        }
    }
    cout<<endl;
    printf("校正值：");
    for(int i = 0;i<=N;i++){
        printf("%.6lf ",Runge_k[i]);
    }
    
}