SIR模型的扰动分析

SIR 模型

模型

SIR 模型把人群分成三类：

$S(t)$ : 易感人群；
$I(t)$ : 感染人群；
$R(t)$ : 痊愈人群。

$[0,1]$ 之间。

再引入两个参数：

$\beta$ : 感染率;
$\gamma$ : 痊愈率。

SIR 模型是下面这个微分方程组

\left\{\begin{array}{l} \frac{{\rm d} S}{{\rm d} t}= -\beta I S \\ \frac{{\rm d} I}{{\rm d} t}= \beta I S - \gamma I\\ \frac{{\rm d} R}{{\rm d} t}= \gamma I. \end{array}\right. \label{SIR}

$S(0) = S_0 \approx 1, I(0) = I_0$ $I_0 + S_0 + R(0) = 1$ $R_0 = \beta/\gamma > 1$ $R$ 的初始值.

$IS$ $\beta$ $\gamma I$ 是痊愈者增加的速度。

$\frac{{\rm d}}{{\rm d} t}(S(t) + I(t) + R(t)) = 0$ $S(t) + I(t) + R(t) = 1$ $R(t)$ 中。

$\eqref{SIR}$ $S$ $R$ $[0,1]$ $t\to \infty$ $S_{\infty}, R_{\infty}$ .

(S, R) 的分析

$t>0$ $S(t) + I(t) + R(t) = 1$ $I$ $S, R$ 的方程组：

\left\{\begin{array}{l}<br>\frac{{\rm d} S}{{\rm d} t}= - \beta(1-S-R) S \\<br>\frac{{\rm d} R}{{\rm d} t}= \gamma (1-S-R).<br>\end{array}\right. \label{SR}

$\eqref{SR}$ $t\to \infty$ $S_{\infty}+ R_{\infty} = 1$ $I_{\infty} = \lim_{t\to \infty} I(t) = 0$ .

$\eqref{SR}$ . 形式上两个方程相除就得到了

\frac{{\rm d}S}{{\rm d} R} = - R_0 S

$S(0) = S_0， R(0)=0$ . 那么通过分离变量法，我们得到

S(t) = S_0 e^{-R_0 R(t)}, \quad R(t) = \frac{1}{R_0}\ln \frac{S_0}{S(t)}.\label{SRrelation}

$R$ $R(0)\ll 1$ ，公式仍然近似成立，

$S$ $R$ $\eqref{SR}$ $R$ 的一个非线性 ODE

\frac{{\rm d} R}{{\rm d} t}= \gamma (1- S_0 e^{- R_0 R} -R).

$t\to \infty$ $R_{\infty} = 1 - S_{\infty}$ $\eqref{SRrelation}$ ，我们得到方程

S_{\infty} = S_0 \exp (-R_0 (1 - S_{\infty})).

$I_0 \ll 1$ $S_0$ $1$ ，那么要求解的方程就简化为

S_{\infty} = e^{-R_0(1-S_{\infty})} \label{Sinf}

$\eqref{Sinf}$ $\eqref{Sinf}$ $S_{\infty}$ 的一个非常好的近似。

I 的分析

$I$ 的方程重写成

\frac{{\rm d} I}{{\rm d} t}= \beta \left (S - \frac{1}{R_0}\right ) I,

$S$ $(0,1)$ $R_0>1$ $S - 1/R_0$ 的符号把时间分成三个阶段.

$S > \frac{1}{R_0}$ $I' > 0$ $I$ 是在增长的。
$t^*$ $S(t^*) = 1/R_0$ .
$I' < 0$ $I$ 开始下降并趋近于零。

扰动分析

$S, I, R \in [0,1]$ $S_0$ 用 1 替代。要求解的非线性方程简化为

\frac{{\rm d} R}{{\rm d} t}= \gamma (1- e^{- R_0 R} -R) \label{nonlinearR}

$f(x) = 1 - e^{-R_0 x} - x$ $R'=\gamma I$ $I(t) =f(R(t))$ .

$\eqref{nonlinearR}$ 的解可以写成

\int \frac{{\rm d}\,R}{f(R)} = \gamma \, t. \label{solution}

$1/f(x)$ 的积分。

$f$ $f(x)$ $\eqref{solution}$ 的近似解。

$\bar t$ $\bar R\in (0,1)$ $\bar R = R(\bar t)$ 的解可以写成下面这种形式

\int_{\bar R}^{R} \frac{{\rm d}\,x}{f(x)} = \gamma \, (t - \bar t). \label{solution2}

$[\bar R, R]$ $f(x)$ $\int_0^1 1/x \, {\rm d}x$ .

线性逼近

$f(x) = a x + b$ 的形式，我们很容易求到线性逼近 ODE的解为

R_L(t) = \left ( \bar R + \frac{b}{a}\right )e^{a\gamma(t - \bar t)} - \frac{b}{a}.

$b=0$ $\bar R = 0$ . 通过求导，我们得到

I_L(t) = a\gamma \left ( \bar R + \frac{b}{a}\right )e^{a\gamma(t - \bar t)} = \bar Ie^{a\gamma(t - \bar t)}.

$\bar I$ $\bar R$ .

二次逼近

$f(x) = ax^2 + b x + c = a (x - \lambda_1)(x - \lambda_2)$ $\lambda_1\leq \lambda_2$ $f(x)$ $(\bar R, R) \subset [\lambda_1, \lambda_2]$ ，这样积分函数的分母就不为零。

使用分解

\frac{\lambda_1 - \lambda_2}{(x - \lambda_1)(x - \lambda_2)} = \frac{1}{x - \lambda_1} + \frac{1}{\lambda_2 - x}

$\eqref{solution2}$ 求积分而得到

\frac{R(t) - \lambda_1}{\lambda_2 - R(t)} = \frac{\bar R - \lambda_1}{\lambda_2 - \bar R} \exp(\rho (t - \bar t)), \label{Rfirst}

$\rho = a\gamma (\lambda_1 - \lambda_2)$ $\eqref{Rfirst}$ $R$ 的二次逼近

R_Q(t) = \frac{\lambda_2C\exp(\rho (t-\bar t)) + \lambda_1}{C\exp(\rho (t-\bar t)) + 1} = \lambda_2 + \frac{\lambda_2 - \lambda_1}{C\exp(\rho (t-\bar t)) + 1},

其中使用到的参数是

\rho = a\gamma (\lambda_1 - \lambda_2), \quad C =\frac{\bar R - \lambda_1}{\lambda_2 - \bar R}.

$R_Q$ $\lambda_1$ $\lambda_2$ ，因为

\lim_{t\to \infty} R_Q(t) = \lambda_2, \quad \lim_{t\to -\infty} R_Q(t) = \lambda_1.

$R_Q'(t^*) = 0$ 的公式是

t^* = \bar t + \frac{1}{\rho}\ln \left (\frac{\lambda_2 - \bar R}{\bar R - \lambda_1}\right ). \label{tbartstar}

$R' = \gamma I$ $I_Q$

I_Q(t) = \frac{(\lambda_2 - \lambda_1)\rho C \exp(\rho (t-\bar t))}{\gamma (C\exp(\rho (t-\bar t)) + 1)^2}.

泰勒展开

$f$ 和它的一阶、二阶导数如下

\begin{align*} f(x) &= 1 - e^{-R_0 x} - x, \\ f'(x) & = R_0e^{-R_0x} - 1, \\ f''(x) & = -R_0^2e^{-R_0x}. \end{align*}

$\bar x$ 处的展开是

f(x) = f(\bar x) + f'(\bar x)(x-\bar x) + \frac{1}{2}f''(\bar x)(x-\bar x)^2 + \mathcal O(|x-\bar x|^3).

$\bar x$ 做展开。然后用推导出的公式来求解相应的线性和二次ODE。我们将看到，使用简单的微积分和ODE知识，就可以得到函数非常好的近似。

初始阶段

$\bar t=0, \bar x = 0$ . 直接计算可得

f(0) = 0, \; f'(0) = R_0 - 1, \; f''(0) = -R_0^2.

$a=R_0-1$ $b=0$ , 所以我们得到

R_L(t) = \bar R e^{(\beta - \gamma) t}.

$R_L$ $y' = \beta y - \gamma y$ $\beta$ $\gamma$ 的速率指数下降。

$\bar R >0$ $I(t) = R'(t)/\gamma$ $\bar R$ $I$ $I_0$ 是给定的，那么

I_L(t) = I_0 e^{(\beta - \gamma) t},

$\bar R = I_0/(R_0-1)$ ，这个初值是很小的正数。

$f$ $x=0$ $(R_0 - 1) x - \frac{1}{2}R_0^2x^2 = x(R_0 - 1 - \frac{1}{2}R_0^2x)$ . 基于此，我们得到

a = -\frac{1}{2}R_0^2, \quad \lambda_1 = 0, \quad \lambda_2 = \frac{2(R_0-1)}{R_0^2}, \quad \rho = \beta - \gamma.

$a, \lambda_1,\lambda_2$ $\rho = \beta -\gamma$ .

earlyapproximationR

earlyapproximationI

我们把通过扰动分析得到的解析解和数值解来做比较。图中蓝色曲线是通过MATLAB 的 ode45 求解的数值解，红色曲线是线性ODE的解，黑色是二次ODE的解。

$0\leq t \leq 30$ $R_L$ $R_Q$ $I_Q$ 会指数下降到 0.

$\eqref{tbartstar}$ $\lambda_2$ $\bar x = 0$ 处的近似，之后的逼近就没有多大的可信度了。

结束阶段

$f(x)$ $\bar x = R_{\infty}$ 的值

f(R_{\infty}) = 0, \quad f'(R_{\infty}) = R_0\exp(-R_0 R_{\infty}) - 1, \quad f''(R_{\infty}) = -R_0^2\exp(-R_0 R_{\infty}).

$f(R_{\infty}) = I(\infty) = 0$ $f(x)$ $\exp(-R_0 R_{\infty}) = 1-R_{\infty} = S_{\infty}$ 从而可以如下简化

f(R_{\infty}) = 0, \quad f'(R_{\infty}) = R_0S_{\infty} - 1, \quad f''(R_{\infty}) = -R_0^2S_{\infty}.

$a = R_0S_{\infty} - 1, \, b/a = - R_{\infty}$ 所以

R_L(t) = R_{\infty} + (\bar R - R_{\infty}) \exp(\rho (\bar t- t))

I_L(t) = (R_{\infty}-\bar R) \rho \exp(\rho (\bar t-t)).

$\rho = \gamma - \beta S_{\infty}$ $\beta - \gamma$ 要大于结束阶段归零的速率。也就是说，一开始，感染人数会以更快的速率指数增长，后面降下来的时候速率要慢一些，但还是指数型的。

$\bar R = R_{\infty}$ $I$ $\bar t, \bar R = R(\bar t)$ , 但不能太靠近已经平缓的尾部.

$f'(R_{\infty})(x - R_{\infty}) + \frac{1}{2} f''(R_{\infty})(x-R_{\infty})^2 = \frac{1}{2} f''(R_{\infty})(x-R_{\infty})\left [2f'(R_{\infty})/f''(R_{\infty})+ x - R_{\infty}\right ].$

直接计算得到

\begin{align*} a &= \frac{1}{2} f''(R_{\infty}) = -\frac{1}{2}R_0^2S_{\infty}, \\ \lambda_1 &= R_{\infty} - 2f'(R_{\infty})/f''(R_{\infty}) = R_{\infty} - \frac{2(1-R_0S_{\infty})}{R_0^2S_{\infty}},\\ \lambda_2 &= R_{\infty},\\ \rho &= \gamma - \beta S_{\infty}. \end{align*}

代入公式我们就得到了二次逼近ODE的解。

endingstateR

endingstageI

$R_Q, I_Q$ $R_L, I_L$ 更好。但是当时间在初始阶段时，两者都差了很远，实际上在拐点附近误差已经很大了。

拐点阶段

$\bar t = t^*, \bar x = \bar R = R^*$ $t^*$ $I$ $S^* = S(t^*), I^* = I(t^*), R^* = R(t^*).$

$I$ 的方程

\frac{{\rm d} I}{{\rm d} t}= \beta \left (S - \frac{1}{R_0}\right ) I,

$I'(t^*) = 0$ 得到

S^* = \frac{1}{R_0},\quad R^* = \frac{\ln R_0}{R_0}, \quad I^* = 1 - \frac{1+\ln R_0}{R_0}.

$I = f(R)$ , 我们得到

f(R^*) = I^*, \; f'(R^*) = 0, \; f''(R^*) = - R_0.

$I_L = I^*$ 只是常数逼近，

R_L(t) = R^* + I^*(t-t^*).

$f$ 的二次逼近是

f_Q(x) = I^* - \frac{1}{2}R_0(x - R^*)^2,

根据此可以算出

a = -\frac{1}{2}R_0, \quad \lambda_1 = R^*-\sqrt{\frac{2I^*}{R_0}},\quad \lambda_2 = R^*+\sqrt{\frac{2I^*}{R_0}}, \quad \rho = \sqrt{2(R_0 - \ln R_0 - 1)}.

根据对称性，函数

R_Q(t) = \lambda_2 + \frac{\lambda_1 - \lambda_2}{\exp(\rho(t-t^*))+ 1},

I_Q(t) = \frac{(\lambda_2 - \lambda_1)\rho \exp(\rho(t-t^*))}{\gamma(\exp(\rho(t- t^*))+ 1)^2}

peaktimeR

peaktimeI

$R_0$ $t^*$ .

全局逼近

$R_Q, I_Q$ 已经提供了很好的逼近，但在时间的两端仍然不是太精确。很自然地，我们可以把三段时间的二次逼近联合起来得到整体更好的逼近。

$t^*$ $t^*$ $[t^* - t_1, t^* + t_2]$ $R_{Q}^1,R_{Q}^2,R_{Q}^3$ , 我们可以使用下面的分片函数定义

R_g = \begin{cases} R_Q^1(t) & 0 \leq t < t^* - t_1\\ R_Q^2(t) & t^* - t_1 \leq t \leq t^* - t_2\\ R_Q^3(t) & t > t^* - t_2. \end{cases}

$I_g$ .

$R_g, I_g$ $\beta,\gamma$ $t^*$ $\eqref{nonlinearR}$ ，但利用简单的微积分和ODE知识，就能得到非常好的逼近。

RIglobal