甲流模型 科技创新

附件： 我们的模型 模型一 1.1模型假设 1. 康复者不再被感染 2. 疾病只在人与人之间传然 3. 患病者治愈几率一定
4. 每个季节人被感染的几率不一样
5. 疫苗一段时间后才出现 1.2符号说明
易感染者占总人数比列s(t) 已感染者占总人数比列i(t)
病愈免疫的移出者占总人数比列r(t) 每个病人的日有效接触率为λ 日治愈率为μ
感染者的日接触率：×λ感染者的日增长数量
＝
感染群日累计数量易感群数量
初期无任何措施SI 模型 N
i
t
∂∂Nsi λ= ()()1s t i t +=
记（t=0）病人比列为0i ， 则
0(1),(0)i
i i i i t
λ∂=-=∂ 解为： 0
1()1
1(1)t
i t e i λ-=
+- (1.1)
对方程（1.1）为Logistic 模型，做出方程的图像如下：
SI 模型的Ｉ（t ）～t 曲线 图1 SI 模型的dI dt
~Ｉ曲线
由图形知： ⎪⎪⎭
⎫ ⎝⎛-=-11ln 01i t m λ （1.2）
当t 取t m 时di /dt 取最大值，说明这个时候是传染病高潮到来时刻,应该引起我们的重视。

从图像中还可以看出当日接触率λ下降的时候，图像中的tm 的位置将会右移增大。

说明控制日接触率可以延缓传染病爆发的到来。

另外当 t 趋近于∞时，i 趋近于1 。

说明经过无限长的时间，所以人都会得病，这显然与实际不符。

产生错误的原因在于没有考虑病人可以治愈这个因素。

说明这个模型仍然有改进的空间。

但是在开始的的时候我们还没有发现传染病的严重性，并没有采取预防措施，所以曲线开始的一小段可以认为和实际是接近的。

考虑病人可以治愈这个因素我们又建立了模型二，另外还要假设这种传染病无免疫性，也就是说，病人治愈成为健康人，健康人可再次被感染。

同时总人数
N 不变，病人和健康人的 比例分别为)(),(t s t i ，每个病人每天有效接触人数为λ, 且使接触的健康人致病，将λ称为日接触率。

病人每天治愈的比例为μ称为日治愈率。

仍然在一段时间t ∆新增的病例数为在这段时间内被传染的人数减去已治愈的人，我们建立数学方程：
被治愈者的变化量：()()Nir t t Nir t Ni t μ+∆-=∆ 被感染者的变化量：()()Ni t t Ni t Nis t Ni t λμ+∆-=∆-∆
i
1 健康者的变化量： ()()Ns t t Ns t Nis t λ+∆-=-∆ 由上述关系式得： 'r i μ= 'i is i λμ=-
,s is λ=-
()()()1r t i t s t ++=，t=0时， ()00i i =,()00s s =,(0)0r =,001i s +=
1ln s i s s μλ=-+
01ln s s s s t s μλλ⎛⎫
∂=--- ⎪∂⎝
⎭ （1.3）
取
/σλμ=，则σ为整个传染期内每个病人有效的接触平均人数，
0()1()[()(1)](0)di t i t i t dt i i λσ⎧=---⎪
⎨
⎪=⎩
（1.4） 根据方程（1.4）di/dt 于i 的图像：
1-1/
i 接着我们做出i~t 的图像，发现当σ取不同值的时候i~t 的图像不同。

如下图所示：
从图像中可以看出：
⎪⎩⎪⎨⎧≤>-=∞1,
01
,1
1)(σσσ
i （1.5）
模型求解的结果：
作出()i t t -曲线，为了方便绘画，取0.25,0.1λμ==， 即λ
σμ
=
=25>1，求出0t →时的()i t 的极限值： y =3/20/(1/4+599/4*exp(-3/20*t))；
limit('3/20/(1/4+599/4*exp(-3/20*t))',t,inf) ans=3/5
图2. 1,σ>()i t t -曲线
模型分析：由此可见接触数σ =1 是一个阈值，当1σ>时，()i t 的增减性取决与0i 的大小，而且当1σ≤时，病人比例()i t 越来越小，最终趋于0 模型二 2.1假设条件
(1)人群只分为健康群(Suscceptible)、病群(Infective)和移出群（Removed ）三类，称为SIR 模型; 时刻t 这三类人群只占群总数的比例分别记作S(t)、I(t)和R(t)。

(2) 病群的日接触率为λ，日移出率为μ，传染期接触数为：λ
σμ
=。

2.2模型求解 由假设1显然有
S(t)+I(t)+R(t)=1 (2.1) 根据2，于是(1.4)仍成立。

对于移出群而言应有 ω=dR
N NI dt
(2.2) 于是(2.2)应修正为
λλωω⎧=-⎪⎪⎪=-⎨⎪⎪=⎪⎩
;;;dS
IS dt dI
SI I dt dR
I dt (2.3)
再记初始时刻的健康群只和病群的比例分别是0S (0S >0)和0I (0I >0),(不妨设移出群只的初始值R(0)= 0R =0)。

注意到条件(2.13)，模型可简化为
λλω⎧=-=⎪⎪⎨
⎪=-=⎪⎩00
,(0);,(0);
dS
IS S S dt dI IS I I I dt (2.4)
方程(2.4)无法求解出解析的解S(t)和I(t)。

我们先做数值计算[3]。

数值计算：
在方程(2.4)中设λ＝1，0.3ω=，I(0)=0.02， S(0)＝0.98；
用MATHLAB 软件编程： Function y=ill(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.02,0.98];
[t,x]=ode45(`ill`,ts,x0);[t,x] Plot(t,x(:;1),t,x(:,2)),grid,pause Plot(x(:,2),x(:,1)),grid [8],
输出的简明计算结果列入表1，I(t) ～S(t)的图形见图3，图4是I(t)～t 的图形称之为相轨线，初值I(0)=0.02，S(0)=0.98相当于图4中的0p 点，随着t 的增加，(S,I)沿轨线自右向左引动。

由表1、图3可以看出，I(t)由初值增长值
至约t=7时达到最大值，然后减少，t→∞，I→0；S(t)则单调减少，t→∞，S→0.0398。

表3 S(t)，I(t)的数值计算结果
00.20.40.60.81
0.050.10.150.20.250.30.35s~i 图形
图3. s~i 比例图
10
20304050
00.10.20.30.40.50.60.70.80.91s ，i 比例图形
图4()()i t t s t t --与的曲线
为了分析I(t)，S(t)的一般变化规律，需要进行相轨线分析。

相轨线分析：
在数值计算和图形观察的基础上，利用相轨线讨论解 I(t)，S(t)的性质。

S →I 平面称为相平面，相轨线在相平面上的定义域(S,I)∈D 为
D={(S,I)|S ≥0, I≥0，S+I≤1} (2.5) 在方程(2.4)中消去dt 并注意到σ的定义域，可得：
σ==-=0
01
1,|S S
dI I I dt S
(2.6) 容易求出方程(2.18)的解为
I =（0S +0I ）-S+
1
σ
㏑
S
S (2.7) 在定义域D 内，(2.7)式表示的曲线即为相轨线，如图5所示。

其中箭头表示了随着时间t 的增加S(t)、R （t ）和I(t)的变化趋势。

分别记S I ∞∞∞，，R 为他们的极限值，下面根据(2.4)和(2.7)分别讨论当t→∞时，S(t)，I(t)，R(t)的变化情况。

(1)：无论初始条件00S I ，如何，感染群只终将在系统中消失，既有I ∞=0。

此时（0，0）是平衡点。

事实上，由(2.4)和
0dS
dt
≤，而S(t)≥0，故S ∞存在。

由0dR I dt ω=≥，且 R(t)≤1，故R ∞存在。

从而I ∞
存在。

若
I ∞=∂>0,则I ∞>2
∂。

对于充分大的t 有I(t)>2
∂。

从而对充分大的t 有
/2dR
dt
μ>∂，着将导致R ∞=∞，这与R ∞存在矛盾。

可得I ∞=0。

因此在相平面上，(2.4)的所有轨线最终将与S 轴相交。

(2)：在(2.7)中令t→∞，取极限，并另令
I ∞=0，则S ∞将是方程
（0S +0I ）-S+
1
σ
㏑
S
S =0 (2.8) 记
1
ρσ
=，在(0,ρ)内的单根。

在图5中S ∞是轨线与S 轴在(0,ρ)内交点的横
坐标。

（3）：对于不同的初始位置（00,s ρ），(2.7)给出了(2.4)的相轨线族，它可写成S+I －
ρ㏑S ＝C(0
,s ρ)的。

易证明，它是[0,1)上的单峰曲线，在*
S ρ=处达到最大值Ｉ(*S )＝1-ρ＋
ρ(
S ρ
),然后单调减少至0。

而S(t)则单调减少到S ∞(0S )。

如图5中由1P 出发的轨线；若0S ≤ρ，则I(t)单调减少至0，S(t)则单调减少到S ∞(0S )。

如图5中由2P 出发的轨线。

综合上面的分析可以看出，如果仅当感染群的比例I(t)有一段时间在增长时才认为发生了甲型H1N1流感的蔓延，那么ρ就是一个阈值。

只有当0S >ρ时才会发生甲型H1N1流感蔓延的现象。

如果设法提高模型的阈值ρ，使的0S ≤ρ，就可以控制甲型H1N1流感的蔓延。

我们注意到在ρ＝μ/λ中，防控水平越高，日接触率λ越小；防疫水平越高，日治愈率μ越大，于是ρ 越大，提高防控水平和防疫水平有助于控制甲型H1N1流感的蔓延。

我们看到在SIR 模型中，参数ρ是重要的，通常称之为相移出率。

ρ可以由观测数据给出估计，因为病群比例的初始值0I 通常很小，可以得到
00
S S LnS LnS ρ∞∞-=- (2.9)
当甲型H1N1流感流行结束后得到S ∞和0S ,由(2.9)就可以给出ρ的估计。

模型改进：
将某个地区整个群总数仍旧为N ，在模型2的基础上又多划分一个潜伏群E （t ）；因为甲型H1N1流感易感禽只由于自身免疫力不同，当易感群和感染群进行有效接触时，一部分直接成为感染群，另一部分经过潜伏期，然后转化为感染群；但在潜伏期只是不发病，却又具有感染健康群的能力。

四个群之间关系如下图
那么模型应进一步修正为
λααθλθωω⎧=--⎪⎪⎪=-⎪⎪⎪=+-⎨
⎪⎪=⎪⎪⎪
====⎪⎩
0000.;;;;(0),(0),(0),(0)dS
IS IS dt dE
IS E dt dI
IS E I dt dR
I dt S S E E I I R R （2.10） 求解时用以下公式代具体的参数解（不详解） 潜伏期的日接触率：×α=
潜伏期的日增长数量
感染群日累计数量易感群数量
感染者的日接触率：×λ感染者的日增长数量
＝
感染群日累计数量易感群数量
退出率＝死亡率+治愈率：θ每日死亡数＋每日治愈数
＝
日累计确诊病例数
潜伏群的发病率：μ每日新确诊病例
＝
每日疑似累计数
其中， α：潜伏日接触率； λ：感染者日接触率；
θ：潜伏群的发病率；μ：传染群的退出率（包括死亡和病愈）。

模型三
季节因素使感染率增加为a
政府防控使感染率减少为b,疫苗使感染率减少d ，个人采取的保护措施使感染率减少为h
医院措施使治愈率增加为c 求解时用以下公式代具体的参数解（不详解）
感染者的日接触率：×λ感染者的日增长数量
＝感染群日累计数量易感群数量
()'1r c i μ=+
()()()
'11'1i is a b d h i c s is a b d h λμλ=+----+=-+---
()()()1r t i t s t ++=，t=0时，()00i i =,()00s s =,(0)0r =,001i s +=
()()00μμλλ取=1+c ,=1+a+b-d-h 模型的分析与变化将于模型二相同，所以直接把模型二中的λ，μ带入即可，
1ln s
i s s μλ=-+
01ln s s s s t s μλλ⎛⎫∂=--- ⎪∂⎝
⎭ 根据我们搜集到数据，首先用MATLAB 画出2009年7月15日到2010年3月31 日，一段时间段内H1N1在我国的发病例，我们分12个阶段
作图：
34
567
8910
4
4567891011
0.511.522.533.544.5
根据得出的图形，可以看出H1N1全球蔓延趋势如下：
在H1N1流行的初期，就是在开始统计的第一二个月时间段病例总数增长并不快，甚至在一定阶段显现出缓慢增加的形式，这是由于初期病例数量小造成的慢速的传染现象。

两月后，总病例数开始急剧上升，增长趋势类似于指数型爆发增长，这有可能是由于患病人数增加增大了传染几率，短时间内有大量个体被传染并且发病，同时由于该病在传染给人后会有持续近一周的潜伏期，大量前期被传染的个体四处走动，造成了总病例数急剧攀升。

大约在三个半月后，总病例增长趋势放缓，开始平稳增长，伴随小幅下降，这时该病毒的传染进入了平稳期，染病人数增长率稳定。

但是，从已知数据和对病毒的研究分析看来，至少从统计开始的200天的时间段内，感染H1N1总人数没有平缓下来的迹象，可见该病毒传染性强，应该增强防控手段。

建议：
1.提高卫生水平。

通过模型结果分析发现，生活卫生情况是决定流感蔓延是否的原因之一。

流感可以通过咳嗽或者打喷嚏传播，因此咳嗽或者打喷嚏时应该掩住口、鼻；越来越多的证据表明，流感病毒可以在一些日常用品表面上存活一段时间，有可能通过手传播到嘴、鼻、眼睛上，因此应勤洗手，还可经常用酒精为日常用品消毒。

2.完善医疗结构，确保健康生活。

根据建模中过程中发现，医疗是帮助确诊病人恢复健康的重要因素。

但由于现阶段很多地方医疗保健措施还非缺乏，导致出现不能及时就医的情况，不利于对流感的预防和治疗，于是出现不断蔓延的趋势。

3.群体免疫，增强人们对流感认识。

导致流感的快速蔓延的情况，是人们对其认识还不够深刻的具体体现。

4.打疫苗，由于传染病有一定的潜伏期，在这段期间我们很有可能接触病毒感染者，这时就只有靠我们自身的免疫。

部分MATLAB的代码：
>> dsolve('Dy=lam*y*(1-y)','y(0)=y0')
ans =
1/(1-exp(-lam*t)*(-1+y0)/y0)
dsolve('Dy=lam*y*(1-y)-miu*y','y(0)=y0')
ans =
(lam-miu)/(lam-exp(-(lam-miu)*t)*(-lam+miu+y0*lam)/y0/(lam-miu)*lam+e xp(-(lam-miu)*t)*(-lam+miu+y0*lam)/y0/(lam-miu)*miu)
y1=(lam-miu)./(lam-exp(-(lam-miu)*t)*(-lam+miu+y01*lam)./y01./(lam-mi u)*lam+exp(-(lam-miu).*t)*(-lam+miu+y01*lam)/y01/(lam-miu)*miu); >>
y2=(lam-miu)./(lam-exp(-(lam-miu)*t)*(-lam+miu+y01*lam)./y02./(lam-mi u)*lam+exp(-(lam-miu).*t)*(-lam+miu+y02*lam)/y02/(lam-miu)*miu);
hold on
>> plot(t,y1,'k-.''linewidth',3);
y=0.7;lam=0.01;miu=0.1;t=0:0.01:100;
>>
y=(lam-miu)./(lam-exp(-(lam-miu)*t)*(-lam+miu+y*lam)./y./(lam-miu)*la m+exp(-(lam-miu).*t)*(-lam+miu+y*lam)/y/(lam-miu)*miu);
>> plot(t,y,'k-','linewidth',2)
>> plot(t,y,'k-','linewidth',2)
>> ai0=0.001;
>> s0=1-ai0;
>> x0=[ai0;s0];
>> [t,x]=ode45('ill',[0,200],x0);
>> [t,x]=ode45('ill',[0,200],x0);
>> plot(t,x(:,1),'k-.',t,x(:,2),'k-');
>> text(40,0.7,'\leftarrowi(t)');
>> text(38,0.13,'\leftarrows(t)');
>> figure
>> plot(x(:,2),x(:,1))
>> plot(t,x(:,1),'k-.',t,x(:,2),'k-'); t=[4 5 6 7 8 11]
>> y=[1667
3748
42009
40486
28779
8221
]
y =
1667
3748
42009
40486
28779
8221
>> plot(t,y)。