solve_ivp 错误：“所需的步长小于数字之间的间距。”

Question

一直在尝试使用来自 scipy 的 RK45 解决牛顿二体问题，但一直遇到 TypeError：'所需的步长小于数字之间的间距。'我尝试了与下面不同的 t_eval 值，但似乎没有任何效果。

from scipy import optimize
from numpy import linalg as LA
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import numpy as np
from scipy.integrate import solve_ivp

AU=1.5e11
a=AU
e=0.5
mss=2E30
ms = 2E30
me = 5.98E24
mv=4.867E24
yr=3.15e7
h=100
mu1=ms*me/(ms+me)
mu2=ms*me/(ms+me)
G=6.67E11
step=24

vi=np.sqrt(G*ms*(2/(a*(1-e))-1/a))
#sun=sphere(pos=vec(0,0,0),radius=0.1*AU,color=color.yellow)
#earth=sphere(pos=vec(1*AU,0,0),radius=0.1*AU)

sunpos=np.array([-903482.12391302, -6896293.6960525, 0.  ])
earthpos=np.array([a*(1-e),0,0])

earthv=np.array([0,vi,0])
sunv=np.array([0,0,0])





def accelerations2(t,pos):
    norme=sum( (pos[0:3]-pos[3:6])**2 )**0.5
    gravit = G*(pos[0:3]-pos[3:6])/norme**3
    sunaa = me*gravit
    earthaa = -ms*gravit
    tota=earthaa+sunaa
    return [*earthaa,*sunaa]

def ode45(f,t,y,h):
        """Calculate next step of an initial value problem (IVP) of an ODE with a RHS described
        by the RHS function with an order 4 approx. and an order 5 approx.
        Parameters:
        t: float. Current time.
        y: float. Current step (position).
        h: float. Step-length.
        Returns:
        q: float. Order 2 approx.
        w: float. Order 3 approx.
        """

        s1 = f(t, y[0],y[1])
        s2 = f(t + h/4.0, y[0] + h*s1[0]/4.0,y[1] + h*s1[1]/4.0)
        s3 = f(t + 3.0*h/8.0, y[0] + 3.0*h*s1[0]/32.0 + 9.0*h*s2[0]/32.0,y[1] + 3.0*h*s1[1]/32.0 + 9.0*h*s2[1]/32.0)
        s4 = f(t + 12.0*h/13.0, y[0] + 1932.0*h*s1[0]/2197.0 - 7200.0*h*s2[0]/2197.0 + 7296.0*h*s3[0]/2197.0,y[1] + 1932.0*h*s1[1]/2197.0 - 7200.0*h*s2[1]/2197.0 + 7296.0*h*s3[1]/2197.0)
        s5 = f(t + h, y[0] + 439.0*h*s1[0]/216.0 - 8.0*h*s2[0] + 3680.0*h*s3[0]/513.0 - 845.0*h*s4[0]/4104.0,y[1] + 439.0*h*s1[1]/216.0 - 8.0*h*s2[1] + 3680.0*h*s3[1]/513.0 - 845.0*h*s4[1]/4104.0)
        s6 = f(t + h/2.0, y[0] - 8.0*h*s1[0]/27.0 + 2*h*s2[0] - 3544.0*h*s3[0]/2565 + 1859.0*h*s4[0]/4104.0 - 11.0*h*s5[0]/40.0,y[1] - 8.0*h*s1[1]/27.0 + 2*h*s2[1] - 3544.0*h*s3[1]/2565 + 1859.0*h*s4[1]/4104.0 - 11.0*h*s5[1]/40.0)
        w1 = y[0] + h*(25.0*s1[0]/216.0 + 1408.0*s3[0]/2565.0 + 2197.0*s4[0]/4104.0 - s5[0]/5.0)
        w2 = y[1] + h*(25.0*s1[1]/216.0 + 1408.0*s3[1]/2565.0 + 2197.0*s4[1]/4104.0 - s5[1]/5.0)
        q1 = y[0] + h*(16.0*s1[0]/135.0 + 6656.0*s3[0]/12825.0 + 28561.0*s4[0]/56430.0 - 9.0*s5[0]/50.0 + 2.0*s6[0]/55.0)
        q2 = y[1] + h*(16.0*s1[1]/135.0 + 6656.0*s3[1]/12825.0 + 28561.0*s4[1]/56430.0 - 9.0*s5[1]/50.0 + 2.0*s6[1]/55.0)

        return w1,w2, q1,q2
t=0
T=10**5
poss=[-903482.12391302, -6896293.6960525, 0. ,a*(1-e),0,0 ]
sol = solve_ivp(accelerations2, [0, 10**5], poss,t_eval=np.linspace(0,10**5,1))
print(sol)

不确定该错误的含义，因为我尝试了许多不同的 t_evl，但似乎没有任何效果。

Answer 1

solve_ivp 中的默认值适用于“正常”情况，其中变量的比例与 0.1 到 100 的范围相差不大。您可以通过重新调整问题来实现这些比例，以便所有长度和相关常数以 AU 为单位，所有时间和相关常数以天为单位。

或者您可以尝试将绝对容差设置为合理的值，例如1e-4*AU。

正如我最近在关于这个主题的另一个问题中告诉你的那样，使用正确的一阶系统也有帮助。在机械系统中，您通常会得到二阶 ODE x''=a(x)。那么传递给 ODE 求解器的一阶系统是[x', v'] = [v, a(x)]，可以实现为

def firstorder(t,state):
    pos, vel = state.reshape(2,-1);
    return [*vel, *accelerations2(t,pos)]

接下来，将地球的加速度应用于地球以及将太阳的加速度应用于太阳总是有帮助的。即，固定对象的顺序。目前初始化首先是太阳，而在加速度计算中，您首先将状态视为地球。先全部切换到太阳

def accelerations2(t,pos):
    pos=pos.reshape(-1,3)
    # pos[0] = sun, pos[1] = earth
    norme=sum( (pos[1]-pos[0])**2 )**0.5
    gravit = G*(pos[1]-pos[0])/norme**3
    sunacc = me*gravit
    earthacc = -ms*gravit
    totacc=earthacc+sunacc
    return [*sunacc,*earthacc]

然后使用正确再现的自然常量就永远不会出错，例如

 G = 6.67E-11

然后求解器调用并打印格式为

state0=[*sunpos, *earthpos, *sunvel, *earthvel]
sol = solve_ivp(firstorder, [0, T], state0, first_step=1e+5, atol=1e-6*a)
print(sol.message)
for t, pos in zip(sol.t, sol.y[[0,1,3,4]].T): 
    print("%.6e"%t, ", ".join("%8.4g"%x for x in pos))

给出短表

The solver successfully reached the end of the integration interval.

       t         x_sun       y_sun    x_earth    y_earth
0.000000e+00 -9.035e+05, -6.896e+06,  7.5e+10,        0
1.000000e+05 -9.031e+05, -6.896e+06, 7.488e+10, 5.163e+09

也就是说，对于这一步，求解器只需要一个内部步骤。