僵尸入侵模型

介绍

import IPython
IPython.display.YouTubeVideo("t1__PGs49jU")

在这个例子中，我们将使用课程中介绍的一种数值方法来求解一阶微分方程组。以下模型用于模拟僵尸入侵：

\[\begin{split}\Left\{\Begin{array}{ll} \displaystyle\frac{\mathm{d}s}{\mathm{d}t}&=\Pi-\displaystyle\beta S Z-\Delta S，\\ [4mm] \displaystyle\frac{\mathm{d}Z}{\mathm{d}t}&=\displaystyle\beta S Z+\zeta R-\alpha Z S，\\ [4mm] \displaystyle\frac{\mathrm{d}R}{\mathrm{d}t}&=\δS+\Alpha\DisplayStyle Z S-\Zeta R \end{array}\右。\end{split}\]

其中这种划分模型包括三种类型的个体 : susceptible individuals \(S\) ，僵尸 \(Z\) ，并删除了个人 \(R\) 不是死于僵尸袭击就是自然死亡。

\(\Pi\) ：出生率(假设不变)
\(\delta\) ：自然死亡率
\(\zeta\) ：死后变成僵尸的复活率
\(\beta\) : rate of becoming a zombie after contact with a zombie \(Z\)
\(\alpha\) ：僵尸的总破坏率

有关更多详细信息，请参阅其他资源：

When Zombies attack !

这是一个比前面看到的SIR模型更复杂的模型

模型图

下图总结了我们研究的僵尸入侵模型。

from IPython.core.display import SVG
import svgutils.transform as sg
from IPython.display import SVG,display
#fig = sg.SVGFigure("20cm", "10cm")
#fig1 = sg.fromfile('zombieDiagram.svg')
#plot1 = fig1.getroot()
#fig.append([plot1])
#fig.save("zombieDiagram2.svg")
#display(SVG(filename='zombieDiagram2.svg'))

display(SVG(filename='zombieDiagram2.svg'))

%matplotlib nbagg
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import integrate

世界末日僵尸的数值

Pi =  0. #birth rate
delta = 1.e-4 #natural death rate
beta = 9.5e-3 #contact rate
zeta = 1.e-4 #resurrection rate
alpha = 1.e-4 #destruction rate

S0 = 1000 #initial population
Z0, R0 = 0, 0
X0 = S0, Z0, R0
tmax = 10
Nt = 160
t = np.linspace(0, tmax, Nt+1) #time grid

第1部分：问题的重新表述

这个问题可以重新制定以符合标准的配方。 \(\dot X = f(X, t)\) ：

\[\begin{split}X=\BEGIN{pMatrix}S\\ Z\\ r\end{pMatrix}\end{split}\]

\[\begin{split}\四边形\点{X}= \BEGIN{pMatrix}\PI-\DisplayStyle\beta S Z-\Delta S\\ \displaystyle\displaystyle\beta S Z+\zeta R-\alpha Z S\\ \增量S+\Alpha\DisplayStyle Z S-\Zeta R \end{pMatrix} =F\Left(X\Right)\end{split}\]

编写函数 \(f\) 在Python中：

def derivative(X, t, Pi, delta, beta, zeta, alpha):
    S, Z, R = X
    dotS = Pi - beta * S * Z - delta * S
    dotZ = beta * S * Z + zeta * R - alpha * Z * S
    dotR = delta * S + alpha * Z * S - zeta * R
    return np.array([dotS, dotZ, dotR])

第2部分：解算

使用解ODE odeint 然后画出来。

X = integrate.odeint(derivative, X0, t, args = (Pi, delta, beta, zeta, alpha))
S, Z, R = X.T #X.T order 3 x (Nt+1)

plt.figure()
plt.grid()
plt.title("Zombie model")
plt.plot(t, S, 'orange', label='Susceptible')
plt.plot(t, Z, 'r', label='Zombie')
plt.plot(t, R, 'g', label='Removed')
plt.xlabel('Time t, [days]')
plt.ylabel('Population')
plt.legend()

plt.show();

<IPython.core.display.Javascript object>

使用不同初始条件和不同数值进行测试：

\(R(0)\) is 1% of the initial population size \(S(0)\)
\(\Pi\) =8(每天8名新生婴儿)

R0 = 0.01 * S0
X0 = S0, Z0, R0
Pi = 8
X = integrate.odeint(derivative, X0, t, args =(Pi, delta, beta, zeta, alpha))
S, Z, R = X.T #X.T order 3 x (Nt+1)

plt.figure()
plt.grid()
plt.title("Zombie model")
plt.plot(t, S, 'orange', label='Susceptible')
plt.plot(t, Z, 'r', label='Zombie')
plt.plot(t, R, 'g', label='Removed')
plt.xlabel('Time t, [days]')
plt.ylabel('Population')
plt.legend()

plt.show();

<IPython.core.display.Javascript object>

现在，用显式Euler和RK4方法求解。比较方法。

def Euler(func, X0, t, Pi, delta, beta, zeta, alpha):
    """
    Euler solver.
    """
    dt = t[1] - t[0]
    nt = len(t)
    X  = np.zeros([nt, len(X0)])
    X[0] = X0
    for i in range(nt-1):
        X[i+1] = X[i] + func(X[i], t[i], Pi, delta, beta, zeta, alpha) * dt
    return X

def RK4(func, X0, t, Pi, delta, beta, zeta, alpha):
    """
    Runge Kutta 4 solver.
    """
    dt = t[1] - t[0]
    nt = len(t)
    X  = np.zeros([nt, len(X0)])
    X[0] = X0
    for i in range(nt-1):
        k1 = func(X[i], t[i], Pi, delta, beta, zeta, alpha)
        k2 = func(X[i] + dt/2. * k1, t[i] + dt/2., Pi, delta, beta, zeta, alpha)
        k3 = func(X[i] + dt/2. * k2, t[i] + dt/2., Pi, delta, beta, zeta, alpha)
        k4 = func(X[i] + dt    * k3, t[i] + dt, Pi, delta, beta, zeta, alpha)
        X[i+1] = X[i] + dt / 6. * (k1 + 2. * k2 + 2. * k3 + k4)
    return X

Xe = Euler(derivative, X0, t, Pi, delta, beta, zeta, alpha)

plt.figure()
plt.grid()
plt.title("Zombie model with explicit Euler")
plt.plot(t, Xe[:,0], 'orange', label='Susceptible')
plt.plot(t, Xe[:,1], 'r', label='Zombie')
plt.plot(t, Xe[:,2], 'g', label='Removed')
plt.xlabel('Time t, [days]')
plt.ylabel('Population')
plt.legend()

plt.show();

<IPython.core.display.Javascript object>

Xrk4 = RK4(derivative, X0, t, Pi, delta, beta, zeta, alpha)

plt.figure()
plt.grid()
plt.title("Zombie model with RK4")
plt.plot(t, Xrk4[:,0], 'orange', label='Susceptible')
plt.plot(t, Xrk4[:,1], 'r', label='Zombie')
plt.plot(t, Xrk4[:,2], 'g', label='Removed')
plt.xlabel('Time t, [days]')
plt.ylabel('Population')
plt.legend()

plt.show();

<IPython.core.display.Javascript object>