在 C++ 中对与 1/r^2 成比例的力进行 Runge-Kutta 4 编码，输出的轨迹与 scipy.optimise.solve_ivp 不同

Question

我正在尝试为力 Fx = -x/r^3 和 Fy = -y/r^3 编写 RK4 解算器，其中 r=sqrt(x^2+y^2) 是原点和位置之间的距离。我将我的结果与 scipy.optimise.solve_ivp 的结果进行比较。

对于初始值 (x, y, vx, vy) = (1, 0, 0, 1)，scipy 输出卫星缓慢落向力源（原点）的轨迹。然而，我的代码输出了稳定的轨迹。有什么解决办法吗？我附上了一张图表，显示了两个输出的轨迹，突出显示了它们的不同之处。

#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;
}

我尝试使用我的 scipy 编写 RK45 方法来测试是否存在差异，但没有任何变化。

Answer 1

solve_ivp 将很好地解决这个问题，但是您必须包含可选参数 rtol （或 atol）来设置预期的相对（或绝对）容差，以便它使用的时间步长的自动调整将是足够的。

在下面的代码中我设置了 rtol=1.0e-6。但是，默认值（如果不包含此参数）是 1.0e-3，如果您尝试一下，您会发现该解决方案是垃圾（并且缓慢螺旋式上升）。

import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp


def deriv( t, Y ):
    x, y, vx, vy = Y
    rcubed = ( x ** 2 + y ** 2 ) ** 1.5      # note: r/r^3 --> 1/r^2, so inverse-square law
    return np.array( [ vx, vy, -x / rcubed, -y / rcubed ] )


r = 1.0
v = 1.0 / math.sqrt( r )               # for a circular orbit
period = 2 * math.pi / ( v / r )       # ditto
tmax = 10 * period
Y0 = np.array( ( r, 0, 0, v ) )        # initial x, y, vx, vy
sol = solve_ivp( deriv, [ 0.0, tmax ], Y0, rtol=1.0e-6 )

x = np.array( sol.y[0] )
y = np.array( sol.y[1] )
plt.plot( x, y )
plt.axis( 'scaled' )
plt.xlabel( 'x' );   plt.ylabel( 'y' )
plt.show()