我在大学的一个项目中，必须使用 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
观察：在我展示代码的第一部分之前，在实现完全正确之前，问题出在积分器函数中，我只是想查看结果，这就是为什么没有计算速度，因为如果我的 r向量是对的，v也是对的。
方程为：r''(向量) = -(GM/r³)*r(vector)
CODE
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)
