1 /*
 2     解数值微分初值问题：
 3     龙格-库塔法求前k个初值 + 亚当姆斯法
 4 */
 5 #include<bits/stdc++.h>
 6 using namespace std;
 7
 8 double f(double x,double y){
 9     //y(0) = 1
10     return (y - 2*x/y);
11 }
12 void getRungeResult(double *Runge_k,double x0,double y0,double h,int N){
13     //求解N个初值，保存在Runge_k[1 to N]中
14     double K1,K2,K3,K4;
15     double x1,y1;
16     for(int i = 1;i<=N;i++){
17         x1 = x0+h;
18         K1 = f(x0,y0);
19         K2 = f(x0+h/2,y0+h/2*K1);
20         K3 = f(x0+h/2,y0+h/2*K2);
21         K4 = f(x1,y0+K4);
22         y1 = y0 + h/6*(K1+2*K2+2*K3+K4);
23         Runge_k[i] = y1;
24         x0 = x1;
25         y0 = y1;
26     }
27     return;
28 }
29
30 //亚当姆斯多步法
31 void Adams(double *Runge_k,double *predict,double x0,double y0,double h,int N){
32     Runge_k[0] = y0;
33     //（0）龙格库塔法求前4个初值
34     getRungeResult(Runge_k,x0,y0,h,3);
35     double y1,y2,y3,dy0,dy1,dy2,dy3;
36     y1 = Runge_k[1];
37     y2 = Runge_k[2];
38     y3 = Runge_k[3];
39     dy0 = f(x0,y0);
40     dy1 = f(x0+h,y1);
41     dy2 = f(x0+2*h,y2);
42     dy3 = f(x0+3*h,y3);
43     double x3 = x0+3*h;
44     double x4,y4,yp,dyp,dy4;
45     for(int i = 4;i<=N;i++){
46         x4 = x3+h;
47         //(1)预测
48         yp = y3 + h/24*(55*dy3-59*dy2+37*dy1-9*dy0);
49         predict[i] = yp;//保存预测值
50         //预测要用dyp
51         dyp = f(x4,yp);
52         //(2)校正
53         y4 = y3 + h/24*(9*dyp + 19*dy3 -5*dy2+dy1);
54         //存起来
55         Runge_k[i] = y4;
56         //求下一次需要用到导
57         dy4 = f(x4,y4);
58         //为下一次循环做准备
59         x3 = x4;
60         y3 = y4;
61         dy0 = dy1;
62         dy1 = dy2;
63         dy2 = dy3;
64         dy3 = dy4;
65     }
66     return;
67 }
68
69
70 /*假设这里保证四阶精度*/
71 int main(){
72     /*说明：x0,y0是初值，h是小区间长度，N是要求的个数*/
73     double x0,y0,h;
74     int N;
75     cout<<"输入初值x0,y0,小区间h,需要的初值个数N：";
76     cin>>x0>>y0>>h>>N;
77     //保存Runge求的4个初始值,龙格法求3个就可以;之后也用这个保存最终的Adams结果
78     double Runge_k[100];
79     //保存预测值，方便以后比较
80     double predict[100];
81     memset(predict,0,sizeof(predict));
82     memset(Runge_k,0,sizeof(Runge_k));
83     Adams(Runge_k,predict,x0,y0,h,N);
84     cout<<endl;
85     printf("预测值：");
86     for(int i = 0;i<=N;i++){
87         if(i<4){
88             printf("%.6lf ",0);
89         }else{
90             printf("%.6lf ",predict[i]);
91         }
92     }
93     cout<<endl;
94     printf("校正值：");
95     for(int i = 0;i<=N;i++){
96         printf("%.6lf ",Runge_k[i]);
97     }
98
99 }