使用四阶龙格库塔方法求解三体问题（解十二元一阶常微分方程组）
问题简化为三个质点的质量相同的情况
输入所求微分方程组的初值 t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3
输入所求微分方程组的微分区间[a,b]
输入所求微分方程组所分解子区间的个数step
输出文件可以直接用Matlab来作图
当时怎么写的我已经完全忘记了，我只记得当时很佩服自己
以上是推导过程，然后给出代码实现：
1 #include<iostream> 2 #include<iomanip> 3 #include<cmath> 4 #include<cstdio> 5 const double G=1; //万有引力常量 6 double m; //三个质点的同一质量 7 double initial[20],result[20]; //迭代数组 8 double a,b,h; //区间和步长 9 double step; //分段数 10 using namespace std; 11 12 //以下六个函数对应六个位移对时间求导的微分方程 13 double fx1(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3) 14 { 15 double dx1; 16 dx1=vx1; 17 return dx1; 18 } 19 double gy1(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3) 20 { 21 double dy1; 22 dy1=vy1; 23 return dy1; 24 } 25 double fx2(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3) 26 { 27 double dx2; 28 dx2=vx2; 29 return dx2; 30 } 31 double gy2(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3) 32 { 33 double dy2; 34 dy2=vy2; 35 return dy2; 36 } 37 double fx3(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3) 38 { 39 double dx3; 40 dx3=vx3; 41 return dx3; 42 } 43 double gy3(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3) 44 { 45 double dy3; 46 dy3=vy3; 47 return dy3; 48 } 49 50 //以下六个函数对应速度对时间求导的微分方程 51 double fvx1(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3) 52 { 53 double dvx1; 54 dvx1=G*m*pow((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1),-1.5)*(x2-x1)+G*m*pow((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1),-1.5)*(x3-x1); 55 return dvx1; 56 } 57 double gvy1(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3) 58 { 59 double dvy1; 60 dvy1=G*m*pow((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1),-1.5)*(y2-y1)+G*m*pow((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1),-1.5)*(y3-y1); 61 return dvy1; 62 } 63 double fvx2(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3) 64 { 65 double dvx2; 66 dvx2=G*m*pow((y1-y2)*(y1-y2)+(x1-x2)*(x1-x2),-1.5)*(x1-x2)+G*m*pow((y3-y2)*(y3-y2)+(x3-x2)*(x3-x2),-1.5)*(x3-x2); 67 return dvx2; 68 } 69 double gvy2(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3) 70 { 71 double dvy2; 72 dvy2=G*m*pow((y1-y2)*(y1-y2)+(x1-x2)*(x1-x2),-1.5)*(y1-y2)+G*m*pow((y3-y2)*(y3-y2)+(x3-x2)*(x3-x2),-1.5)*(y3-y2); 73 return dvy2; 74 } 75 double fvx3(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3) 76 { 77 double dvx3; 78 dvx3=G*m*pow((y2-y3)*(y2-y3)+(x2-x3)*(x2-x3),-1.5)*(x2-x3)+G*m*pow((y1-y3)*(y1-y3)+(x1-x3)*(x1-x3),-1.5)*(x1-x3); 79 return dvx3; 80 } 81 double gvy3(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3) 82 { 83 double dvy3; 84 dvy3=G*m*pow((y2-y3)*(y2-y3)+(x2-x3)*(x2-x3),-1.5)*(y2-y3)+G*m*pow((y1-y3)*(y1-y3)+(x1-x3)*(x1-x3),-1.5)*(y1-y3); 85 return dvy3; 86 } 87 //龙格库塔方法 88 void RK4( 89 double (*fx1)(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3), 90 double (*gy1)(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3), 91 double (*fx2)(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3), 92 double (*gy2)(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3), 93 double (*fx3)(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3), 94 double (*gy3)(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3), 95 96 double (*fvx1)(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3), 97 double (*gvy1)(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3), 98 double (*fvx2)(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3), 99 double (*gvy2)(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3), 100 double (*fvx3)(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3), 101 double (*gvy3)(double t,double x1,double y1,double x2,double y2,double x3,double y3,double vx1,double vy1,double vx2,double vy2,double vx3,double vy3) 102 ) 103 { 104 //前六个方程的中间值 105 double x1f1,x1f2,x1f3,x1f4,y1g1,y1g2,y1g3,y1g4; 106 double x2f1,x2f2,x2f3,x2f4,y2g1,y2g2,y2g3,y2g4; 107 double x3f1,x3f2,x3f3,x3f4,y3g1,y3g2,y3g3,y3g4; 108 //后六个方程的中间值 109 double vx1f1,vx1f2,vx1f3,vx1f4,vy1g1,vy1g2,vy1g3,vy1g4; 110 double vx2f1,vx2f2,vx2f3,vx2f4,vy2g1,vy2g2,vy2g3,vy2g4; 111 double vx3f1,vx3f2,vx3f3,vx3f4,vy3g1,vy3g2,vy3g3,vy3g4; 112 //迭代值 113 double t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3; 114 double nx1,ny1,nx2,ny2,nx3,ny3,nvx1,nvy1,nvx2,nvy2,nvx3,nvy3; 115 //初始化 116 t0=initial[0]; 117 x1=initial[1]; 118 y1=initial[2]; 119 x2=initial[3]; 120 y2=initial[4]; 121 x3=initial[5]; 122 y3=initial[6]; 123 vx1=initial[7]; 124 vy1=initial[8]; 125 vx2=initial[9]; 126 vy2=initial[10]; 127 vx3=initial[11]; 128 vy3=initial[12]; 129 //计算k1 130 x1f1=fx1(t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3); 131 y1g1=gy1(t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3); 132 x2f1=fx2(t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3); 133 y2g1=gy2(t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3); 134 x3f1=fx3(t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3); 135 y3g1=gy3(t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3); 136 137 vx1f1=fvx1(t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3); 138 vy1g1=gvy1(t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3); 139 vx2f1=fvx2(t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3); 140 vy2g1=gvy2(t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3); 141 vx3f1=fvx3(t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3); 142 vy3g1=gvy3(t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3); 143 //计算k2 144 x1f2=fx1(t0+h/2,x1+h*x1f1/2,y1+h*y1g1/2,x2+h*x2f1/2,y2+h*y2g1/2,x3+h*x3f1/2,y3+h*y3g1/2,vx1+h*vx1f1/2,vy1+h*vy1g1/2,vx2+h*vx2f1/2,vy2+h*vy2g1/2,vx3+h*vx3f1/2,vy3+h*vy3g1/2); 145 y1g2=gy1(t0+h/2,x1+h*x1f1/2,y1+h*y1g1/2,x2+h*x2f1/2,y2+h*y2g1/2,x3+h*x3f1/2,y3+h*y3g1/2,vx1+h*vx1f1/2,vy1+h*vy1g1/2,vx2+h*vx2f1/2,vy2+h*vy2g1/2,vx3+h*vx3f1/2,vy3+h*vy3g1/2); 146 x2f2=fx2(t0+h/2,x1+h*x1f1/2,y1+h*y1g1/2,x2+h*x2f1/2,y2+h*y2g1/2,x3+h*x3f1/2,y3+h*y3g1/2,vx1+h*vx1f1/2,vy1+h*vy1g1/2,vx2+h*vx2f1/2,vy2+h*vy2g1/2,vx3+h*vx3f1/2,vy3+h*vy3g1/2); 147 y2g2=gy2(t0+h/2,x1+h*x1f1/2,y1+h*y1g1/2,x2+h*x2f1/2,y2+h*y2g1/2,x3+h*x3f1/2,y3+h*y3g1/2,vx1+h*vx1f1/2,vy1+h*vy1g1/2,vx2+h*vx2f1/2,vy2+h*vy2g1/2,vx3+h*vx3f1/2,vy3+h*vy3g1/2); 148 x3f2=fx3(t0+h/2,x1+h*x1f1/2,y1+h*y1g1/2,x2+h*x2f1/2,y2+h*y2g1/2,x3+h*x3f1/2,y3+h*y3g1/2,vx1+h*vx1f1/2,vy1+h*vy1g1/2,vx2+h*vx2f1/2,vy2+h*vy2g1/2,vx3+h*vx3f1/2,vy3+h*vy3g1/2); 149 y3g2=gy3(t0+h/2,x1+h*x1f1/2,y1+h*y1g1/2,x2+h*x2f1/2,y2+h*y2g1/2,x3+h*x3f1/2,y3+h*y3g1/2,vx1+h*vx1f1/2,vy1+h*vy1g1/2,vx2+h*vx2f1/2,vy2+h*vy2g1/2,vx3+h*vx3f1/2,vy3+h*vy3g1/2); 150 151 vx1f2=fvx1(t0+h/2,x1+h*x1f1/2,y1+h*y1g1/2,x2+h*x2f1/2,y2+h*y2g1/2,x3+h*x3f1/2,y3+h*y3g1/2,vx1+h*vx1f1/2,vy1+h*vy1g1/2,vx2+h*vx2f1/2,vy2+h*vy2g1/2,vx3+h*vx3f1/2,vy3+h*vy3g1/2); 152 vy1g2=gvy1(t0+h/2,x1+h*x1f1/2,y1+h*y1g1/2,x2+h*x2f1/2,y2+h*y2g1/2,x3+h*x3f1/2,y3+h*y3g1/2,vx1+h*vx1f1/2,vy1+h*vy1g1/2,vx2+h*vx2f1/2,vy2+h*vy2g1/2,vx3+h*vx3f1/2,vy3+h*vy3g1/2); 153 vx2f2=fvx2(t0+h/2,x1+h*x1f1/2,y1+h*y1g1/2,x2+h*x2f1/2,y2+h*y2g1/2,x3+h*x3f1/2,y3+h*y3g1/2,vx1+h*vx1f1/2,vy1+h*vy1g1/2,vx2+h*vx2f1/2,vy2+h*vy2g1/2,vx3+h*vx3f1/2,vy3+h*vy3g1/2); 154 vy2g2=gvy2(t0+h/2,x1+h*x1f1/2,y1+h*y1g1/2,x2+h*x2f1/2,y2+h*y2g1/2,x3+h*x3f1/2,y3+h*y3g1/2,vx1+h*vx1f1/2,vy1+h*vy1g1/2,vx2+h*vx2f1/2,vy2+h*vy2g1/2,vx3+h*vx3f1/2,vy3+h*vy3g1/2); 155 vx3f2=fvx3(t0+h/2,x1+h*x1f1/2,y1+h*y1g1/2,x2+h*x2f1/2,y2+h*y2g1/2,x3+h*x3f1/2,y3+h*y3g1/2,vx1+h*vx1f1/2,vy1+h*vy1g1/2,vx2+h*vx2f1/2,vy2+h*vy2g1/2,vx3+h*vx3f1/2,vy3+h*vy3g1/2); 156 vy3g2=gvy3(t0+h/2,x1+h*x1f1/2,y1+h*y1g1/2,x2+h*x2f1/2,y2+h*y2g1/2,x3+h*x3f1/2,y3+h*y3g1/2,vx1+h*vx1f1/2,vy1+h*vy1g1/2,vx2+h*vx2f1/2,vy2+h*vy2g1/2,vx3+h*vx3f1/2,vy3+h*vy3g1/2); 157 //计算k3 158 x1f3=fx1(t0+h/2,x1+h*x1f2/2,y1+h*y1g2/2,x2+h*x2f2/2,y2+h*y2g2/2,x3+h*x3f2/2,y3+h*y3g2/2,vx1+h*vx1f2/2,vy1+h*vy1g2/2,vx2+h*vx2f2/2,vy2+h*vy2g2/2,vx3+h*vx3f2/2,vy3+h*vy3g2/2); 159 y1g3=gy1(t0+h/2,x1+h*x1f2/2,y1+h*y1g2/2,x2+h*x2f2/2,y2+h*y2g2/2,x3+h*x3f2/2,y3+h*y3g2/2,vx1+h*vx1f2/2,vy1+h*vy1g2/2,vx2+h*vx2f2/2,vy2+h*vy2g2/2,vx3+h*vx3f2/2,vy3+h*vy3g2/2); 160 x2f3=fx2(t0+h/2,x1+h*x1f2/2,y1+h*y1g2/2,x2+h*x2f2/2,y2+h*y2g2/2,x3+h*x3f2/2,y3+h*y3g2/2,vx1+h*vx1f2/2,vy1+h*vy1g2/2,vx2+h*vx2f2/2,vy2+h*vy2g2/2,vx3+h*vx3f2/2,vy3+h*vy3g2/2); 161 y2g3=gy2(t0+h/2,x1+h*x1f2/2,y1+h*y1g2/2,x2+h*x2f2/2,y2+h*y2g2/2,x3+h*x3f2/2,y3+h*y3g2/2,vx1+h*vx1f2/2,vy1+h*vy1g2/2,vx2+h*vx2f2/2,vy2+h*vy2g2/2,vx3+h*vx3f2/2,vy3+h*vy3g2/2); 162 x3f3=fx3(t0+h/2,x1+h*x1f2/2,y1+h*y1g2/2,x2+h*x2f2/2,y2+h*y2g2/2,x3+h*x3f2/2,y3+h*y3g2/2,vx1+h*vx1f2/2,vy1+h*vy1g2/2,vx2+h*vx2f2/2,vy2+h*vy2g2/2,vx3+h*vx3f2/2,vy3+h*vy3g2/2); 163 y3g3=gy3(t0+h/2,x1+h*x1f2/2,y1+h*y1g2/2,x2+h*x2f2/2,y2+h*y2g2/2,x3+h*x3f2/2,y3+h*y3g2/2,vx1+h*vx1f2/2,vy1+h*vy1g2/2,vx2+h*vx2f2/2,vy2+h*vy2g2/2,vx3+h*vx3f2/2,vy3+h*vy3g2/2); 164 165 vx1f3=fvx1(t0+h/2,x1+h*x1f2/2,y1+h*y1g2/2,x2+h*x2f2/2,y2+h*y2g2/2,x3+h*x3f2/2,y3+h*y3g2/2,vx1+h*vx1f2/2,vy1+h*vy1g2/2,vx2+h*vx2f2/2,vy2+h*vy2g2/2,vx3+h*vx3f2/2,vy3+h*vy3g2/2); 166 vy1g3=gvy1(t0+h/2,x1+h*x1f2/2,y1+h*y1g2/2,x2+h*x2f2/2,y2+h*y2g2/2,x3+h*x3f2/2,y3+h*y3g2/2,vx1+h*vx1f2/2,vy1+h*vy1g2/2,vx2+h*vx2f2/2,vy2+h*vy2g2/2,vx3+h*vx3f2/2,vy3+h*vy3g2/2); 167 vx2f3=fvx2(t0+h/2,x1+h*x1f2/2,y1+h*y1g2/2,x2+h*x2f2/2,y2+h*y2g2/2,x3+h*x3f2/2,y3+h*y3g2/2,vx1+h*vx1f2/2,vy1+h*vy1g2/2,vx2+h*vx2f2/2,vy2+h*vy2g2/2,vx3+h*vx3f2/2,vy3+h*vy3g2/2); 168 vy2g3=gvy2(t0+h/2,x1+h*x1f2/2,y1+h*y1g2/2,x2+h*x2f2/2,y2+h*y2g2/2,x3+h*x3f2/2,y3+h*y3g2/2,vx1+h*vx1f2/2,vy1+h*vy1g2/2,vx2+h*vx2f2/2,vy2+h*vy2g2/2,vx3+h*vx3f2/2,vy3+h*vy3g2/2); 169 vx3f3=fvx3(t0+h/2,x1+h*x1f2/2,y1+h*y1g2/2,x2+h*x2f2/2,y2+h*y2g2/2,x3+h*x3f2/2,y3+h*y3g2/2,vx1+h*vx1f2/2,vy1+h*vy1g2/2,vx2+h*vx2f2/2,vy2+h*vy2g2/2,vx3+h*vx3f2/2,vy3+h*vy3g2/2); 170 vy3g3=gvy3(t0+h/2,x1+h*x1f2/2,y1+h*y1g2/2,x2+h*x2f2/2,y2+h*y2g2/2,x3+h*x3f2/2,y3+h*y3g2/2,vx1+h*vx1f2/2,vy1+h*vy1g2/2,vx2+h*vx2f2/2,vy2+h*vy2g2/2,vx3+h*vx3f2/2,vy3+h*vy3g2/2); 171 //计算k4 172 x1f4=fx1(t0+h,x1+h*x1f3,y1+h*y1g3,x2+h*x2f3,y2+h*y2g3,x3+h*x3f3,y3+h*y3g3,vx1+h*vx1f3,vy1+h*vy1g3,vx2+h*vx2f3,vy2+h*vy2g3,vx3+h*vx3f3,vy3+h*vy3g3); 173 y1g4=gy1(t0+h,x1+h*x1f3,y1+h*y1g3,x2+h*x2f3,y2+h*y2g3,x3+h*x3f3,y3+h*y3g3,vx1+h*vx1f3,vy1+h*vy1g3,vx2+h*vx2f3,vy2+h*vy2g3,vx3+h*vx3f3,vy3+h*vy3g3); 174 x2f4=fx2(t0+h,x1+h*x1f3,y1+h*y1g3,x2+h*x2f3,y2+h*y2g3,x3+h*x3f3,y3+h*y3g3,vx1+h*vx1f3,vy1+h*vy1g3,vx2+h*vx2f3,vy2+h*vy2g3,vx3+h*vx3f3,vy3+h*vy3g3); 175 y2g4=gy2(t0+h,x1+h*x1f3,y1+h*y1g3,x2+h*x2f3,y2+h*y2g3,x3+h*x3f3,y3+h*y3g3,vx1+h*vx1f3,vy1+h*vy1g3,vx2+h*vx2f3,vy2+h*vy2g3,vx3+h*vx3f3,vy3+h*vy3g3); 176 x3f4=fx3(t0+h,x1+h*x1f3,y1+h*y1g3,x2+h*x2f3,y2+h*y2g3,x3+h*x3f3,y3+h*y3g3,vx1+h*vx1f3,vy1+h*vy1g3,vx2+h*vx2f3,vy2+h*vy2g3,vx3+h*vx3f3,vy3+h*vy3g3); 177 y3g4=gy3(t0+h,x1+h*x1f3,y1+h*y1g3,x2+h*x2f3,y2+h*y2g3,x3+h*x3f3,y3+h*y3g3,vx1+h*vx1f3,vy1+h*vy1g3,vx2+h*vx2f3,vy2+h*vy2g3,vx3+h*vx3f3,vy3+h*vy3g3); 178 179 vx1f4=fvx1(t0+h,x1+h*x1f3,y1+h*y1g3,x2+h*x2f3,y2+h*y2g3,x3+h*x3f3,y3+h*y3g3,vx1+h*vx1f3,vy1+h*vy1g3,vx2+h*vx2f3,vy2+h*vy2g3,vx3+h*vx3f3,vy3+h*vy3g3); 180 vy1g4=gvy1(t0+h,x1+h*x1f3,y1+h*y1g3,x2+h*x2f3,y2+h*y2g3,x3+h*x3f3,y3+h*y3g3,vx1+h*vx1f3,vy1+h*vy1g3,vx2+h*vx2f3,vy2+h*vy2g3,vx3+h*vx3f3,vy3+h*vy3g3); 181 vx2f4=fvx2(t0+h,x1+h*x1f3,y1+h*y1g3,x2+h*x2f3,y2+h*y2g3,x3+h*x3f3,y3+h*y3g3,vx1+h*vx1f3,vy1+h*vy1g3,vx2+h*vx2f3,vy2+h*vy2g3,vx3+h*vx3f3,vy3+h*vy3g3); 182 vy2g4=gvy2(t0+h,x1+h*x1f3,y1+h*y1g3,x2+h*x2f3,y2+h*y2g3,x3+h*x3f3,y3+h*y3g3,vx1+h*vx1f3,vy1+h*vy1g3,vx2+h*vx2f3,vy2+h*vy2g3,vx3+h*vx3f3,vy3+h*vy3g3); 183 vx3f4=fvx3(t0+h,x1+h*x1f3,y1+h*y1g3,x2+h*x2f3,y2+h*y2g3,x3+h*x3f3,y3+h*y3g3,vx1+h*vx1f3,vy1+h*vy1g3,vx2+h*vx2f3,vy2+h*vy2g3,vx3+h*vx3f3,vy3+h*vy3g3); 184 vy3g4=gvy3(t0+h,x1+h*x1f3,y1+h*y1g3,x2+h*x2f3,y2+h*y2g3,x3+h*x3f3,y3+h*y3g3,vx1+h*vx1f3,vy1+h*vy1g3,vx2+h*vx2f3,vy2+h*vy2g3,vx3+h*vx3f3,vy3+h*vy3g3); 185 //计算最终结果 186 nx1=x1+h*(x1f1+2*x1f2+2*x1f3+x1f4)/6; 187 ny1=y1+h*(y1g1+2*y1g2+2*y1g3+y1g4)/6; 188 nx2=x2+h*(x2f1+2*x2f2+2*x2f3+x2f4)/6; 189 ny2=y2+h*(y2g1+2*y2g2+2*y2g3+y2g4)/6; 190 nx3=x3+h*(x3f1+2*x3f2+2*x3f3+x3f4)/6; 191 ny3=y3+h*(y3g1+2*y3g2+2*y3g3+y3g4)/6; 192 193 nvx1=vx1+h*(vx1f1+2*vx1f2+2*vx1f3+vx1f4)/6; 194 nvy1=vy1+h*(vy1g1+2*vy1g2+2*vy1g3+vy1g4)/6; 195 nvx2=vx2+h*(vx2f1+2*vx2f2+2*vx2f3+vx2f4)/6; 196 nvy2=vy2+h*(vy2g1+2*vy2g2+2*vy2g3+vy2g4)/6; 197 nvx3=vx3+h*(vx3f1+2*vx3f2+2*vx3f3+vx3f4)/6; 198 nvy3=vy3+h*(vy3g1+2*vy3g2+2*vy3g3+vy3g4)/6; 199 //迭代过程 200 result[0]=t0+h; 201 result[1]=nx1; 202 result[2]=ny1; 203 result[3]=nx2; 204 result[4]=ny2; 205 result[5]=nx3; 206 result[6]=ny3; 207 result[7]=nvx1; 208 result[8]=nvy1; 209 result[9]=nvx2; 210 result[10]=nvy2; 211 result[11]=nvx3; 212 result[12]=nvy3; 213 } 214 int main() 215 { 216 freopen("input.txt","r",stdin); 217 freopen("ouput.txt","w",stdout); 218 //cout<<"输入三个质点的相同质量m:"<<endl; 219 //cin>>m; 220 m=1; 221 //cout<<"输入所求微分方程组的初值 t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3:"<<endl; 222 initial[0]=0; 223 initial[1]=-1; 224 initial[2]=0; 225 initial[3]=1; 226 initial[4]=0; 227 initial[5]=0; 228 initial[6]=0; 229 cin>>initial[7]>>initial[8]; 230 initial[9]=initial[7]; 231 initial[10]=initial[8]; 232 initial[11]=-2*initial[9]; 233 initial[12]=-2*initial[10]; 234 //cout<<"输入所求微分方程组的微分区间[a,b]:"<<endl; 235 cin>>a>>b; 236 237 //cout<<"输入所求微分方程组所分解子区间的个数step:"<<endl; 238 cin>>step; 239 240 cout<<setiosflags(ios::right)<<setiosflags(ios::fixed)<<setprecision(10); 241 cout<<setw(15)<<"t"<<setw(15)<<"x1"<<setw(15)<<"y1"<<setw(15)<<"x2"<<setw(15)<<"y2"<<setw(15)<<"x3"<<setw(15)<<"y3"; 242 cout<<setw(15)<<"Vx1"<<setw(15)<<"Vx2"<<setw(15)<<"Vy1"<<setw(15)<<"Vy2"<<setw(15)<<"Vx3"<<setw(15)<<"Vy3"<<endl; 243 h=(b-a)/step; //计算步长 244 for(int i=0;i<=12;i++) 245 cout<<setw(15)<<initial[i]; //输出初始值 246 cout<<endl; 247 for(int i=0;i<step;i++) 248 { 249 RK4(fx1,gy1,fx2,gy2,fx3,gy3,fvx1,gvy1,fvx2,gvy2,fvx3,gvy3); 250 for(int i=0;i<=12;i++) 251 cout<<setw(15)<<result[i]; //输出每一次的迭代结果 252 cout<<endl; 253 254 for(int i=0;i<=12;i++) 255 initial[i]=result[i]; //迭代 256 } 257 return 0; 258 }
数值计算算是一个比较有意思的方向
