typeerror - 解决_ivp错误: 'Required step size is less than spacing between numbers.'

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

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