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 (将#修改为@）