Python中非线性二阶常微分方程的Rk4积分器 [英] Rk4 Integrator of a nonlinear second order ODE in Python
问题描述
我在我大学的一个项目中，我必须使用 Python 实现 Runge-Kutta 4 阶积分器.我知道我可以使用例如 Sympy，但这里的目标是实现该方法，代码已准备好用 Fortran 语言编写，所以基本上我有一个包含正确解决方案值的数据库，我必须在我的代码中获得类似的解决方案.但是，我们有一些问题；我使用线性方程(一阶和二阶)做了几次相同的事情，但这是牛顿万有引力定律中的二阶非线性方程.代码没有错误，我的问题是我的代码做错了什么，无法得到正确的结果.
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.:在我展示代码之前，它的第一部分在实现之前是完全正确的，问题在于积分器函数，我只是想看看结果，这就是为什么没有计算速度，因为如果我的 r 向量是正确的，v 也会正确.等式是:r''(vector) = -(GM/r³)*r(vector)
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)
推荐答案
Runge-Kutta 方法的发明者实际上是 Martin Wilhelm Kutta.(Runge 1895 做了一些奇怪的事情，Heun 1900 让它变得不那么奇怪，Kutta 1901 让它变得完全灵活和系统.)
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时尤其明显且令人震惊.这个值需要为每个阶段重新计算.如果导数向量取决于状态的函数，则该值不能为常数.
- You need to treat a coupled system as a coupled system, you can not de-couple the integration of the components. At best you will get a first order integrator this way.
- This is especially visible and egregious in your use of
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.
参见 无法让 RK4 求解Python 中的轨道体 和 是有更好的方法来整理这堆代码吗?这是 Runge-Kutta 4 阶，带有 4 个 ODE，用于工作代码示例.
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.
这篇关于Python中非线性二阶常微分方程的Rk4积分器的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
