Scipy odeint, дающий предупреждение lsoda

Question

Я совершенно не знаком с кодированием, и я хочу решить эти 5 дифференциальных уравнений численно. Я взял python template и применил его к моему делу. Вот упрощенная версия того, что я писал:

import numpy as np 
from math import * 
from matplotlib import rc, font_manager 
import matplotlib.pyplot as plt 
from scipy.integrate import odeint 

#Constants and parameters 
alpha=1/137. 
k=1.e-9  
T=40.  
V= 6.e-6 
r = 6.9673e12 
u = 1.51856e7 

#defining dy/dt's 
def f(y, t): 
     A = y[0] 
     B = y[1] 
     C = y[2] 
     D = y[3] 
     E = y[4] 
     # the model equations 
     f0 = 1.519e21*(-2*k/T*(k - (alpha/pi)*(B+V))*A) 
     f1 = (3*B**2 + 3*C**2 + 6*B*C + 2*pi**2*B*T + pi**2*T**2)**-1*(-f0*alpha/(3*pi**3) - 2*r*(B**3 + 3*B*C**2 + pi**2*T**2*B) - u*(D**3 - E**3)) 
     f2 = u*(D**3 - E**3)/(3*C**2) 
     f3 = -u*(D**3 - E**3)/(3*D**2) 
     f4 = u*(D**3 - E**3)/(3*E**2) + r*(B**3 + 3*B*C**2 + pi**2*T**2*B)/(3*E**2) 
     return [f0, f1, f2, f3, f4] 

# initial conditions 
A0 = 2.e13 
B0 = 0. 
C0 = 50.   
D0 = 50. 
E0 = C0/2. 

y0 = [A0, B0, C0, D0, E0]  # initial condition vector 
t = np.linspace(1e-15, 1e-10, 1000000) # time grid 

# solve the DEs 
soln = odeint(f, y0, t, mxstep = 5000) 
A = soln[:, 0] 
B = soln[:, 1] 
C = soln[:, 2] 
D = soln[:, 3] 
E = soln[:, 4] 

y2 = [A[-1], B[-1], C[-1], D[-1], E[-1]] 
t2 = np.linspace(1.e-10, 1.e-5, 1000000) 
soln2 = odeint(f, y2, t2, mxstep = 5000) 
A2 = soln2[:, 0] 
B2 = soln2[:, 1] 
C2 = soln2[:, 2] 
D2 = soln2[:, 3] 
E2 = soln2[:, 4] 

y3 = [A2[-1], B2[-1], C2[-1], D2[-1], E2[-1]]  
t3 = np.linspace(1.e-5, 1e1, 1000000) 
soln3 = odeint(f, y3, t3) 
A3 = soln3[:, 0] 
B3 = soln3[:, 1] 
C3 = soln3[:, 2] 
D3 = soln3[:, 3] 
E3 = soln3[:, 4] 

#Plot 
rc('text', usetex=True) 
plt.subplot(2, 1, 1) 
plt.semilogx(t, B, 'k') 
plt.semilogx(t2, B2, 'k') 
plt.semilogx(t3, B3, 'k') 

plt.subplot(2, 1, 2) 
plt.loglog(t, A, 'k') 
plt.loglog(t2, A2, 'k') 
plt.loglog(t3, A3, 'k') 

plt.show() 

Я получаю следующее сообщение об ошибке:

lsoda-- warning..internal t (=r1) and h (=r2) are   
    such that in the machine, t + h = t on the next step 
    (h = step size). solver will continue anyway   
    In above, R1 = 0.2135341098625E-06 R2 = 0.1236845248713E-22 

Для некоторого набора параметров, при воспроизведении с mxstep в odeint (также пытались hmin и hmax но Ждут» t замечают какую-либо разницу), хотя ошибка сохраняется, мои графики выглядят хорошо и не затрагиваются, но в большинстве случаев они есть. Иногда ошибка я получаю просит меня работать с опцией odeint full_output=1 и при этом я получаю:

A2 = soln2[:, 0] 
TypeError: tuple indices must be integers, not tuple 

Я не понимаю, что это означает, что при поиске для него.

Я хотел бы понять, где лежит проблема и как ее решить. Является ли odeint даже подходящим для того, что я пытаюсь сделать?

Answer 1

В этом lsoda предупреждение, t относится к текущему значению времени, а h относится к текущему размеру шага интеграции. Размер шага стал настолько близким к нулю, что текущее время плюс размер шага оценивается как равное текущему времени из-за ошибки округления (то есть r1 + r2 == r1). Такая проблема обычно возникает, когда ODEs, которые вы интегрируете, плохо себя ведет.

На моей машине предупреждающее сообщение появляется только при вычислении soln2. Здесь я построил значения каждого из параметров в области, где происходят предупреждения (обратите внимание, что я для ясности переключился на линейные оси). Я пометил шаг по времени, где впервые появилась ошибка lsoda (в r1 = 0.2135341098625E-06) с красной линией:

Это не случайно, что появление предупреждающего сообщения совпадает с «изломом» в E!

Я подозреваю, что происходит то, что изгиб представляет особенность в градиенте E, что вынуждает решателя делать все меньше и меньше шагов, пока размер шага не достигнет пределов точности с плавающей запятой. Фактически, вы можете видеть другую точку перегиба в D, которая совпадает с «колебанием» в B, вероятно, вызванным тем же явлением.

Для получения общих советов о том, как иметь дело с особенностями при интеграции ОДУ, взгляните на section 5.1.2 here (или на любой достойный учебник по ОДУ).

---

Вы получали ошибку с full_output=True, потому что в этом случае odeint возвращает кортеж, содержащий раствор и dict содержащий дополнительную выходную информацию. Попробуйте распаковать кортеж следующим образом:

soln, info = odeint(..., full_output=True)