实验报告
1 实验题目
儿童麻疹流行蔓延的数学模型
2 实验问题陈述
试组建一个能描述儿童麻疹流行蔓延的数学模型,我们将考虑在接种疫苗成为有效的防疫手段之前的麻疹的流行。下表一给出了英国伦敦在1647年---1660年间每年麻疹病的死亡人数。
表一:伦敦每年麻疹死亡人数(1647---1660)
年代 47 48 49 50 51 52 53 54 55 56 57 58 59 60 人数 5 92 3 33 33 62 8 52 11 153 15 80 6 74
可以看出,它是以2年为周期的周期性流行。已知麻疹的潜伏期是0.5周,在这段时期内一个被感染的孩子表面看来是正常的,但却会传染给别人。过了这段时间后,患病的孩子一直被隔离到痊愈为止。痊愈后孩子是免疫的。假设每个感染者随机地与他人接触。证明你的模型有某种周期性质。如果不然,就修改你的模型。因为麻疹的流行肯定是趋于周期式地出现的。估计你组建的模型中的参数,以拟合0.5周的潜伏期及2年周期性流行的观测结果。判断估计出的参数是否实际。
3 实验目的通过表中数据,建立麻疹流行蔓延模型,以拟合0.5周的潜伏期及2年周期性流行的观测结果,判断估计出的参数是否实际。
4 实验内容
模型假设:
(1)除感病特征外,人群中的个体间没有差异,感病者与易感者的个体在人群中混合是均匀的。
(2)人群的数量足够大,只考虑传染过程的平均效应。
(3)易感者感病的机会与他接触感病者的机会成正比。
(4)疾病的传染率为常数。
(5)一般的麻疹爆发在几十天,我们不考虑在一次麻疹爆发时间内某地区的出生
人口和死亡人口,以及人口的迁入和迁出。
(6)感病痊愈者(即移出者)移出模型,而不再成为易感者人群中的成员。
变量说明:
S(t):易感者在人群中所占的比例
I(t):感病者在人群中所占的比例
R(t):移出者在人群中所占的比例
K:疾病的传染率
h:单位时间内痊愈的百分数
一个传染期内每个病人有效接触易感者的平均人数,成为接触数---
初始时刻
问题分析:
对麻疹流行蔓延的周期性质进行说明。
通过对SIR模型及麻疹流行的机理分析,在一次麻疹爆发以后绝大多数人体内具有了麻疹免疫抗体,因此绝大多数新生婴儿体内具有抗体,考虑到引起流行周期的原因是易感人群的积累,易感人群来源于新生儿因母体抗体逐渐消失而易感、既往没有患过麻疹的儿童和成人。
通常认为在自然感染状态下,这些易感者积累到一个以上出生队列时,就达到爆发的“临界”。若将同一年出生的人群组定义为一个出生队列,出生队列出现的周期性在一定程度上可以说明麻疹流行的周期性。
所以我们用积累一个出生队列的时间来表示一次爆发的临界。
模型建立
通过对问题的分析,模型可以表示为:
,其中
考虑到初始条件,可知上述三个方程是相容的,因此可以化简为:
由于方程组无法求出解析解,故可以在S-I的像平面上讨论解的性质,相轨线的定义域为:
由以上方程可知轨线的方程为:
,
其解为:
。
5 实验结果分析与讨论
由题目中表格给出的麻疹死亡人数与年份的对应关系,用MATLAB编程画出的曲线图如图一所示:
【图一】
由于对S(t)和I(t)的求1解非常困难,所以先用数值计算的方法来预估计S(t)和I(t)的一般变化规律。
在方程(1)中设k=1,h=0.3,I(0)=0.02,S(0)=0.98。
编写MATLAB程序并运行得到如图二,
【图二】
【注】图中蓝色曲线为I(t),即病人比例;绿色曲线为S(t),即健康人比例。
从时间流程图中可以看出,随着时间的增加,S(t)单调递减,I(t)在时达到峰值以后会随时间减小,当时,S(t)值很小,而I(t)=0。说明在一次麻疹疫情爆发以后绝大多数的人体内已经具有麻疹免疫抗体,被移出除传染系统。
c. 对结果的分析
参数中取h为0.3,则潜伏期为1/h=3.33天,约等于0.5周。同时我们可以看到把h取为0.3得到的曲线符合实际情况,说明潜伏期为麻疹病毒的潜伏期是0.5周是正确的观点。
另外,当时,,,说明在一次麻疹疫情爆发以后绝大多数的人体内已经具有麻疹免疫抗体,被移出除传染系统。
查资料知,绝大多数的婴儿在9个月时血内的母亲抗体已测不出
【图三】
【注】:
表示第一年的年初;
表示第一年年末;
表示第二年年初;
表示第二年年末;
在时刻出生的婴儿到时刻抗体消失,时刻出生的婴儿到时刻抗体消失。
易感人群从时刻开始积累,在时刻易感人群刚好积累一个出生队列,因此易感人群积累一个出生队列的时间为2年。
当易感人群积累到一个出生队列时,就是第二次麻疹爆发的“临界”,因此可以说麻疹流行的周期为2年。
综上所述,儿童麻疹流行蔓延的模型具有周期性,且以0.5周的潜伏期和2年的周期性流行。
故模型所得结果与题目要求是一致的。
6 实验程序(Matlab或者其它软件语言陈述)
用MATLAB编写程序如下:
(1) 画图一:
x=47:60;
y=[5 92 3 33 33 62 8 52 11 153 15 80 6 74];
plot(x,y,'rp-.')
xlabel('年份'); ylabel('');
(2) S(t)和I(t)的变化规律及画图三
M文件为:(chuanran.m)
function y=chuanran(t,x)
a=1;b=0.3;
y=[a*x(1)*x(2)-b*x(1);-a*x(1)*x(2)];
命令框中输入:
>> ts=0:50;
>> x0=[0.020,0.98];
>> [t,x]=ode45('chuanran',ts,x0);
>> plot(t,x(:,1),'r*-');
>> hold on
>> plot(t,x(:,2),'b.-');
>> legend('病人的比例','健康人的比例');
第 1 页 共 5 页