4

编辑： 所以我发现 NDSolve for ODE 正在使用 Runge Kutta 来求解方程。如何在我的 python 代码上使用 Runge Kutta 方法来解决下面的 ODE？

从我关于带有浮动条目的文本文件的帖子中，我能够确定这一点，python并mathematica立即开始以 10 到负 6 的容差开始发散。

结束编辑

在过去的几个小时里，我一直试图弄清楚为什么我在 Mathematica 和 Python 中的解决方案会有所5000不同km。

我被引导相信一个程序在模拟超过数百万秒的飞行时间时具有更高的容错性。

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

有了，我离我所在的地方Mathematica不到 10公里。L4Python5947

代码如下：

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"];*)
4

3 回答 3

2

如果此时您的问题已减少到只想使用 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])

（+显然用 ode 替换 odeint 导入）。

请注意，对于这种类型的 ODE，这要慢得多。

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

于 2013-05-05T07:41:10.713 回答
1

我重新编写了 ode 的 def 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
于 2013-05-02T04:58:35.463 回答
-1

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

花一个时间，使 2 个轨迹彼此足够不同（例如，轨迹中的最后一个时间/点）。将此作为新的起点并及时整合到初始时间（或者，如果您愿意，可以在时间上向前整合，但反转模拟中所有物体的速度）。此模拟的最后一点应该（接近）您在原始模拟中的初始点。通过比较您与原始初始点的实际接近程度，您可以判断 python 与 mathematica 求解器例程的优劣。我的猜测是 Mathematica 例程要好得多，所以所有这些 Python 的东西都是浪费时间。

另一种方式可能是这样：使用 Mathematica，您得到的是一个插值函数，您可以用符号方式推导它，然后对其导数进行数值计算。将这些值插入 ODE 的不同点，看看它们是否满足 ODE。在 python 中做类似的事情并比较结果。

于 2013-04-27T15:21:25.593 回答