我正在尝试为力 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 方法来测试是否存在差异,但没有任何变化。
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()