具有时滞的传染病动力学模型Read.docx

具有时滞的传染病动力学模型Read.docx

《具有时滞的传染病动力学模型Read.docx》由会员分享，可在线阅读，更多相关《具有时滞的传染病动力学模型Read.docx（17页珍藏版）》请在冰豆网上搜索。

具有时滞的传染病动力学模型Read

具有时滞的传染病动力学模型

数值仿真

专业：

信息与计算科学

班级：

信息042

姓名：

张婧怡

学号：

04052225

指导教师：

牛大田老师

摘要：

:

具有时滞的传染病模型能较好反映传染病的潜伏期、免疫期等问题，对其研究越来越受到重视。

利用计算机模拟方法对具有时滞的传染病动力学模型进行了分析，阐述了数值仿真的基本原则。

采用时滞微分方程的数值解法对模型进行了定性分析，给出了数值仿真的应用实例。

结果表明该方法是有效的，并具有潜在的应用价值。

关键词:

传染病模型；时滞微分方程；数值仿真

一、问题背景

近年来随着环境污染、生态破坏和频繁的国际交流，国内外重人传染病爆发时有发生，如SARS,禽流感和艾滋病等，这使得对传染病的研究越来越重要。

在对传染病的诸多研究中，利用数学模型对传染病进行定性研究是个重要课题。

由于具有时滞的传染病模型能更好与实际情况相符，所以近年来对其研究受到了人们的广泛重视。

随着计算机的发展，在研究方法上，除了传统的理论分析外，计算机模拟也是研究的重要手段，一此重大发现就是通过数值仿真得到的。

虽然如此，但迄今为止利用计算机仿真方法研究时滞传染病模型的工作还很少。

我们试图在此方面对此尝试，将时滞微分方程数值算法应用于传染病模型研究，取得了较好的结果，总结出了数值仿真的基本原则，说明该方法其有潜在的应用价值。

二、模型建立

“时滞”在传染病中是个基本的因素，并在传染病的传播过程中起着重要作用，它可以反映传染病的潜伏期，患者对疾病的感染期和恢复者对玖病的免疫期等实际现象，因此使用时滞”模型更贴近实际。

如Busenberg和Cooke将时滞因素引入到由媒介传播疾病的SEIS模型中用时滞项来反映传染病的潜伏期，建立了如图1所小的仓室框图。

此模型中，把传染病地区的人群分为三类：

用S（t）,E（t）,I（t}分别表示易感者、在潜伏期的感染者和染病者。

箭头所指方向可以清楚的显小出各类人群流动的情况，T>0是模型的时滞项，代表疾病在人群中的潜伏期，r>0表小感染者被治愈后返回到易感人群中的速率，是易感者和传病媒介间的有效接触系数。

由仓库框图容易得到对应数学表达的动力学模型：

（1）

上述传染病动力学模型实质上是个具有时滞的微分方程组，对该模型进行数值仿真，就是对方程组

（1）求解，通过研究该方程组解的变化，从而得到如传染病的发展趋势等相关内容。

三、模型分析

通过对实际模型的多次数值仿真实验，得到了数值仿真的基本原则

1、运算精度优先：

由于一般传染病模型对应的方程组维数不高，对非刚性的时滞传染病模型应使用至少四阶精度的数值方法如果采用精度不高的数值方法，会造成较大的误差，从而对理论分析造成很大的偏差。

下面我们从实际的程序数值运行和图示方法加以比较。

2、模型分析程序：

一阶Euler方法

#include

#defineN10

floatf（floata,floatb）

{floatc;

c=b-2*a/b;

returnc;

}

main（）

{inti;

floath=0.2;

floatX[N]={0.0};

floatY[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

#defineN10

#defineM30

floatf（floata,floatb）

{floatc;

c=b-2*a/b;

returnc;}

main（）

{inti;

floath=0.2;

floatX[N]={0.0};

floatY[N]={1.0};

floatk[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所示的相图，

可以看出当算法精度不高时，在每步处下滑到另条解曲线、容易得出错误的分析结果所以我们采用经典四阶RK方法对模型进行数值仿真。

四、时滞微分方程的数值仿真实现

由于传染病模型多是非线性系统，因此，给出如下的时滞微分方程的一般形式:

（2）

这里

是给定向量函数。

其中

是方程的时滞项，对时滞微分方程的数值求解关键是在对时滞项的处理上下而用Runge-Kutta法给出求解上述问题的数值解法形式。

令（a,b,c）表示给定的Runge-Kutta方法其中

，

，

对应的Runge—Kutta方法具体形式如下

（3）

取步长h=

，m为正整数，将式（3}应用于方程

（2），得

（4）

如果用“~”符号表示精确解和数值解的近似，则

，

代表第n次迭代后时滞项的近似值，其中

又因为

，所以有

，则

当

时，有

（5）

以上算法的重点就是当计算

的值时，高要用到前面

的值，由于

的值也在网格点上，所以将其代替时滞项的值进行迭代运算，如图4所示。

图中双箭头表小每次迭代计算要用到前面的值。

圈4时滞徽分方程迭代过程示意图

五、数值仿真应用实例

根据出生死亡率、有无潜伏期、患病后是否有治愈或康复后是否有免疫力等因素，传染病模型还可以分为SI模型、SIS模型、SIR模型、SIRS模型等等，这里我们建立了一种带时滞SIQR模型，其简化形式如下:

（8）

其中I（t），Q（t）R（t）分别表不感染者、隔离者和恢复者，

是模型的时滞项，代表隔离者的隔离时间，

是种群规模的极限值，r>0表不隔离率，

代表种群的死亡率。

下面利用数值仿真方法来解决这个问题。

程序

#include

#include

#defineM30

#defineN120

doublef（doublet,doublei,doubleq,doubler）

{doublec;

c=-0.6*i+0.8*（1-（i+r）/（1000-q））*i+t*0;

returnc;}

doubleg（doublet,doublei,doubleq,doubler）

{doublec;

c=-0.2*q+0.4*i-0.4*i*（t-1）*exp（-0.2）+r*0;

returnc;}

doubleh（doublet,doublei,doubleq,doubler）

{doublec;

c=-0.2*r+0.4*i*（t-1）*exp（-0.2）+q*0;

returnc;}

voidmain（）

{/*doublef（doublet,doublei,doubleq,doubler）;

doubleg（doublet,doublei,doubleq,doubler）;

doubleh（doublet,doublei,doubleq,doubler）;*/

inti;

doublej=1.0;

doublek[M],l[M],m[M];

doubleT[N]={0.0000001};

doubleI[N]={0.001};

doubleQ[N]={0.001};

doubleR[N]={0.001};

for（i=1;i<=99;i++）

T[i]=T[i-1]+j;

for（i=0;i<99;i++）

{k[i]=f（T[i],I[i],Q[i],R[i]）;

m[i]=h（T[i],I[i],Q[i],R[i]）;

k[i+1]=f（T[i]+j/2.0,I[i]+j/2.0*k[i],Q[i]+j/2.0*l[i],R[i]+j/2.0*m[i]）;

m[i+1]=h（T[i]+j/2.0,I[i]+j/2.0*k[i],Q[i]+j/2.0*l[i],R[i]+j/2.0*m[i]）;

k[i+2]=f（T[i]+j/2.0,I[i]+j/2.0*k[i+1],Q[i]+j/2.0*l[i+1],R[i]+j/2.0*m[i+1]）;

m[i+2]=h（T[i]+j/2.0,I[i]+j/2.0*k[i+1],Q[i]+j/2.0*l[i+1],R[i]+j/2.0*m[i+1]）;

k[i+3]=f（T[i+1],I[i]+j*k[i+2],Q[i]+j*l[i+2],R[i]+j*m[i+2]）;

I[i+1]=I[i]+j/6.0*（k[i]+2*k[i+1]+2.0*k[i+2]+k[i+3]）;

l[i]=g（T[i],I[i],Q[i],R[i]）;

l[i+1]=g（T[i]+j/2.0,I[i]+j/2.0*k[i],Q[i]+j/2.0*l[i],R[i]+j/2.0*m[i]）;

l[i+2]=g（T[i]+j/2.0,I[i]+j/2.0*k[i+1],Q[i]+j/2.0*l[i+1],R[i]+j/2.0*m[i+1]）;

l[i+3]=g（T[i+1],I[i]+j*k[i+2],Q[i]+j*l[i+2],R[i]+j*m[i+2]）;

Q[i+1]=Q[i]+j/6.0*（l[i]+2*l[i+1]+2.0*l[i+2]+l[i+3]）;

m[i+3]=h（T[i+1],I[i]+j*k[i+2],Q[i]+j*l[i+2],R[i]+j*m[i+2]）;

R[i+1]=R[i]+j/6.0*（m[i]+2*m[i+1]+2.0*m[i+2]+m[i+3]）;

}

for（i=0;i<=99;i++）

printf（"I[%d]=%fQ[%d]=%fR[%d]=%f\n",i,I[i],i,Q[i],i,R[i]）;

}

数据

I[1]=0.122121Q[1]=0.134221R[1]=0.069789

I[2]=0.149134Q[2]=0.133561R[2]=0.082639

I[3]=0.182112Q[3]=0.089096R[3]=0.147959

I[4]=0.222359Q[4]=-0.011815R[4]=0.279218

I[5]=0.271454Q[5]=-0.186444R[5]=0.494890

I[6]=0.331302Q[6]=-0.457876R[6]=0.819675

I[7]=0.404198Q[7]=-0.856490R[7]=1.286036

I[8]=0.492883Q[8]=-1.421805R[8]=1.936112

I[9]=0.600614Q[9]=-2.204732R[9]=2.824049

I[10]=0.731220Q[10]=-3.270275R[10]=4.018781

I[11]=0.889147Q[11]=-4.700674R[11]=5.607276

I[12]=1.079468Q[12]=-6.598955R[12]=7.698154

I[13]=1.307837Q[13]=-9.092697R[13]=10.425529

I[14]=1.580343Q[14]=-12.337714R[14]=13.952698

I[15]=1.903229Q[15]=-16.521046R[15]=18.475082

I[16]=2.282420Q[16]=-21.862366R[16]=24.221431

I[17]=2.722799Q[17]=-28.612454R[17]=31.451940

I[18]=3.227217Q[18]=-37.047024R[18]=40.451451

I[19]=3.795253Q[19]=-47.453908R[19]=51.515734

I[20]=4.421852Q[20]=-60.111823R[20]=64.929024

I[21]=5.096090Q[21]=-75.259860R[21]=80.932030

I[22]=5.800441Q[22]=-93.058833R[22]=99.681703

I[23]=6.510974Q[23]=-113.548563R[23]=121.207123

I[24]=7.198820Q[24]=-136.608452R[24]=145.369235

I[25]=7.832927Q[25]=-161.930961R[25]=171.834407

I[26]=8.383717Q[26]=-189.017272R[26]=200.071364

I[27]=8.826868Q[27]=-217.200742R[27]=229.377039

I[28]=9.146315Q[28]=-245.697221R[28]=258.930121

I[29]=9.335788Q[29]=-273.674360R[29]=287.863828

I[30]=9.398646Q[30]=-300.327286R[30]=315.344654

I[31]=9.346282Q[31]=-324.947200R[31]=340.643147

I[32]=9.195679Q[32]=-346.972483R[32]=363.186038

I[33]=8.966740Q[33]=-366.017086R[33]=382.584525

I[34]=8.679919Q[34]=-381.876293R[34]=398.639018

I[35]=8.354427Q[35]=-394.513829R[35]=411.324666

I[36]=8.007094Q[36]=-404.036284R[36]=420.763955

I[37]=7.651819Q[37]=-410.660968R[37]=427.192772

I[38]=7.299473Q[38]=-414.682392R[38]=430.925305

I[39]=6.958090Q[39]=-416.441076R[39]=432.321583

I[40]=6.633239Q[40]=-416.296905R[40]=431.759903

I[41]=6.328445Q[41]=-414.607976R[41]=429.615059

I[42]=6.045624Q[42]=-411.715000R[42]=426.242397

I[43]=5.785465Q[43]=-407.930742R[43]=421.967115

I[44]=5.547762Q[44]=-403.533696R[44]=417.077961

I[45]=5.331688Q[45]=-398.765115R[45]=411.824387

I[46]=5.136002Q[46]=-393.828555R[46]=406.416277

I[47]=4.959213Q[47]=-388.891174R[47]=401.025468

I[48]=4.799698Q[48]=-384.086192R[48]=395.788445

I[49]=4.655793Q[49]=-379.516036R[49]=390.809709

I[50]=4.525847Q[50]=-375.255816R[50]=386.165465

I[51]=4.408273Q[51]=-371.356885R[51]=381.907375

I[52]=4.301566Q[52]=-367.850312R[52]=378.066203

I[53]=4.204330Q[53]=-364.750175R[53]=374.655248

I[54]=4.115282Q[54]=-362.056599R[54]=371.673515

I[55]=4.033257Q[55]=-359.758525R[55]=369.108582

I[56]=3.957212Q[56]=-357.836186R[56]=366.939171

I[57]=3.886221Q[57]=-356.263312R[57]=365.137413

I[58]=3.819474Q[58]=-355.009052R[58]=363.670841

I[59]=3.756270Q[59]=-354.039660R[59]=362.504106

I[60]=3.696010Q[60]=-353.319917R[60]=361.600447

I[61]=3.638190Q[61]=-352.814347R[61]=360.922923

I[62]=3.582393Q[62]=-352.488216R[62]=360.435435

I[63]=3.528281Q[63]=-352.308330R[63]=360.103538

I[64]=3.475586Q[64]=-352.243665R[64]=359.895071

I[65]=3.424101Q[65]=-352.265825R[65]=359.780624

I[66]=3.373671Q[66]=-352.349350R[66]=359.733840

I[67]=3.324185Q[67]=-352.471894R[67]=359.731594

I[68]=3.275571Q[68]=-352.614293R[68]=359.754046

I[69]=3.227785Q[69]=-352.760520R[69]=359.784601

I[70]=3.180804Q[70]=-352.897576R[70]=359.809780

I[71]=3.134626Q[71]=-353.015305R[71]=359.819029

I[72]=3.089258Q[72]=-353.106161R[72]=359.804479

I[73]=3.044716Q[73]=-353.164949R[73]=359.760675

I[74]=3.001022Q[74]=-353.188538R[74]=359.684278

I[75]=2.958196Q[75]=-353.175567R[75]=359.573767

I[76]=2.916261Q[76]=-353.126157R[76]=359.429139

I[77]=2.875233Q[77]=-353.041628R[77]=359.251618

I[78]=2.835129Q[78]=-352.924236R[78]=359.043388

I[79]=2.795958Q[79]=-352.776928R[79]=358.807344

I[80]=2.757726Q[80]=-352.603125R[80]=358.546864

I[81]=2.720433Q[81]=-352.406531R[81]=358.265621

I[82]=2.684076Q[82]=-352.190968R[82]=357.967409

I[83]=2.648645Q[83]=-351.960239R[83]=357.656006

I[84]=2.614128Q[84]=-351.718019R[84]=357.335062

I[85]=2.580506Q[85]=-351.467765R[85]=357.008008

I[86]=2.547761Q[86]=-351.212655R[86]=356.677996

I[87]=2.515870Q[87]=-350.955545R[87]=356.347851

I[88]=2.484808Q[88]=-350.698941R[88]=356.020044

I[89]=2.454549Q[89]=-350.444987R[89]=355.696686

I[90]=2.425065Q[90]=-350.195471R[90]=355.379521

I[91]=2.396328Q[91]=-349.951832R[91]=355.069948

I[92]=2.368311Q[92]=-349.715181R[92]=354.769032

I[93]=2.340986Q[93]=-349.486328R[93]=354.477537

I[94]=2.314324Q[94]=-349.265809R[94]=354.195950

I[95]=2.288299Q[95]=-349.053921R[95]=353.924517

I[96]=2.262885Q[96]=-348.850751R[96]=353.663276

I[97]=2.238058Q[97]=-348.656212R[97]=353.412090

I[98]=2.213792Q[98]=-348.470075R[98]=353.170678

I[99]=2.190066Q[99]=-348.291996R[99]=352.938649

将数值形式（4）应用到模型（8），在Matlab6.5环境下，采用四阶Runge-Kutta方法仿真求解。

分别取

时，可以验证参数值满足条件（9）,其仿真结果如图5所不，图6显示当I,Q,R取不同的初始值时，仍然有I,Q,R逐渐趋近

。

因此，可以看出随着时间增大时，（I,Q,R）逐渐趋近

。

用数值方法还可验证，在满足条件（9）的情况下任意改变参数，系统仍然保证有（I,Q,R）仍趋向

。

因此，条件（9}是否为地方病平衡点的全局渐进稳定性条件。

六、误差分析

由以上分析可以看出，模型存在误差

七、模型改进及结论

由于具有时滞的传染病模型的复杂性，还有许多的问题尚待解决。

如种群规模变动的具有时滞的传染病模型的动力