我正在尝试评估在三个维度中绕着一个中心运行的 2 个天体之间的扰动幅度。为了做到这一点，我需要对错误进行估计，我使用 Richardson 外推法进行了估计，如此处所述。我正在使用 Runge-Kutta 4 和velocity Verlet。的时间步长分别为 Verlet 和 Runge-Kutta获得了错误和$O(10^{-7})$$O(10^{-12})$$0.001$. 当我考虑扰动时，只有速度 Verlet 仍然给我预期的误差，而 Runge-Kutta 给我一个 Verlet 相同阶数的误差，使用相同的先前时间步长。即使我减少或增加时间步长，也会发生同样的情况。实际上，对于受干扰的问题，我已经预料到与错误有关的一些不好的事情，因为这是我为 Runge-Kutta 的理查森分数获得的

可以看出，图中的水平应该是左右，但实际上是。这种尖峰也出现在速度 Verlet 的同一图中，但对于后者，它们的幅度远小于上图中的幅度，并且值大约为，正如 Verlet 所预期的那样。我检查了所有的 Runge-Kutta 代码和 ode 系统代码，但没有发现任何错误。此外，我使用更简单的函数测试了这两种方法，获得了正确的结果，因此我无法理解 Runge-Kutta 出现这种错误估计的原因。$2^p = 16$$2$$F_h$$4$

编辑 1

功能代码如下

vector<double> F(double t, vector<double> x, vector<CB> objs, vector<double> s) {
    //objs store the masses of the objects
    //s store the positions of the objects in the order [x1,y1,z1,x2,y2,z2,...];
    //x store the positions and velocity of the target (perturbed) object in the order [x,Vx,y,Vy,z,Vz]
    vector<double> p(x.size(),0);
    double m, d;
    vector<double> y(x.size()/2,0); 

    for (int i = 0; i < objs.size(); i++) {
        //temporary store the mass of the i-th object
        m = objs[i].CB::getmass(); 
        //temporary store the position of the i-th object from the s vector
        y = {s[i*x.size()/2], s[i*x.size()/2 + 1], s[i*x.size()/2 + 2]}; 

        d = sqrt(pow((x[0] - y[0]),2) + pow((x[2] - y[1]),2) + pow((x[4] - y[2]),2));

        //x system
        p[0] = x[1];
        p[1] += -G*m*(x[0] - y[0])/pow(d,3); 
        //y system                           
        p[2] = x[3];
        p[3] += -G*m*(x[2] - y[1])/pow(d,3);
        //z system
        p[4] = x[5];
        p[5] += -G*m*(x[4] - y[2])/pow(d,3);
    }

    return p;
}