Ошибка solve_ivp: 'Требуемый размер шага меньше, чем интервал между числами.' - PullRequest
Новая
       22

Ошибка solve_ivp: 'Требуемый размер шага меньше, чем интервал между числами.'

1 голос
/ 07 января 2020

Я пытался решить ньютоновскую проблему двух тел, используя RK45 от scipy, но продолжаю сталкиваться с ошибкой 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 и, похоже, ничего не работает.

1 Ответ

2 голосов
/ 08 января 2020

Значения по умолчанию в solve_ivp сделаны для «нормальной» ситуации, когда шкалы переменных не слишком отличаются от диапазона от 0,1 до 100. Вы можете достичь этих шкал, изменив масштаб задачи так, чтобы все длины и связанные константы находятся в AU, а все время и связанные константы - в днях.

Или вы можете попытаться установить абсолютное допуск на что-то разумное, например 1e-4*AU.

. Это также помогает использовать правильная система первого заказа, как я недавно говорил вам в другом вопросе по этой теме c. В механической системе вы обычно получаете 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

, то есть для этого шага решателю требуется только один внутренний шаг.

...