<p>我认为对于您选择的参数,您遇到了<a href="http://en.wikipedia.org/wiki/Stiff_equation" rel="noreferrer">stiffness</a>的问题-由于数值不稳定,解算器的步长在解曲线的斜率实际上相当浅的区域变得非常小。由<code>scipy.integrate.odeint</code>包装的Fortran解算器<code>lsoda</code>试图在适合“stiff”和“non stiff”系统的方法之间自适应地切换,但在这种情况下,似乎无法切换到stiff方法。</p>
<p>非常粗略地说,你只需大幅度增加最大允许步数,解算器最终就会达到:</p>
<pre><code>SIR = spi.odeint(eq_system, PopIn, t_interval,mxstep=5000000)
</code></pre>
<p>更好的选择是使用面向对象的ODE解算器<code>scipy.integrate.ode</code>,它允许您显式地选择是使用stiff方法还是非stiff方法:</p>
<pre><code>import numpy as np
from pylab import *
import scipy.integrate as spi
def run():
#Parameter Values
S0 = 99.
I0 = 1.
R0 = 0.
PopIn= (S0, I0, R0)
beta= 0.50
gamma=1/10.
mu = 1/25550.
t_end = 15000.
t_start = 1.
t_step = 1.
t_interval = np.arange(t_start, t_end, t_step)
#Solving the differential equation. Solves over t for initial conditions PopIn
def eq_system(t,PopIn):
'''Defining SIR System of Equations'''
#Creating an array of equations
Eqs= np.zeros((3))
Eqs[0]= -beta * (PopIn[0]*PopIn[1]/(PopIn[0]+PopIn[1]+PopIn[2])) - mu*PopIn[0] + mu*(PopIn[0]+PopIn[1]+PopIn[2])
Eqs[1]= (beta * (PopIn[0]*PopIn[1]/(PopIn[0]+PopIn[1]+PopIn[2])) - gamma*PopIn[1] - mu*PopIn[1])
Eqs[2]= gamma*PopIn[1] - mu*PopIn[2]
return Eqs
ode = spi.ode(eq_system)
# BDF method suited to stiff systems of ODEs
ode.set_integrator('vode',nsteps=500,method='bdf')
ode.set_initial_value(PopIn,t_start)
ts = []
ys = []
while ode.successful() and ode.t < t_end:
ode.integrate(ode.t + t_step)
ts.append(ode.t)
ys.append(ode.y)
t = np.vstack(ts)
s,i,r = np.vstack(ys).T
fig,ax = subplots(1,1)
ax.hold(True)
ax.plot(t,s,label='Susceptible')
ax.plot(t,i,label='Infected')
ax.plot(t,r,label='Recovered')
ax.set_xlim(t_start,t_end)
ax.set_ylim(0,100)
ax.set_xlabel('Time')
ax.set_ylabel('Percent')
ax.legend(loc=0,fancybox=True)
return t,s,i,r,fig,ax
</code></pre>
<p>输出:</p>
<p><img src="https://i.stack.imgur.com/llAGh.png" alt="enter image description here"/></p>