回答此问题可获得 20 贡献值,回答如果被采纳可获得 50 分。
<p>我正在为一个空闲的项目实现一个非常简单的易受感染恢复模型和一个稳定的种群-通常是一个非常琐碎的任务。但是我使用PysCeS或SciPy遇到解算器错误,这两种方法都使用lsoda作为底层解算器。这种情况只发生在参数的特定值上,我不明白为什么。我使用的代码如下:</p>
<pre><code>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)
</code></pre>
<p>这会产生以下错误:</p>
<pre><code> 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.
</code></pre>
<p>通常,当我遇到这样的问题时,我建立的方程组有点严重的问题,但我俩都看不出有什么问题。奇怪的是,如果您将mu更改为类似<code>1/15550</code>的值,它也可以工作。如果系统出了问题,我<em>也</em>用R实现了模型,如下所示:</p>
<pre><code>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))
</code></pre>
<p>这个<em>也</em>使用lsoda,但似乎进展顺利。有人能看到Python代码中出了什么问题吗?</p>