I'm in a project in my University and I have to implement the Runge-Kutta 4-order integrator using Python. I know I can use for example Sympy, but the goal here is implement the method, the code is ready written in Fortran language, so basically I have a data base with the correct values of solutions and I have to get similar solutions in my code. However, we have some problems; I did the same several times using linear equations ( first and second order ), however this is a second order nonlinear equation from Newton universal law of gravitation. The code has no error, my problem is what my code is doing wrong that I cannot get the right results.
Below I'll show some of the expected values and the ones i'm getting, after them I'll show the code.
I'd be really pleasured if someone could help me.
RIGHT RESULTS (expected results)
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
WRONG RESULTS (from the code below)
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.: Before I show the code the first part of it before the implementation is completely right, the problem is in the integrator function, I'm just trying to see the results, that's why the velocity is not being computed since if my r vector is right, v will be as well. The equation is : r''(vector) = -(GM/r³)*r(vector)
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)
The inventor of the Runge-Kutta methods was really named Martin Wilhelm Kutta. (Runge 1895 did something strange, Heun 1900 made it less strange, Kutta 1901 made it fully flexible and systematic.)
You have a severe conceptual error in your implementation.
R3. This value needs to be re-computed for every stage. If the derivatives vector depends on a function of the state, then this value can not be constant.See 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 for working code examples.