问题简化为三个质点的质量相同的情况

输入所求微分方程组的初值 t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3

输入所求微分方程组的微分区间[a,b]

输入所求微分方程组所分解子区间的个数step

输出文件可以直接用Matlab来作图

当时怎么写的我已经完全忘记了，我只记得当时很佩服自己

以上是推导过程，然后给出代码实现：

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 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:"<
         
          >m;220     m=1;221     //cout<<"输入所求微分方程组的初值 t0,x1,y1,x2,y2,x3,y3,vx1,vy1,vx2,vy2,vx3,vy3:"<
          
           >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]:"<
           
            >a>>b;236 237 //cout<<"输入所求微分方程组所分解子区间的个数step:"<
            
             >step;239 240 cout<
             
              <
              
               <
              
             
            
           
          
         
        
       
      
     

数值计算算是一个比较有意思的方向