1
具有时滞的传染病动力学模型
数值仿真
专业：信息与计算科学
班级：信息 042
姓名：张婧怡
学号：04052225
指导教师：牛大田老师
2
摘 要：:
具有时滞的传染病模型能较好反映传染病的潜伏期、免疫期等问题，对其研究越
来越受到重视。利用计算机模拟方法对具有时滞的传染病动力学模型进行了分析，
阐述了数值仿真的基本原则。采用时滞微分方程的数值解法对模型进行了定性分
析，给出了数值仿真的应用实例。结果表明该方法是有效的，并具有潜在的应用
价值。
关键词:
传染病模型；时滞微分方程；数值仿真
一、问题背景
近年来随着环境污染、生态破坏和频繁的国际交流，国内外重人传染病爆发
时有发生，如 SARS,禽流感和艾滋病等，这使得对传染病的研究越来越重要。在
对传染病的诸多研究中，利用数学模型对传染病进行定性研究是个重要课题。由
于具有时滞的传染病模型能更好与实际情况相符，所以近年来对其研究受到了人
们的广泛重视。随着计算机的发展，在研究方法上，除了传统的理论分析外，计
算机模拟也是研究的重要手段，一此重大发现就是通过数值仿真得到的。虽然如
此，但迄今为止利用计算机仿真方法研究时滞传染病模型的工作还很少。我们试
图在此方面对此尝试，将时滞微分方程数值算法应用于传染病模型研究，取得了
较好的结果，总结出了数值仿真的基本原则，说明该方法其有潜在的应用价值。
二、模型建立
“时滞”在传染病中是个基本的因素，并在传染病的传播过程中起着重要
作用，它可以反映传染病的潜伏期，患者对疾病的感染期和恢复者对玖病的免疫
期等实际现象，因此使用时滞”模型更贴近实际。如 Busenberg 和 Cooke 将时滞
因素引入到由媒介传播疾病的 SEIS 模型中用时滞项来反映传染病的潜伏期，建
立了如图 1 所小的仓室框图。
3
1
g
此模型中，把传染病地区的人群分为三类：用 S(t), E (t),I(t}分别表示易感者、在
潜伏期的感染者和染病者。箭头所指方向可以清楚的显小出各类人群流动的情况，
T>0 是模型的时滞项，代表疾病在人群中的潜伏期，r> 0 表小感染者被治愈后返
回到易感人群中的速率， 是易感者和传病媒介间的有效接触系数。由仓库框图
容易得到对应数学表达的动力学模型：
S
( ) ( ) ( )
( ) ( ) ( )
( ) ( ) ( )
d
I t S t I t
dt
dE
S t S t T I t T
dt
dI
S t T I t T I t
dt
g b
b b
b g
ì
= -
ï
ï
ï
= - - -
í
ï
ï
= - - -
ï
î
（1）
上述传染病动力学模型实质上是个具有时滞的微分方程组，对该模型进行数值
仿真，就是对方程组(1)求解，通过研究该方程组解的变化，从而得到如传染
病的发展趋势等相关内容。
三、模型分析
通过对实际模型的多次数值仿真实验，得到了数值仿真的基本原则
1、 运算精度优先：
由于一般传染病模型对应的方程组维数不高，对非刚性的时滞传染病模型
应使用至少四阶精度的数值方法如果采用精度不高的数值方法，会造成较大的误
差，从而对理论分析造成很大的偏差。下面我们从实际的程序数值运行和图示方
法加以比较。
2、模型分析程序：
一阶 Euler 方法
#include<stdio.h>
#define N 10
4
float f(float a,float b)
{float c;
c=b-2*a/b;
return c;
}
main()
{int i;
float h=0.2;
float X[N]={0.0};
float Y[N]={1.0};
for(i=1;i<=5;i++)
{ X[i]=X[i-1]+h;
Y[i]=Y[i-1]+h*f(X[i-1],Y[i-1] );
}
for(i=0;i<6;i++)
printf("%f\n",Y[i]);
}
四阶 R –K 方法
#include<stdio.h>
#define N 10
#define M 30
float f(float a,float b)
{float c;
c=b-2*a/b;
return c; }
main()
{int i;
float h=0.2;
float X[N]={0.0};
float Y[N]={1.0};
5
float k[M];
for(i=1;i<=5;i++)
X[i]=X[i-1]+h;
for(i=0;i<5;i++)
{k[i]=f(X[i],Y[i]);
k[i+1]=Y[i]+h/k[i]/2.0-2.0*(X[i]+h/2.0)/(Y[i]+h*k[i]/2.0);
k[i+2]=Y[i]+h*k[i+1]/2.0-2.0*(X[i]+h/2.0)/(Y[i]+h*k[i+1]/2.0);
k[i+3]=Y[i]+h*k[i]-2.0*(X[i+1])/(Y[i]+h*k[i]);
Y[i+1]=Y[i]+h*(k[i]+2*k[i+1]+2*k[i+2]+k[i+3])/6.0;
}
for(i=0;i<6;i++)
printf("%f\n",Y[i]);
}
3、图例解释：
某模型，采用四阶 R -K 方法可得到如图 2 所示的闭合轨线图，其中闭轨反
映了系统解的周期变化如果用一阶 Euler 方法得到如图 3 所示的相图，
