我正在上初学者编程课程，我们有一个小项目。我选择使用欧拉方法来模拟三体问题。即使系统是混乱的，也有一些特殊情况会产生对称的轨迹。我得到了大部分，但最重要的一个没有出来。它被称为图 8，应该是稳定的。
我确信初始条件很好。程序也应该如此，因为我得到了其他预期的结果。所以问题是我的程序是否有错误，或者欧拉方法不适用于这种特定情况。（然后我会更改为 Runge-Kutta 方法）
我很感激任何帮助！
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;
double const ms=1.0000000000000000;
double const me=1.0000000000000000;
double const mm=1.0000000000000000;
double const G=1.000000000000000000;
double const dt=0.0001000000000000;
/* 1=Sonne,2=Erde, 3=Mond */
double fsx(double x1, double y1, double x2, double y2, double x3, double y3){
double dv;
dv=G*(me/((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))*((x2-x1)/sqrt((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))+G*(mm/((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)))*((x3-x1)/sqrt((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)));
/*Projektion der Gravitationskraft auf Körper 1 auf x-Richtung */
return dv;
}
double fsy(double x1, double y1, double x2, double y2, double x3, double y3){
double dv=G*(me/((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))*((y2-y1)/sqrt((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))+G*(mm/((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)))*((y3-y1)/sqrt((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1))); /*Projektion der Gravitationskraft auf Körper 1 auf y-Richtung */
return dv;
}
double fex(double x1, double y1, double x2, double y2, double x3, double y3){
double dv;
dv=G*(ms/((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))*((x2-x1)/sqrt((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))+G*(mm/((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)))*((x3-x1)/sqrt((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)));
/*Projektion der Gravitationskraft auf Körper 1 auf x-Richtung */
return dv;
}
double fey(double x1, double y1, double x2, double y2, double x3, double y3){
double dv=G*(ms/((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))*((y2-y1)/sqrt((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))+G*(mm/((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)))*((y3-y1)/sqrt((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1))); /*Projektion der Gravitationskraft auf Körper 1 auf y-Richtung */
return dv;
}
double fmx(double x1, double y1, double x2, double y2, double x3, double y3){
double dv;
dv=G*(me/((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))*((x2-x1)/sqrt((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))+G*(ms/((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)))*((x3-x1)/sqrt((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)));
/*Projektion der Gravitationskraft auf Körper 1 auf x-Richtung */
return dv;
}
double fmy(double x1, double y1, double x2, double y2, double x3, double y3){
double dv=G*(me/((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))*((y2-y1)/sqrt((y2-y1)*(y2-y1)+(x2-x1)*(x2-x1)))+G*(ms/((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1)))*((y3-y1)/sqrt((y3-y1)*(y3-y1)+(x3-x1)*(x3-x1))); /*Projektion der Gravitationskraft auf Körper 1 auf y-Richtung */
return dv;
}
double euler(double (f)(double,double,double,double,double,double), double x1, double y1, double x2, double y2, double x3, double y3, double vi){
double vf=vi + f(x1, y1, x2, y2, x3, y3)*(dt);
return vf;
}
int main() {
float xf3, yf3, tmax;
ofstream out("project1.dat");
ti=0.00000000000000000, xi1=-0.97000436, yi1=0.24308753, vxi1=0.466203685, vyi1=0.43236573, tmax=3.000000000000, xi2=-xi1, vxi2=vxi1, yi2=-yi1, vyi2=vyi1, xi3=0.0000000000000, vxi3=-2*vxi1, yi3=0.0000000000000000, vyi3=-2*vyi1;
for(tf=0.0;tf<=tmax;tf=tf+dt){
vxf1=euler(fsx,xi1,yi1,xi2,yi2,xi3,yi3,vxi1);
vyf1=euler(fsy,xi1,yi1,xi2,yi2,xi3,yi3,vyi1);
vxf2=euler(fex,xi2,yi2,xi1,yi1,xi3,yi3,vxi2); /* Gravitationskraft und damit auch fx und fy sind symmetrisch -> es reichen 2 Funktionen */
vyf2=euler(fey,xi2,yi2,xi1,yi1,xi3,yi3,vyi2);
vxf3=euler(fmx,xi3,yi3,xi2,yi2,xi1,yi1,vxi3);
vyf3=euler(fmy,xi3,yi3,xi2,yi2,xi1,yi1,vyi3);
xf1=xi1+(vxf1)*(dt);
yf1=yi1+(vyf1)*(dt);
xf2=xi2+(vxf2)*(dt);
yf2=yi2+(vyf2)*(dt);
xf3=xi3+(vxf3)*(dt);
yf3=yi3+(vyf3)*(dt);
out << tf << " "<< xf1 << " " << yf1 <<" "<< vxf1 << " " << vyf1<< " "<< xf2 << " " << yf2 <<" "<< vxf2 << " " << vyf2<< " "<< xf3 << " " << yf3 <<" "<< vxf3 << " " << vyf3 << endl;
ti=tf;
vxi1=vxf1;
vyi1=vyf1;
vxi2=vxf2;
vyi2=vyf2;
vxi3=vxf3;
vyi3=vyf3;
xi1=xf1;
yi1=yf1;
xi2=xf2;
yi2=yf2;
xi3=xf3;
yi3=yf3;
}
out.close();
}
grafic（它应该是一个完美的8）：
