Числовая проблема в функции odeint в зависимости от начальных значений

Question

Я попытался отладить мой код (который решает связанные ODE). Я нашел Функция odeint имеет определенные пределы принятия десятичных значений из-за этой проблемы (возможно) мои графики отличаются для разных начальных условий.

проблема заключается в строке 20: x = np.linspace(-30., -50., 5)
проблема заключается в строке 25: state0 = [gi,Li,60.]

Когда у меня есть начальное условие в х (-30) и в state0 (60) код работать должным образом:

Но когда я изменить начальные условия х (-25) и в state0 (50) код не работает:

Я решил использовать тот же код в Fortran, он работает для обоих пределов. Может ли кто-нибудь предложить, как решить эту проблему.

import numpy as np 
import matplotlib.pyplot as plt 
import scipy.integrate as si 
#function for solving differential equations: 
def func1(state, x): 
    g = state[0] 
    L = state[1] 
    u = state[2] 
    phi1 = 1.645 
    phi2 = 2.* 1.202 
    N = (-2.*g/(6.*np.pi + 5.*g))*(18./(1. - 2.*L) + 5.*np.log(1.- 2.*L) - phi1 + 6.) 
    B = (-(2. - N)*L) - ((g/np.pi)* (5.*np.log(1.-2.*L) - phi2 + (5.*N/4.))) 
    udx = -2.*(1.- ((3.* 0.5* L)/(3.* 0.330))) #5.1 
    gdx = (udx) * (2. + N)*g #5.1a 
    Ldx = (udx) *(B) #5.1b 
    print gdx, Ldx, udx 
    return gdx, Ldx, udx 
gt = 10.**(-60) 
Lt = (1.202*gt)/np.pi 
x = np.linspace(-30., -50., 500) #initial conditions 
for i in np.arange(len(x)): 
    gi= np.float64(gt * np.exp(-4.*x[i])) 
    Li = np.float64(Lt * np.cosh(4.*x[i])) 
    break 
state0 = [gi,Li,60.] #initial conditions 
G= si.odeint(func1, state0, x) 
a = np.transpose(G) # a[0]=g; a[1]=lambda (L); a[2]=u 
plt.plot(a[1],a[0]) 
plt.title('RGIII Tragectory') 
plt.xlabel('L ---->') 
plt.ylabel('g ----->') 
plt.show() 

Answer 1

Действительно, odeint решатель борется с этой ОДУ системы. По-видимому, движение, описываемое этим ОДУ, очень неоднородно по времени. Исправление должен сказать решатель использовать более жесткий контроль ошибок:

G = si.odeint(func1, state0, x, rtol=1e-10) 

Это дает спиральный участок для обоих ваших наборов параметров.

---

Вообще, когда odeint действует странно, первое, что нужно сделать, это попросить его infodict:

G, info = si.odeint(func1, state0, x, full_output=True) 

infodict имеет много информации в нем, но уже беглый взгляд на info['hu'] (с настройками, вызвавшими вашу проблему) показывает смешные значения используемого временного шага; это как если бы решатель решил, что это решение никуда не денется и покрывает запрошенный интервал за один раз; рисуя крошечный сегмент линии.