Scipy odeint给lsoda警告 [英] Scipy odeint giving lsoda warning

问题描述

我是编码的新手，我想用数值方法求解这5个微分方程.我使用了 python模板并将其应用于我的案例.这是我写的内容的简化版本:

I am totally new to coding and I want to solve these 5 differential equations numerically. I took a python template and applied it to my case. Here's the simplified version of what I wrote:

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

对于某些参数集，在odeint中使用mxstep时(也尝试过hmin和hmax，但没有发现任何差异)，尽管错误仍然存​​在，但我的图形看起来不错并且没有受到影响，但大多数时候都是这样. 有时，我得到的错误要求我使用odeint选项full_output=1运行，并这样做:

For some set of parameters, upon playing with mxstep in odeint (also tried hmin and hmax but didn't notice any difference), although the error persists, my graphs look good and are not impacted, but most of the times they are. Sometimes the error I get asks me to run with the odeint option full_output=1 and in doing so I obtain:

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

搜索时我不明白这是什么意思.

I don't get what this means when searching for it.

我想了解问题出在哪里以及如何解决. odeint甚至适合我想要做的事情吗?

I would like to understand where the problem lies and how to solve it. Is odeint even suitable for what I'm trying to do?

推荐答案

在lsoda警告中，t表示当前时间值，而h表示当前积分步长.步长变得非常接近零，由于舍入误差(即r1 + r2 == r1)，当前时间加上步长的值等于当前时间.当您集成的ODE表现不佳时，通常会发生这种问题.

In that lsoda warning, t refers to the current time value and h refers to the current integration step size. The step size has become so close to zero that the current time plus the step size evaluates as equal to the current time due to rounding error (i.e. r1 + r2 == r1). This sort of problem usually occurs when the ODEs you are integrating are badly behaved.

在我的机器上，警告消息似乎仅在计算soln2时出现.在这里，我绘制了发生警告的区域中每个参数的值(请注意，为清楚起见，我已切换到线性轴).我已经用红线标记了第一次出现lsoda错误的时间步长(在r1 = 0.2135341098625E-06):

On my machine the warning message only seems to occur when computing soln2. Here I've plotted the values of each of the parameters in the region where the warnings are occurring (note that I've switched to linear axes for clarity). I've marked the time step where the lsoda error first appeared (at r1 = 0.2135341098625E-06) with a red line:

警告消息的出现与E中的扭结"相吻合并非偶然！

It's no coincidence that the appearance of the warning message coincides with the 'kink' in E!

我怀疑正在发生的情况是，扭结表示E梯度中的奇异点，这迫使求解器采取越来越小的步长，直到步长达到浮点精度的极限为止.实际上，您可以在D中看到另一个拐点，该拐点与B中的摆动"重合，这可能是由相同的现象引起的.

I suspect that what's happening is that the kink represents a singularity in the gradient of E, which is forcing the solver to take smaller and smaller steps until the step size reaches the limits of floating point precision. In fact, you can see another inflection point in D which coincides with a 'wobble' in B, probably caused by the same phenomenon.

有关在集成ODE时如何处理奇点的一些一般性建议，请参阅第5.1.2节(或有关ODE的任何体面教科书).

For some general advice on how to deal with singularities when integrating ODEs, take a look at section 5.1.2 here (or any decent textbook on ODEs).

您遇到了full_output=True错误，因为在这种情况下，odeint返回一个包含解决方案的元组和一个dict包含其他输出信息的元组.尝试像这样打开元组的包装:

You were getting an error with full_output=True because in this case odeint returns a tuple containing the solution and a dict containing additional output information. Try unpacking the tuple like this:

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

这篇关于Scipy odeint给lsoda警告的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！