matlab ode45三体问题,使用四阶龙格库塔方法求解三体问题（解十二元一阶常微分方程组）...

1 #include

2 #include

3 #include

4 #include

5 const double G=1; //万有引力常量

6 double m; //三个质点的同一质量

7 double initial[20],result[20]; //迭代数组

8 double a,b,h; //区间和步长

9 double step; //分段数

10 using namespacestd;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,doublevy3)14 {15 doubledx1;16 dx1=vx1;17 returndx1;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,doublevy3)20 {21 doubledy1;22 dy1=vy1;23 returndy1;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,doublevy3)26 {27 doubledx2;28 dx2=vx2;29 returndx2;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,doublevy3)32 {33 doubledy2;34 dy2=vy2;35 returndy2;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,doublevy3)38 {39 doubledx3;40 dx3=vx3;41 returndx3;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,doublevy3)44 {45 doubledy3;46 dy3=vy3;47 returndy3;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,doublevy3)52 {53 doubledvx1;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 returndvx1;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,doublevy3)58 {59 doubledvy1;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 returndvy1;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,doublevy3)64 {65 doubledvx2;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 returndvx2;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,doublevy3)70 {71 doubledvy2;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 returndvy2;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,doublevy3)76 {77 doubledvx3;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 returndvx3;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,doublevy3)82 {83 doubledvy3;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 returndvy3;86 }87 //龙格库塔方法

88 voidRK4(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,doublevy3),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,doublevy3),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,doublevy3),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,doublevy3),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,doublevy3),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,doublevy3),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,doublevy3),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,doublevy3),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,doublevy3),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,doublevy3),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,doublevy3),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,doublevy3)102 )103 {104 //前六个方程的中间值

105 doublex1f1,x1f2,x1f3,x1f4,y1g1,y1g2,y1g3,y1g4;106 doublex2f1,x2f2,x2f3,x2f4,y2g1,y2g2,y2g3,y2g4;107 doublex3f1,x3f2,x3f3,x3f4,y3g1,y3g2,y3g3,y3g4;108 //后六个方程的中间值

109 doublevx1f1,vx1f2,vx1f3,vx1f4,vy1g1,vy1g2,vy1g3,vy1g4;110 doublevx2f1,vx2f2,vx2f3,vx2f4,vy2g1,vy2g2,vy2g3,vy2g4;111 doublevx3f1,vx3f2,vx3f3,vx3f4,vy3g1,vy3g2,vy3g3,vy3g4;112 //迭代值

113 doublet0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3;114 doublenx1,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 intmain()215 {216 freopen("input.txt","r",stdin);217 freopen("ouput.txt","w",stdout);218 //cout<>m;

220 m=1;221 //cout<

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<

235 cin>>a>>b;236

237 //cout<

238 cin>>step;239

240 cout<

244 for(int i=0;i<=12;i++)245 cout<

246 cout<

252 cout<

254 for(int i=0;i<=12;i++)255 initial[i]=result[i]; //迭代

256 }257 return 0;258 }