4

編集: ODE の NDSolve が Runge Kutta を使用して方程式を解いていることがわかりました。Python コードで Runge Kutta メソッドを使用して、以下の ODE を解くにはどうすればよいですか?

フロート エントリを含むテキスト ファイルに関する私の投稿から、それを判断しpython、mathematicaすぐに 10 からマイナス 6 までの許容範囲で発散を開始することができました。

編集を終了

ここ数時間、Mathematica と Python での私のソリューションが5000何か違う理由を理解しようとしてきましたkm。

数百万秒以上の飛行時間をシミュレートする場合、1 つのプログラムのエラー許容度が高いと私は信じています。

私の質問は、どちらのプログラムがより正確かということです。それが python でない場合、どうすれば精度を調整できますか?

では、私は km 離れている場所Mathematicaから 10 km 未満です。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])

(+明らかに odeint インポートを ode に置き換えます)。

このタイプの ODE では、これはかなり遅いことに注意してください。

dop853 を使用するには、solver.set_integrator('dop853') を使用します。

于 2013-05-05T07:41:10.713 に答える
-1

精度を確認する 1 つの方法は次のとおりです。

2 つの軌跡が互いに十分に異なる時間を取ります (たとえば、軌跡の最後の時間/ポイント)。これを新しい開始点として使用し、最初の時間に戻るように積分します (または、必要に応じて、時間的に前方に積分しますが、シミュレーションに存在するすべてのボディの速度を逆にします)。このシミュレーションの最終ポイントは、元のシミュレーションの最初のポイント (に近い) である必要があります。元の初期点に実際にどれだけ近づいたかを比較することで、Python と Mathematica のソルバー ルーチンがどれだけ優れているかを判断できます。私の推測では、Mathematica ルーチンの方がはるかに優れているので、この Python のことはすべて時間の無駄です。

もう 1 つの方法は次のとおりです。Mathematica を使用すると、補間関数が得られます。この補間関数を記号的に導出し、その導関数を数値的に評価できます。これらの値を ODE のさまざまな点に挿入し、ODE を満たすかどうかを確認します。Python で同様のことを行い、結果を比較します。

于 2013-04-27T15:21:25.593 に答える