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 }
轉載于:https://www.cnblogs.com/duye/p/9172043.html