Python与Mathematica中的ODE集成结果比较

Question

编辑： 我发现ODE的NDSolve使用了龙格库塔算法来解决方程。 我如何在我的Python代码中使用龙格库塔方法来解决下面的ODE？

从我在文本文件与浮点条目上发布的内容中，我能够确定python和mathematica立即以10的负6次方的公差开始发散。

结束编辑

在过去几个小时里，我一直在试图找出为什么我的Mathematica和Python中的解决方案会相差5000多公里。

我认为一个程序在模拟数百万秒的飞行时间时具有更高的误差容限。

我的问题是哪个程序更精确，如果不是Python，我该如何调整精度？

使用Mathematica，我距离L4不到10km，而使用Python，我距离L4有5947km。

以下是代码：

Python

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy import linspace
from scipy.optimize import brentq

me = 5.974 * 10 ** (24)  #  mass of the earth                                     
mm = 7.348 * 10 ** (22)  #  mass of the moon                                      
G = 6.67259 * 10 ** (-20)  #  gravitational parameter                             
re = 6378.0  #  radius of the earth in km                                         
rm = 1737.0  #  radius of the moon in km                                          
r12 = 384400.0  #  distance between the CoM of the earth and moon                 
d = 300 #  distance the spacecraft is above the Earth                             
pi1 = me / (me + mm)
pi2 = mm / (me + mm)
mue = 398600.0  #  gravitational parameter of earth km^3/sec^2                    
mum = G * mm  #  grav param of the moon                                           
mu = mue + mum
omega = np.sqrt(mu / (r12 ** 3))

nu = -np.pi / 4  #  true anomaly  pick yourself                                   

xl4 = r12 / 2 - 4671  #  x location of L4                                         
yl4 = np.sqrt(3) / 2 * r12  #  y                                                  

print("The location of L4 is", xl4, yl4)

#  Solve for Jacobi's constant                                                    
def f(C):
    return (omega ** 2 * (xl4 ** 2 + yl4 ** 2) + 2 * mue / r12 + 2 * mum / r12
            + 2 * C)


c = brentq(f, -5, 0)

print("Jacobi's constant is",c)

x0 = (re + 200) * np.cos(nu) - pi2 * r12  #  x location of the satellite          
y0 = (re + 200) * np.sin(nu)  #  y location                                       

print("The satellite's initial position is", x0, y0)

vbo = (np.sqrt(omega ** 2 * (x0 ** 2 + y0 ** 2) + 2 * mue /
               np.sqrt((x0 + pi2 * r12) ** 2 + y0 ** 2) + 2 * mum /
               np.sqrt((x0 - pi1 * r12) ** 2 + y0 ** 2) + 2 * -1.21))

print("Burnout velocity is", vbo)

gamma = 0.4678 * np.pi / 180  #  flight path angle pick yourself                  

vx = vbo * (np.sin(gamma) * np.cos(nu) - np.cos(gamma) * np.sin(nu))
#  velocity of the bo in the x direction                                          
vy = vbo * (np.sin(gamma) * np.sin(nu) + np.cos(gamma) * np.cos(nu))
#  velocity of the bo in the y direction                                          

print("The satellite's initial velocity is", vx, vy)

#  r0 = [x, y, 0]                                                                 
#  v0 = [vx, vy, 0]                                                               
u0 = [x0, y0, 0, vx, vy, 0]


def deriv(u, dt):
return [u[3],  #  dotu[0] = u[3]                                                 
        u[4],  #  dotu[1] = u[4]                                                 
        u[5],  #  dotu[2] = u[5]                                                 
        (2 * omega * u[4] + omega ** 2 * u[0] - mue * (u[0] + pi2 * r12) /
         np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum *
         (u[0] - pi1 * r12) /
         np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
        #  dotu[3] = that                                                        
        (-2 * omega * u[3] + omega ** 2 * u[1] - mue * u[1] /
         np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum * u[1] /
         np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
        #  dotu[4] = that                                                        
        0]  #  dotu[5] = 0                                                       


dt = np.linspace(0.0, 6.0 * 86400.0, 2000000.0)  #  secs to run the simulation    
u = odeint(deriv, u0, dt)
x, y, z, x2, y2, z2 = u.T

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z, color = 'r')
#  adding the Lagrange point                                                      
phi = np.linspace(0, 2 * np.pi, 100)
theta = np.linspace(0, np.pi, 100)
xm = 2000 * np.outer(np.cos(phi), np.sin(theta)) + xl4
ym = 2000 * np.outer(np.sin(phi), np.sin(theta)) + yl4
zm = 2000 * np.outer(np.ones(np.size(phi)), np.cos(theta))
ax.plot_surface(xm, ym, zm, color = '#696969', linewidth = 0)
ax.auto_scale_xyz([-8000, 385000], [-8000, 385000], [-8000, 385000])
#  adding the earth                                                               
phi = np.linspace(0, 2 * np.pi, 100)
theta = np.linspace(0, np.pi, 100)
xm = 2000 * np.outer(np.cos(phi), np.sin(theta))
ym = 2000 * np.outer(np.sin(phi), np.sin(theta))
zm = 2000 * np.outer(np.ones(np.size(phi)), np.cos(theta))
ax.plot_surface(xm, ym, zm, color = '#696969', linewidth = 0)
ax.auto_scale_xyz([-8000, 385000], [-8000, 385000], [-8000, 385000])

plt.show()

#  The code below finds the distance between path and l4                          
my_x, my_y, my_z = (xl4, yl4, 0.0)

delta_x = x - my_x
delta_y = y - my_y
delta_z = z - my_z
distance = np.array([np.sqrt(delta_x ** 2 + delta_y ** 2 + delta_z ** 2)])

minimum = np.amin(distance)

print("Closet approach to L4 is", minimum)

Mathematica

ClearAll["Global`*"];
me = 5.974*10^(24);
mm = 7.348*10^(22);
G = 6.67259*10^(-20);
re = 6378;
rm = 1737;
r12 = 384400;

\[Pi]1 = me/(me + mm);
\[Pi]2 = mm/(me + mm);
M = me + mm;
\[Mu]1 = 398600;
\[Mu]2 = G*mm;
\[Mu] = \[Mu]1 + \[Mu]2;
\[CapitalOmega] = Sqrt[\[Mu]/r12^3];
\[Nu] = -\[Pi]/4;

xl4 = 384400/2 - 4671;
yl4 = Sqrt[3]/2*384400 // N;

Solve[\[CapitalOmega]^2*(xl4^2 + yl4^2) + 2 \[Mu]1/r12 + 
   2 \[Mu]2/r12 + 2*C == 0, C]
x = (re + 200)*Cos[\[Nu]] - \[Pi]2*r12 // N
y = (re + 200)*Sin[\[Nu]] // N

{{C -> -1.56824}}

-19.3098

-4651.35

vbo = Sqrt[\[CapitalOmega]^2*((x)^2 + (y)^2) + 
   2*\[Mu]1/Sqrt[(x + \[Pi]2*r12)^2 + (y)^2] + 
   2*\[Mu]2/Sqrt[(x - \[Pi]1*r12)^2 + (y)^2] + 2*(-1.21)]

10.8994

\[Gamma] = 0.4678*Pi/180;
vx = vbo*(Sin[\[Gamma]]*Cos[\[Nu]] - Cos[\[Gamma]]*Sin[\[Nu]]);
vy = vbo*(Sin[\[Gamma]]*Sin[\[Nu]] + Cos[\[Gamma]]*Cos[\[Nu]]);

r0 = {x, y, 0};
v0 = {vx, vy, 0}

{7.76974, 7.64389, 0}

s = NDSolve[{x1''[t] - 
      2*\[CapitalOmega]*x2'[t] - \[CapitalOmega]^2*
       x1[t] == -\[Mu]1/((Sqrt[(x1[t] + \[Pi]2*r12)^2 + 
             x2[t]^2])^3)*(x1[t] + \[Pi]2*
          r12) - \[Mu]2/((Sqrt[(x1[t] - \[Pi]1*r12)^2 + 
             x2[t]^2])^3)*(x1[t] - \[Pi]1*r12), 
    x2''[t] + 
      2*\[CapitalOmega]*x1'[t] - \[CapitalOmega]^2*
       x2[t] == -\[Mu]1/(Sqrt[(x1[t] + \[Pi]2*r12)^2 + x2[t]^2])^3*
       x2[t] - \[Mu]2/(Sqrt[(x1[t] - \[Pi]1*r12)^2 + x2[t]^2])^3*
       x2[t], x3''[t] == 0, x1[0] == r0[[1]], x1'[0] == v0[[1]], 
    x2[0] == r0[[2]], x2'[0] == v0[[2]], x3[0] == r0[[3]], 
    x3'[0] == v0[[3]]}, {x1, x2, x3}, {t, 0, 1000000}];

ParametricPlot3D[
 Evaluate[{x1[t], x2[t], x3[t]} /. s], {t, 0, 10*24*3600}, 
 PlotStyle -> {Red, Thick}]

g1 = ParametricPlot3D[
   Evaluate[{x1[t], x2[t], x3[t]} /. s], {t, 0, 5.75*3600*24}, 
   PlotStyle -> {Red}, 
   PlotRange -> {{-10000, 400000}, {-10000, 400000}}];
g2 = Graphics3D[{Blue, Opacity[0.6], Sphere[{-4671, 0, 0}, re]}];
g3 = Graphics3D[{Green, Opacity[0.6], Sphere[{379729, 0, 0}, rm]}];
g4 = Graphics3D[{Black, Sphere[{xl4, yl4, 0}, 2000]}];
Show[g2, g1, g3, g4, Boxed -> False]


(*XYdata=Flatten[Table[Evaluate[{x1[t],x2[t],x3[t]}/.s],{t,5.5*24*\
3600,5.78*24*3600,1}],1];
X1Y1data=Flatten[Table[Evaluate[{x1'[t],x2'[t],x3'[t]}/.s],{t,5.5*24*\
3600,5.78*24*3600,1}],1];
SetDirectory[NotebookDirectory[]];
Export["OrbitData.txt",XYdata,"CSV"];
Export["OrbVeloc.txt",X1Y1data,"CSV"];*)

- dustin

---

我的建议是找一个具有精确分析解的问题，并用两种方法解决它。 - agentp

你可以尝试调整 odeint 的 atol 和 rtol 参数（参见 http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html#scipy.integrate.odeint）。但这并不能告诉你 odeint 解决方案的精度是比 Mathematica 解决方案更高还是更低。 - Warren Weckesser

你可能想尝试的一件事是明确告诉NDSolve使用什么方法。odeint和Mathematica都是动态/自适应选择解决方案方法，因此你很难得到完全相同的结果。如果没有仔细研究问题，很难说你看到的差异是否显著或者只是噪声。 - agentp

@WarrenWeckesser 我将 atol 和 rtol 改为 1e-14，但轨迹仍然保持在 5xxx 公里之外。 - dustin

@george 我发现 NDSolve 使用了龙格-库塔算法。我该如何在 Python 代码中使用龙格-库塔算法来解决常微分方程？ - dustin

Answer 1

如果您现在的问题只是想使用Runge-Kutta方法，您可以例如替换这一行代码：

u = odeint(deriv, u0, dt)

使用类似以下内容的方式：

#reverse the order of arguments
def deriv2(t,u):
    return deriv(u,t)

# initialize a 4th order Runge-Kutta ODE solver
solver = ode(deriv2).set_integrator('dopri5')
solver.set_initial_value(u0)
u = np.empty((len(dt), 6))
u[0,:] = u0
for ii in range(1,len(dt)):
    u[ii] = solver.integrate(dt[ii])

(+显然要将odeint导入替换为ode)。

请注意，对于这种类型的ODE，这种方法的速度明显较慢。

要使用dop853，请使用solver.set_integrator('dop853')。

- HenriV

---

也许有一种方法可以让Cython运行重型计算以加快速度？ - dustin

是的，那可能有所帮助，但我只是开始研究性能调优。请记住，RK4是一个表现良好的求解器（至少对于非刚性问题），您可能不需要那么多时间步长。 - HenriV

Answer 2

我重新编写了 ODE 的 deriv 部分，现在它可以工作了！因此，Mathematica 绘图和 Python 得出的结果一致。

def deriv(u, dt):
    return [u[3],  #  dotu[0] = u[3]                                                 
            u[4],  #  dotu[1] = u[4]                                                 
            u[5],  #  dotu[2] = u[5]                                                 
            (2 * omega * u[4] + omega ** 2 * u[0] - mue * (u[0] + pi2 * r12) /
             np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum *
             (u[0] - pi1 * r12) /
             np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
            #  dotu[3] = that                                                        
            (-2 * omega * u[3] + omega ** 2 * u[1] - mue * u[1] /
             np.sqrt(((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) - mum * u[1] /
             np.sqrt(((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)),
            #  dotu[4] = that                                                        
            0]  #  dotu[5] = 0      

- dustin

Answer 3

检查准确性的一种方法可能是这样的：

选择两条轨迹在某个时间点或位置上有足够的差异（例如轨迹的最后一个时间/点）。将此作为新的起点，并向初始时间反向积分（或者，如果您愿意，向前积分，但要反转模拟中存在的所有物体的速度）。此模拟的最终点应该接近于原始模拟中的初始点。通过比较实际到达原始初始点的距离，您可以判断 Python 与 Mathematica 求解器程序的好坏。我猜想 Mathematica 程序要好得多，所以所有这些 Python 的东西都是浪费时间。

另一种方法可能是：

使用 Mathematica 可以获得一个插值函数，您可以对其进行符号导数并在不同点数值地评估其导数。将这些值插入 ODE 中的不同点并查看它们是否满足 ODE。在 Python 中执行类似操作并比较结果。

- magma

---

第三个选项：将这个问题移动到（几乎）新的 http://mathematica.stackexchange.com/ 姐妹站点，那里每天都有很酷的 MMA 想法交流。 - magma

我找到了问题所在。NDSolve使用龙格-库塔算法进行积分。昨天我编辑了我的帖子，询问如何在Python代码中实现此功能。 - dustin

你是如何确定的？在我看来，默认值很可能是Adams（与odeint相同）。坦白地说，这是Mathematica的一个恼人之处，因为通常没有简单的方法可以知道它选择了什么自动设置。但是，你可以手动指定方法并查看是否产生了显着不同的结果。 - agentp