奇怪的SciPy ODE集成错误

奇怪的SciPy ODE集成错误

我认为，对于您选择的参数，您会遇到刚度问题-
由于数值的不稳定性，求解器的步长在求解曲线的斜率实际上很浅的区域变得越来越小。

lsoda

由包裹的Fortran求解器

scipy.integrate.odeint

尝试在适用于“刚性”和“非刚性”系统的方法之间进行自适应切换，但是在这种情况下，它似乎无法转换为刚性方法。

非常粗略地讲，您可以大幅增加允许的最大步数，而求解器将最终达到目标：

SIR = spi.odeint(eq_system, PopIn, t_interval,mxstep=5000000)

更好的选择是使用面向对象的ODE求解器

scipy.integrate.ode

，它使您可以显式选择是使用刚性方法还是非刚性方法：

import numpy as npfrom pylab import *import scipy.integrate as spidef 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