Python 中非线性二阶 ODE 的 Rk4 积分器

Question

我在我大学的一个项目中，我必须使用 Python 实现 Runge-Kutta 4 阶积分器。我知道我可以使用 Sympy，但这里的目标是实现方法，代码已经用 Fortran 语言编写，所以基本上我有一个包含正确解决方案值的数据库，我必须在我的代码中获得类似的解决方案.但是，我们有一些问题；我使用线性方程（一阶和二阶）做了几次相同的操作，但这是牛顿万有引力定律的二阶非线性方程。代码没有错误，我的问题是我的代码做错了什么，导致我无法得到正确的结果。

下面我将展示一些预期值和我得到的值，在它们之后我将展示代码。

如果有人能帮助我，我将非常高兴。

正确的结果（预期结果）

  r           t (days)

-12912.5186     .0000
-13135.2914     .0023
-13342.8424     .0046
-13534.9701     .0069
-13711.4971     .0093
-13872.2704     .0116
-14017.1611     .0139
-14146.0643     .0162
-14258.8997     .0185
-14355.6106     .0208
-14436.1641     .0231
-14500.5505     .0255
-14548.7833     .0278
-14580.8984     .0301
-14596.9536     .0324
-14597.0282     .0347
-14581.2222     .0370
-14549.6560     .0394
-14502.4692     .0417
-14439.8201     .0440
-14361.8851     .0463
-14268.8576     .0486
-14160.9475     .0509
-14038.3802     .0532
-13901.3958     .0556
-13750.2482     .0579
-13585.2046     .0602
-13406.5442     .0625
-13214.5576     .0648
-13009.5461     .0671
-12791.8207     .0694
-12561.7015     .0718
-12319.5167     .0741
-12065.6021     .0764
-11800.2999     .0787
-11523.9589     .0810
-11236.9327     .0833
-10939.5799     .0856
-10632.2630     .0880
-10315.3480     .0903
-9989.2038      .0926
-9654.2014      .0949
-9310.7139      .0972
-8959.1154      .0995

错误的结果（来自下面的代码）

  r            t (seconds)

-12912.518615   0.000000
-10894.236220   3600.000000
-8051.384478    7200.000000
-2829.162198    10800.000000
39786.739120    14400.000000
39564.796772    18000.000000
39340.531265    21600.000000
39113.878351    25200.000000
38884.770893    28800.000000
38653.138691    32400.000000
38418.908276    36000.000000
38182.002705    39600.000000
37942.341331    43200.000000
37699.839549    46800.000000
37454.408529    50400.000000
37205.954917    54000.000000
36954.380518    57600.000000
36699.581939    61200.000000
36441.450207    64800.000000
36179.870344    68400.000000
35914.720909    72000.000000
35645.873482    75600.000000
35373.192107    79200.000000
35096.532668    82800.000000
34815.742202    86400.000000

Obs.: 在我展示代码之前，它的第一部分在实现完全正确之前，问题出在积分器函数中，我只是想看看结果，这就是为什么速度不是计算是因为如果我的 r 向量是正确的，v 也将是正确的。等式是： r''（向量）= -（GM/r³）*r（向量）

代码

import numpy as np

# alternative to not typing all the time

TINTE = 5           #days 
a = 26551.0         #kilometers
e = 0.1             
i = 55              #degrees
OM = 102            #degrees
w = 32              #degrees
f = 12              #degrees

# Mass of central body

Mc = 5.97240e+24     #kg (Earth = 7.97240D+24    Sol = 1.98850D+30)
M2 = 5.97240e+24     #kg (Earth = 7.97240D+24    Sol = 1.98850D+30)
M3 = 7.34600e+22     #kg Mass of the Moon
G = 6.67408e-20      #Value prepared for km
#Mi = Mc/(M2+M3)      #G*Mc - alternatively
#PI = math.acos(-1.0) 
TN = 27.321660       #Time converter

# Dados do Sistema

tempo = list()
xc = list()
yc = list()
zc = list()

#Transformation of orbital elements in position and velocity in the ECI coordinate system

P = a*(1-e**2)
R = P/(1+e*(np.cos(np.deg2rad(f))))

X = list()
X.append(R*((np.cos(np.deg2rad(OM)))*(np.cos(np.deg2rad(w+f))) - (np.sin(np.deg2rad(OM)))*(np.cos(np.deg2rad(i)))*(np.sin(np.deg2rad(w+f)))))
X.append(R*((np.sin(np.deg2rad(OM)))*(np.cos(np.deg2rad(w+f))) + (np.cos(np.deg2rad(OM)))*(np.cos(np.deg2rad(i)))*(np.sin(np.deg2rad(w+f)))))
X.append(R*(np.sin(np.deg2rad(i)))*(np.sin(np.deg2rad(w+f))))

V = list()
V.append((-(Mi/P)**0.5)*((np.cos(np.deg2rad(OM)))*((np.sin(np.deg2rad(w+f)))+e*(np.sin(np.deg2rad(w)))) + (np.sin(np.deg2rad(OM)))*(np.cos(np.deg2rad(i)))*((np.cos(np.deg2rad(w+f))) + e*(np.cos(np.deg2rad(w))))))
V.append((-(Mi/P)**0.5)*((np.sin(np.deg2rad(OM)))*((np.sin(np.deg2rad(w+f)))+e*(np.sin(np.deg2rad(w)))) - (np.cos(np.deg2rad(OM)))*(np.cos(np.deg2rad(i)))*((np.cos(np.deg2rad(w+f))) + e*(np.cos(np.deg2rad(w))))))
V.append(((Mi/P)**0.5)*((np.sin(np.deg2rad(i)))*(np.cos(np.deg2rad(w+f)))+e*(np.cos(np.deg2rad(w)))))

Vp = (V[0]**2 + V[1]**2 + V[2]**2)**0.5

xc.append(X[0])
yc.append(X[1])
zc.append(X[2])

Vx = V[0]
Vy = V[1]
Vz = V[2]

def RUNGE_KUTAH_4(X,V):
    
    #variables
    RT = 6370                  #km
    G = 6.67408e-20            #Value prepared for km
    p = X
    ç = V
    R = ( p[0]**2 + p[1]**2 + p[2]**2 )**0.5
    R3 = R*R*R
    Ve = Vp
        
    # initial state
    tempo.append(0)
    t = 0
    r1 = p[0]
    r2 = p[1]
    r3 = p[2]
    u1 = ç[0]
    u2 = ç[1]
    u3 = ç[2]
    
    #step
    delta_t = 3600
    
    def rk4(r,u,R3):
        m1 = u
        k1 = -((G*Mc)/(R3))*r
        m2 = u + 0.5*delta_t*k1
        t_2 = t + 0.5*delta_t
        r_2 = r + 0.5*delta_t*m1
        u_2 = m2
        k2 = -((G*Mc)/(R3))*r
        m3 = u + 0.5*delta_t*k2
        t_3 = t + 0.5*delta_t
        r_3 = r + 0.5*delta_t*m2
        u_3 = m3
        k3 = -((G*Mc)/(R3))*r
        m4 = u + 0.5*delta_t*k3
        t_4 = t + delta_t
        r_4 = r + delta_t*m3
        u_4 = m4
        k4 = -((G*Mc)/(R3))*r
        
        r = r + (delta_t/6)*(m1+2*(m2+m3)+m4)
        u = u + (delta_t/6)*(k1+2*(k2+k3)+k4)
        
        return [r,u]
        

    # step by step solution 
    lim = 86400*TINTE
    while t < lim:
        r1 = rk4(r1,u1,R3)[0]
        r2 = rk4(r2,u2,R3)[0]
        r3 = rk4(r3,u3,R3)[0]
        
        R = (r1**2 + r2**2 + r3**2)**0.5
        R3 = R*R*R
        t += delta_t
        
        tempo.append(t)
        xc.append(r1)

#-------------------------------------------------------------------------------------------------------------------------------

RUNGE_KUTAH_4(X,V)

Answer 1

龙格-库塔方法的发明者真名叫马丁·威廉·库塔。（Runge 1895 做了一些奇怪的事情，Heun 1900 让它变得不那么奇怪，Kutta 1901 让它变得完全灵活和系统。）

您的实施存在严重的概念错误。

你需要把一个耦合系统当作一个耦合系统来对待，你不能解耦组件的集成。充其量您将通过这种方式获得一阶积分器。
这在您使用 R3 时尤为明显和令人震惊。每个阶段都需要重新计算该值。如果导数向量依赖于一个状态函数，那么这个值就不能是常量。

有关工作代码示例，请参阅 Cannot get RK4 to solve for position of orbiting body in Python and Is there a better way to tidy this chuck of code? It is the Runge-Kutta 4th order with 4 ODEs。

Python 中非线性二阶 ODE 的 Rk4 积分器

Rk4 Integrator of a nonlinear second order ODE in Python

python

ode

runge-kutta

orbital-mechanics

代码