如何解决scipy.integrate.solve_ivp 有问题
我正在尝试使用 scipy.integrate.solve_ivp 为我的 n 体模拟计算牛顿引力方程的解,但是我很困惑如何将函数传递到 solve_ivp。我有以下代码:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
G = 6.67408e-11
m_sun = 1988500e24
m_jupiter = 1898.13e24
m_earth = 5.97219e24
au = 149597870.700e3
v_factor = 1731460
year = 31557600.e0
init_s = np.array([-6.534087946884256E-03*au,6.100454846284101E-03*au,1.019968145073305E-04*au,-6.938967653087248E-06*v_factor,-5.599052606952444E-06*v_factor,2.173251724105919E-07*v_factor])
init_j = np.array([2.932487231769548E+00*au,-4.163444383137574E+00*au,-4.833604407653648E-02*au,6.076788230491844E-03*v_factor,4.702729516645153E-03*v_factor,-1.554436340872727E-04*v_factor])
variables_s = init_s
variables_j = init_j
N = 2
tStart = 0e0
t_End = 25*year
Nt = 2000
dt = t_End/Nt
temp_end = dt
t=tStart
domain = (t,temp_end)
planetsinit = np.vstack((init_s,init_j))
planetspos = planetsinit[:,0:3]
mass = np.vstack((1988500e24,1898.13e24))
def weird_division(n,d):
return n / d if d else 0
variables_save = np.zeros((N,6,Nt))
variables_save[:,:,0] = planetsinit
pos_s = planetspos[0]
pos_j = planetspos[1]
while t < t_End:
t_index = int(weird_division(t,dt))
for index in range(len(planetspos)):
for otherindex in range(len(planetspos)):
if index != otherindex:
x1_p1,x2_p1,x3_p1 = planetsinit[index,0:3]
x1_p2,x2_p2,x3_p2 = planetsinit[otherindex,0:3]
m = mass[otherindex]
def f_grav(t,y):
x1_p1,x3_p1,v1_p1,v2_p1,v3_p1 = y
x1_diff = x1_p1 - x1_p2
x2_diff = x2_p1 - x2_p2
x3_diff = x3_p1 - x3_p2
dydt = [v1_p1,v3_p1,-(x1_diff)*G*m/((x1_diff)**2+(x2_diff)**2+(x3_diff)**2)**(3/2),-(x2_diff)*G*m/((x1_diff)**2+(x2_diff)**2+(x3_diff)**2)**(3/2),-(x3_diff)*G*m/((x1_diff)**2+(x2_diff)**2+(x3_diff)**2)**(3/2)]
return dydt
solution = solve_ivp(fun=f_grav,t_span=domain,y0=planetsinit[index])
planetsinit[index] = solution['y'][0:6,-1]
variables_save[index,t_index] = solution['y'][0:6,-1]
planetspos[index] = planetsinit[index][0:3]
t += dt
temp_end += dt
domain = (t,temp_end)
pos_s = variables_save[0,0:3,:]
pos_j = variables_save[1,:]
plt.plot(variables_save[0,:][0],variables_save[0,:][1])
plt.plot(variables_save[1,variables_save[1,:][1])
上面的代码工作得很好,并产生了一个稳定的轨道。但是,当我计算函数外的加速度并将其馈入 f_grav 函数时,出现问题并产生不再稳定的轨道。但是我很困惑,因为我不知道为什么数据不同,因为我似乎已经通过了完全相同的输入。这让我认为这可能是函数 f_grav 被solve_ivp 积分器插值的方式?为了计算加速度,我所做的就是将循环中的以下代码更改为:
x1_p1,0:3]
x1_p2,0:3]
m = mass[otherindex]
x1_diff = x1_p1 - x1_p2
x2_diff = x2_p1 - x2_p2
x3_diff = x3_p1 - x3_p2
ax = -(x1_diff)*G*m/((x1_diff)**2+(x2_diff)**2+(x3_diff)**2)**(3/2)
ay = -(x2_diff)*G*m/((x1_diff)**2+(x2_diff)**2+(x3_diff)**2)**(3/2)
az = -(x3_diff)*G*m/((x1_diff)**2+(x2_diff)**2+(x3_diff)**2)**(3/2)
def f_grav(t,y):
x1_p1,v3_p1 = y
dydt = [v1_p1,ax,ay,az]
return dydt
solution = solve_ivp(fun=f_grav,y0=planetsinit[index])
planetsinit[index] = solution['y'][0:6,-1]
variables_save[index,-1]
planetspos[index] = planetsinit[index][0:3]
正如我所说,我不知道为什么会产生不同的轨道,如下所示,任何有关为什么或如何解决它的提示,我都将不胜感激。为了澄清为什么我不能按原样使用工作代码,因为当涉及更多的物体时,我的目标是对所有其他行星的加速度贡献求和,这是在函数本身中计算加速度的方式不可能的。
对不起,代码块很大,但我确实觉得它是合适的,因为它可以运行并且问题本身更清晰。
两者的时间周期相同,dt,但是左边的轨道是稳定的,右边的不是
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。