solve_ivp error: 'Required step size is less than spacing between numbers.'

Question

一直在尝试使用 scipy 中的 RK45 解决牛顿二体问题，但是将运行保留在 TypeError 中：'Required step size is less than spacing between numbers.' 我尝试了 [=16= 的不同值] 比下面的要好，但似乎没有任何效果。

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 中的默认值适用于 "normal" 情况，其中变量的比例与 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))

给出短table

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

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

solve_ivp error: 'Required step size is less than spacing between numbers.'

solve_ivp error: 'Required step size is less than spacing between numbers.'

physics

typeerror

runge-kutta