C++实现RK4求解1/r²中心力轨迹与scipy.solve_ivp结果不符

我正在编写RK4求解器，用于求解满足Fx = -x/r³、Fy = -y/r³的中心力问题，其中**r=√(x²+y²)**是原点到位置的距离。我将自己的求解结果与scipy.optimise.solve_ivp的输出做对比。

当初始值为(x, y, vx, vy) = (1, 0, 0, 1)时，scipy输出的轨迹显示卫星缓慢落向力源（原点），但我的C++代码输出的是稳定轨迹。我尝试过实现scipy使用的RK45方法测试，但结果没有变化。

我的C++代码如下：

#include <iostream>
#include <vector>
#include <tuple>
#include <cmath>
#include <fstream>

//double const constantForce_x = 1e-3, constantForce_y = 1e-3;
double const constantForce_x = 0., constantForce_y = 0.;

// Function to compute the central force components
std::pair<double, double> central_force(double x, double y) {
    double r = sqrt(x*x + y*y);
    double Fx = -x / (r*r*r) + constantForce_x;
    double Fy = -y / (r*r*r) + constantForce_y;
    return std::make_pair(Fx, Fy);
}

// Function to perform one step of the RK4 method
std::tuple<double, double, double, double> rk4_step(double x, double y, double vx, double vy, double dt) {
    auto k1 = central_force(x, y);
    double k1x = k1.first, k1y = k1.second;
    double k1vx = vx, k1vy = vy;

    auto k2 = central_force(x + k1vx * dt / 2, y + k1vy * dt / 2);
    double k2x = k2.first, k2y = k2.second;
    double k2vx = vx + k1x * dt / 2, k2vy = vy + k1y * dt / 2;

    auto k3 = central_force(x + k2vx * dt / 2, y + k2vy * dt / 2);
    double k3x = k3.first, k3y = k3.second;
    double k3vx = vx + k2x * dt / 2, k3vy = vy + k2y * dt / 2;

    auto k4 = central_force(x + k3vx * dt, y + k3vy * dt);
    double k4x = k4.first, k4y = k4.second;
    double k4vx = vx + k3x * dt, k4vy = vy + k3y * dt;

    x += (k1vx + 2*k2vx + 2*k3vx + k4vx) * dt / 6;
    y += (k1vy + 2*k2vy + 2*k3vy + k4vy) * dt / 6;
    vx += (k1x + 2*k2x + 2*k3x + k4x) * dt / 6;
    vy += (k1y + 2*k2y + 2*k3y + k4y) * dt / 6;

    return std::make_tuple(x, y, vx, vy);
}

// Function to simulate trajectory
std::vector<std::pair<double, double>> simulate_trajectory(double x0, double y0, double vx0, double vy0, double dt, int steps) {
    double x = x0, y = y0, vx = vx0, vy = vy0;
    std::vector<std::pair<double, double>> trajectory;
    trajectory.push_back(std::make_pair(x, y));

    for (int i = 0; i < steps; ++i) {
        std::tie(x, y, vx, vy) = rk4_step(x, y, vx, vy, dt);
        trajectory.push_back(std::make_pair(x, y));
    }

    return trajectory;
}

int main() {
    // Parameters
    double x0 = 1.0;    // initial x position
    double y0 = 0.0;    // initial y position
    double vx0 = 0.0;   // initial x velocity
    double vy0 = 1.0;   // initial y velocity
    double dt = 0.001;   // time step
    int steps = 50000;   // number of steps

    // Simulate trajectory
    auto trajectory = simulate_trajectory(x0, y0, vx0, vy0, dt, steps);

    // Write trajectory to file
    std::ofstream outfile("trajectory.dat");
    for (auto point : trajectory) {
        outfile << point.first << " " << point.second << "\n";
    }
    outfile.close();

    return 0;
}

内容的提问来源于stack exchange，提问作者JamesRWR