我在大学的一个项目中，必须使用Python实现Runge-Kutta四阶积分器。我知道我可以使用例如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)
龙格-库塔方法的发明者实际上名叫马丁·威廉·库塔。(Runge 1895做了一些奇怪的事情，Heun 1900让它不那么奇怪，Kutta 1901让它完全灵活和系统。(
您的实现中存在严重的概念错误。
- 您需要将耦合系统视为耦合系统，不能对组件的集成进行解耦。这样最多只能得到一阶积分器
- 在使用
R3时，这一点尤其明显且令人震惊。每个阶段都需要重新计算此值。如果导数向量取决于状态的函数，那么这个值不可能是常数
请参阅在Python中无法让RK4求解轨道体的位置，以及有更好的方法来整理这堆代码吗？它是龙格-库塔四阶，有4个ODE用于工作代码示例。