10

私は非常に単純な Susceptible-Infected-Recovered モデルを、アイドル状態のサイド プロジェクト (通常はかなり些細な作業) のために安定した人口で実装しています。しかし、PysCeS または SciPy のいずれかを使用してソルバー エラーが発生しています。どちらも、基礎となるソルバーとして lsoda を使用しています。これは、パラメータの特定の値に対してのみ発生し、その理由について困惑しています。私が使用しているコードは次のとおりです。

import numpy as np
from pylab import *
import scipy.integrate as spi

#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(PopIn,t):
    '''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

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

これにより、次のエラーが発生します。

 lsoda--  at current t (=r1), mxstep (=i1) steps   
       taken on this call before reaching tout     
      In above message,  I1 =       500
      In above message,  R1 =  0.7818108252072E+04
Excess work done on this call (perhaps wrong Dfun type).
Run with full_output = 1 to get quantitative information.

通常、私がそのような問題に遭遇した場合、私が設定した方程式系に最終的に何か問題がありますが、私はどちらもそれに問題があるとは思いません. 奇妙なことに、 mu を のように変更しても機能します1/15550。システムに何か問題があった場合に備えて、次のようにモデルを R に実装しました。

require(deSolve)

sir.model <- function (t, x, params) {
  S <- x[1]
  I <- x[2]
  R <- x[3]
  with (
    as.list(params),
{
    dS <- -beta*S*I/(S+I+R) - mu*S + mu*(S+I+R)
    dI <- beta*S*I/(S+I+R) - gamma*I - mu*I
    dR <- gamma*I - mu*R
  res <- c(dS,dI,dR)
  list(res)
}
  )
}

times <- seq(0,15000,by=1)
params <- c(
 beta <- 0.50,
 gamma <- 1/10,
 mu <- 1/25550
)

xstart <- c(S = 99, I = 1, R= 0)

out <- as.data.frame(lsoda(xstart,times,sir.model,params))

これもlsoda を使用していますが、問題なく動作しているようです。Python コードで何が問題になっているのか、誰にもわかりますか?

4

2 に答える 2

23

選択したパラメータでは、剛性の問題が発生していると思います。数値が不安定なため、解曲線の勾配が実際には非常に浅い領域では、ソルバーのステップ サイズが非常に小さくなっています。lsodaによってラップされているFortran ソルバーは、 scipy.integrate.odeint「堅い」システムと「堅くない」システムに適した方法を適応的に切り替えようとしますが、この場合、堅い方法への切り替えに失敗しているようです。

非常に大雑把に、最大許容ステップを大幅に増やすだけで、ソルバーは最終的にそこに到達します。

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

より良いオプションは、オブジェクト指向の ODE ソルバーscipy.integrate.odeを使用することです。これにより、スティッフ メソッドを使用するか非スティッフ メソッドを使用するかを明示的に選択できます。

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

出力：

ここに画像の説明を入力

于 2013-06-06T23:58:58.047 に答える
1

感染人口PopIn[1]はゼロに減少します。明らかに、(通常の) 数値の不正確さにより、PopIn[1]t=322.9 付近で負 (約 -3.549e-12) になります。その後、最終的に解は t=7818.093 付近で爆発し、PopIn[0]+無限大に向かって、-無限大にPopIn[1]向かっていきます。

編集：「クイックフィックス」の以前の提案を削除しました。それは疑わしいハックでした。

于 2013-06-06T23:59:23.510 に答える