Rk4 Интегратор нелинейного ОДУ второго порядка на Python
Я участвую в проекте в своем университете, и мне нужно реализовать интегратор 4-го порядка Runge-Kutta с использованием Python. Я знаю, что могу использовать, например, Sympy, но цель здесь - реализовать метод, код готов, написан на языке Фортран, поэтому в основном у меня есть база данных с правильными значениями решений, и я должен получить аналогичные решения в моем коде. Однако у нас есть некоторые проблемы; Я делал то же самое несколько раз, используя линейные уравнения (первого и второго порядка), однако это нелинейное уравнение второго порядка из универсального закона всемирного тяготения Ньютона. В коде нет ошибок, моя проблема в том, что мой код делает неправильно, что я не могу получить правильные результаты.
Ниже я покажу некоторые ожидаемые значения и те, которые я получаю, после них я покажу код.
Я был бы очень рад, если бы мне кто-нибудь помог.
ПРАВИЛЬНЫЕ РЕЗУЛЬТАТЫ (ожидаемые результаты)
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(вектор)
КОД
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)
1 ответ
Изобретателем методов Рунге-Кутта действительно был Мартин Вильгельм Кутта. (Runge 1895 сделал что-то странное, Heun 1900 сделал его менее странным, Kutta 1901 сделал его полностью гибким и систематическим.)
У вас серьезная концептуальная ошибка в вашей реализации.
- Вам нужно рассматривать связанную систему как связанную систему, вы не можете разъединить интеграцию компонентов. Так в лучшем случае вы получите интегратора первого порядка.
- Это особенно заметно и вопиюще при использовании вами
R3. Это значение необходимо пересчитывать для каждого этапа. Если вектор производных зависит от функции состояния, то это значение не может быть постоянным.
См. Не удается заставить RK4 определить положение вращающегося тела в Python и Есть ли лучший способ привести в порядок этот кусок кода? Это 4-й порядок Рунге-Кутты с 4 ODE для рабочих примеров кода.