【建模教程】-传染病模型建模

----宋停云与您分享--------宋停云与您分享----对传染病的传播的研究摘 要本文以常见传染病的传播为研究方向，并结合微分方程的知识建立传染病的传播与控制模型。在模型的基础上，运用MATLAB软件拟合出患者人数与时间的关系曲线，从而能够从图中直观地对该病的传播作出分析并提出应对措施。在问题二、三中，利用所建立的微分方程模型代入给出的数据，从而用MATLAB拟合出各项人数随时间的变化曲线，分析所得图形及其合理性，得到有关该传染病的信息。Th传染病 微分方程模型 MATLAB 曲线拟合 应对措施----宋停云与您分享--------宋停云与您分享----目录对传染病的传播的研究 2摘 要 2一、问题重述 4相关情况 4问题的提出 4二、模型假设 4三、符号的约定和说明 5四、对问题一的解答 5问题分析 5模型准备 6模型的建立 7五、对问题二的解答 8问题分析 8模型的建立 8结果分析 9六、对问题三的解答 9问题分析 9模型的建立 9结果分析 10对隔离强度p的灵敏度分析 11七、对问题四的解答 11八、模型的评价及推广 12模型的优缺点 12模型的推广 13九、附录： 14----宋停云与您分享--------宋停云与您分享----一、问题重述.相关情况2013----宋停云与您分享--------宋停云与您分享----毒潜伏期为a

a天，患病者的治愈时间为a

天，假设该病毒可以通过人与人----宋停云与您分享--------宋停云与您分享----1 2 3之间的直接接触，患者每天接触的人数为r，因接触被感染的概率为 (为感染率)．为了控制疾病的传播与扩散，将人群分成五类，患者、疑似患者、治愈者、死亡者、正常人．潜伏期内的患者被隔离的比例为p（为潜伏期内患者被隔离的百分数）。.问题的提出问题一：在合理的假设下建立该病毒扩散与传播的控制模型。问题一：在合理的假设下建立该病毒扩散与传播的控制模型。问题二：利用你所建立的模型对如下数据进行模拟：a3,a 7,a 60,r10,p30%,50%210，12322变化的曲线，并分析所给结果的合理性。3080，患者人数将有何变化问题四：请据此模型，给出控制此传染病传播的建议。二、模型假设1、在该传染病的考察期内，该考察地区的总人数为常数，不考虑人口流动。2、将病毒的所有传播途径都视为与病原体的直接接触造成。3、忽略该考察时间内人口的自然死亡率和出Th率。4、被隔离的人群完全断绝与外界的接触，因此不具有传染性。5、被治愈者获得抗体，不考虑其二次传染患病。----宋停云与您分享--------宋停云与您分享----6、将治愈者和死亡者定义为非参与者，即退出研究的传染病传播体系。7、疑似患者即为潜伏期的患者，是被有效接触后具有传染性且传染概率也为，经过隔离治疗可转为治愈者（非参与者，治愈时间为a天。38、潜伏期患者和确诊患者接触传染的均为易感病正常人，且均将其传染为潜伏期患者。三、符号的约定和说明I：确诊患者E：潜伏期患者（即疑似患者)R：非参与者（痊愈和死亡的患者)S：普通易感的正常人：潜伏期患者和确诊患者的传染概率----宋停云与您分享--------宋停云与您分享----a~a1 2

：传染性病毒的潜伏期----宋停云与您分享--------宋停云与您分享----a：潜伏期患者和确诊患者被治愈的时间3N：该地区总人数r：该人群的人均每天接触人数p：潜伏期内患者被隔离的百分数四、对问题一的解答问题分析ThI：确诊患者E：潜伏期患者（被病毒有效接触后有传染性的人)R：非参与者（痊愈和死亡的患者）S：普通易感者（正常人）建立上述五种情况的人数在单位时间变化的微分方程模型。----宋停云与您分享--------宋停云与您分享----由上述五类得到以下关系图：图表1（非参与者被未隔离潜伏期患者传染正常人SE潜伏期确诊患者I被未隔离确诊患者传染隔离治疗（非参与者模型准备（1）单位时间内正常人数变化：易感正常人与未隔离潜伏期病人及确诊患者接触后均变为潜伏期患者，结合以上所给信息故 NdSN IpSNE(1p)Sdt----宋停云与您分享--------宋停云与您分享----即 dSdt

I1pSE(p)S ········（1）----宋停云与您分享--------宋停云与您分享----（2）单位时间内潜伏期患者（疑似患者）人数变化：潜伏期患者的数量变化为正常人被感染为潜伏期患者人数减去潜伏期患者被治愈和转为确诊患者的人数，结合以上信息----宋停云与您分享--------宋停云与您分享----II

pS

2 EE(1p)S

1Ep···（2）----宋停云与您分享--------宋停云与您分享----dt aa a1 2 3（3）单位时间内确诊患者人数变化：确诊患者人数为潜伏期患者转变人数减去被治愈人数----宋停云与您分享--------宋停云与您分享--------宋停云与您分享--------宋停云与您分享----dIdt2aa1a即 dIdt2aa1a

··············（3）----宋停云与您分享--------宋停云与您分享----1 2 3（4）单位时间内非参与者的人数变化：非参与者人数为确诊患者被治愈人数或死亡数----宋停云与您分享--------宋停云与您分享----1即 dR I1

················（4）----宋停云与您分享--------宋停云与您分享----dt a3（5）总人数： SREIN·············（5）----宋停云与您分享--------宋停云与您分享----对模型的部分说明：121222aa2aa

aa

，即单位时间内潜伏期病人以比例常数----宋停云与您分享--------宋停云与您分享----11 21----宋停云与您分享--------宋停云与您分享----2、确诊病人平均死亡或痊愈的疗程为a3

,即单位时间感染者的治愈率为a131

0；----宋停云与您分享--------宋停云与您分享--------宋停云与您分享--------宋停云与您分享----3、潜伏期患者平均疗程为a3

，即单位时间内潜伏期患者的治愈率为a3

0；----宋停云与您分享--------宋停云与您分享----4、单位时间内每个易感者与病人的接触率参数为r0；模型的建立----宋停云与您分享--------宋停云与您分享----dSdt

I pSE(1 p)S----宋停云与您分享--------宋停云与您分享----dE I pS 2dt a a 1 2

EE(1 p)S 1 Epa3----宋停云与您分享--------宋停云与您分享----dI

2 E 1 I----宋停云与您分享--------宋停云与您分享----dt a1a2 a3dR 1 Idt a 3SEIR N----宋停云与您分享--------宋停云与您分享----I(0)I0

，S(0)S0

，R(0)R0

，A(0)A0

，E(0)E0

N1107为----宋停云与您分享--------宋停云与您分享----系统中各类的初始值。五、对问题二的解答问题分析MATLAB模型的建立----宋停云与您分享--------宋停云与您分享----a3,a1 2

7,a3

60,r10,p30%,50%,N1107,初始发病人数 100----宋停云与您分享--------宋停云与您分享----（即I 100），疑似患者210（即E 210），患者2天后入院，疑似患者20 0天后被隔离。这样可以得到患者人数随时间变化的曲线（如下图）：图表2----宋停云与您分享--------宋停云与您分享----结果分析第13.31天）患者人数达到最大值78170001002006000六、对问题三的解答问题分析该问题是建立在问题二的基础上，利用问题二建立的模型，提高p值得到新的患模型的建立----宋停云与您分享--------宋停云与您分享----问题三中a1

3,a2

7,a3

60,r10,p80%,50%,N1107----宋停云与您分享--------宋停云与您分享----100（I0100），210（E0210），22天后被隔离。这样可得到与问题二的对比图（绿线为p80%的图，蓝线为p30%的图）：----宋停云与您分享--------宋停云与您分享----结果分析分析绿线可以看到病毒传染初期患者人数依然急剧升高，最高点（第12.65天）患者人数达到最大值7555000人，100天后下降到1921000人。(天）大值（人）度p患病入院天数n（天）人）问题二13.317.817×10^630210问题三12.657.555×10^680210p80p即绿线与蓝线p80p30相比，达到最高点100p30时人数。从对比中可以看出，提高隔离强度可以更好地控制疫情，减少患病人数。----宋停云与您分享--------宋停云与您分享----p的灵敏度分析对于该问题中的隔离强度p由30变为80，通过问题三的图像中两条曲线pp此政府和意愿应尽量增大隔离强度p。七、对问题四的解答根据题目我们建立出模型：----宋停云与您分享--------宋停云与您分享----dSdt

I pSE(1 p)S----宋停云与您分享--------宋停云与您分享----dE I pS 2dt a a 1 2

EE(1 p)S 1 Epa3----宋停云与您分享--------宋停云与您分享----1dI1dt

2a a1 2

E Ia3----宋停云与您分享--------宋停云与您分享----dR 1 Idt a 3SEIR N从模型中我们可以看出正常人的减少是由于被潜伏期患者以及确诊患者的----宋停云与您分享--------宋停云与您分享----的接触。模型中p越大，S越少，而E增加也越来越少，因此对疫情的控制有很好的效果。Th84大隔离防治措施力度，且改进医疗手段，使治愈时间a3减小来控制疫情。八、模型的评价及推广模型的优缺点模型的优点：(1)、将医学领域的问题转化到数学领域进行分析和讨论，可以清楚地定量地得出传染病的发展趋势和高峰以及未来的预测，具有很强的可靠性和实用性。(2)、模型中各个变量的关系明确，易于模型的求解。(3)、本文的数学模型是以连续的微分方程为基础，不会得出准确的解析解，本模型的缺点：(1)、采用微分方程方法建立数学模型，易受外界因素变化的影响，其稳定性具有相对性。----宋停云与您分享--------宋停云与您分享----(2)、模型中的参数变量有其自身的随机性，虽然本文对已知数据进行统计平均的处理方法，但在计算过程中存在误差。(3)、模型中涉及的参数较多，在实际Th活中很难确定各参数，因此模型具有理想化。模型的推广SIRThTh：[1].杨启凡，《数学建模》，浙江：浙江大学出版社，2006.6[2]MATLAB[3].姜启源，数学建模案例选集，北京：高等教育出版社，2006.7[4]（第三版[5]./p-692135025731.html[6].梁国业等，《数学建模》，北京：冶金工业出版社，2004.9[7].韩中庚等，《数学建模方法及其应用》，北京：高等教育出版社，2005.6[8]MATLAB2009.4----宋停云与您分享--------宋停云与您分享----九、附录：问题二新建m文件夹:functionx=pencil(t,x)%s=x(1)e=x(2)i=x(3)r=x(4);a1=3;a2=7;a3=60;p=0.3;m=10;b=0.5;x=[-b*x(3)*(1-p)*x(1)-b*x(2)*(1-p)*x(1),b*x(3)*(1-p)*x(1)-2/(a1+a2)*x';在命令窗口内输入：s0=[10000000,100,0,210];[t,x]=ode23s(@pencil,[0,100],s0)plot(t,x(:,3));holdontext(0,100,'(0,100)','color','r')text(13.31,7.817e+006,'(13.31,7.817e+006)','color','r')text(100,2.006e+006,'(100,2.006e+006)','color','r')plot(0,100,'g+',13.31,7.817e+006,'g+',100,2.006e+006,'g+')问题三新建m文件夹:functionx=pencil(t,x)%s=x(1)e=x(2)i=x(3)r=x(4);a1=3;a2=7;a3=60;p=0.3;m=10;b=0.5;x=[-b*x(3)*(1-p)*x(1)-b*x(2)*(1-p)*x(1),b*x(3)*(1-p)*x(1)-2/(a1+a2)*x';----宋停云与您分享--------宋停云与您分享----functionx=pen(t,x)%s=x(1)e=x(2)i=x(3)r=x(4);a1=3;a2=7;a3=60;p=0.8;m=10;b=0.5;x=[-b*x(3)*(1-p)*x(1)-b*x(2)*(1-p)*x(1),b*x(3)*(1-p)*x(1)-2/(a1+a2)*x';命令窗口输入：s0=[10000000,100,0,210];[t,x]=ode23s(@pencil,[0,100],s0)plot(t,x(:,3));holdontext(0,100,'(0,100)','color','r')text(13.31,7.817e+006,'(13.31,7.817e+006)','color','r')text(100,2.006e+006,'(100,2.006e+006)','color','r')plot(0,100,'g+',13.31,7.817e+006,'g+',100,2.006e+006,'g+')s0=[10000000,100,0,210];[t,x]=ode23s(@pen,[0,100],s0)plot(t,x(:,3),'g');holdontext(0,100,'(0,100)','color','r')text(12.65,7.555e+006,'(12.65,7.555e+006)','color','r')text(100,1.921e+006,'(100,1.921e+006)','color','r')plot(0,100,'g+',12.65,7.555e+006,'g+',100,1.921e+006,'g+')----宋停云与您分享--------宋停云与您分享----Per3D技术文件Modelingphase中的工作流程和细节sapetabspkpmcad中导入几何信息。（1）注意连梁一定要用梁单元模拟，不要用壳单元。从pkpm中导入的模型如对于深连梁的处理办法：（2）所有几何尺寸不同或配筋不同的构件都有分组（前期较粗可以只按照几何尺寸分组，以便后续分组导入。（3）用etabs分组导处perform文本时一定要计算一下恒载，然后导出文本，并注意修改文本前面不要有空行。（4）perform下建模菜单中的“Elements”new，定义需要导入组中构（可与etabs中的分组名一致（5）perform下建模菜单中的“Import/Exportstructuredata”选单元类型，按（vistC盘或其他盘的根目录下以便寻找skip行填实际文本的跳过的数据行，形成几何模型。perform下建模菜单中的“Nodes”完成以下工作：（1）在supports中完成边界约束。（2）massetabsD+0.5L每层数值按节点数分配到（N-m-k。（3）在slaving中定义刚性隔壁：注意要每层一个刚性隔壁名。perform下建模菜单中的“Loadpatterns”完成以下工作：1elementsloadetabs中的梁荷载。定了Dead和Live两种工况。（2）自重按照单元组来定义其具有自重（componentproperties中可以定义各----宋停云与您分享--------宋停云与您分享----组单元的密度。注意延梁不需要定义自重。（3）其他荷载工况是也在此处定义，如静力风、雪等，在后续分析结果中也可以进行组合。perform下建模菜单中的“Addordeleteframe”完成以下工作：perform下建模菜单中的“componentproperties”完成以下工作：（1）定义柱单元：CISSectionBuildingP-MMoment-Curvature关系曲线（kN-m的单位）和基本几何特性。CrossSects中选“ColumnReinforcedconcreteSection”new“Sect_组名”并ok，右侧出现三个菜单，需填写前两个：在“StiffnessDimension”中填写柱的几何尺寸和几何特性（型钢柱需从Sn”中抄写过来，一般柱可让其自己计算；在“InelasticStrength”中只填写“Axial/Bending”CISSectionBuildingP-MPTPCPB2轴（0度）3轴（90度）PB0轴力对应的弯矩。plotP-M关系曲线。其他选项可不定义。Inelastic中选“P-M2-M3HingeConcreteCurvaturenew菜单，输入“Hinge_组名”okUseCrossSectionDeformationCapabilities点选yes，右侧需填写三个菜单：在“SectionandDimension”中选“Column,ReinforcedconcreteSection”及定义好的“Sect_组名”在“BasicF-DRelationship”PB2（0度3（90度）Moment-CurvaturePBDX；在“Deformation Capabilities”填写三水准下的曲率，读取Moment-Curvature关系曲线中对应的三个曲率值。----宋停云与您分享--------宋停云与您分享----Compound中选“FrameMemberCompoundComponentnew菜单，输入“Com_组名”ok，右侧需填写两个菜单：在“BasicCompounds”CompoundTypeP-M2-M3HingeConcreteCurvatureCompoundName中选已经定义好的“Hinge_组名”，在“Tributarylength，proportion”0.02ADD添加；然后再在CompoundTypeColumn,ReinforcedconcreteSectionCompoundName中选已经定义好的“Sect_组名”，在“Tributarylength，proportion”ADD（该过程重复一次CompoundTypeP-M2-M3HingeConcreteCurvatureTypeCompoundName中选已经定义好的_组名yproportion0.98，ADD添加。完成柱的塑性铰的定义。SelfWeight中定义单位长度的密度，注意长度单位。（2）定义梁单元：CISSectionBuildingMoment-Curvature（一kN-m的单位）和基本几何特性。CrossSects中选“BeamReinforcedconcreteSection”new菜单，输入“Sect_组名”并ok，右侧出现三个菜单，需填写第一个：在“StiffnessDimension”中填写梁的几何尺寸和几何特性（型钢梁需从Sn中抄写过来，一般柱可让其自己计算；同时要主要连梁、边梁和中梁的放大系数。Inelastic中选“MomentHinge,Curvaturenew“Hinge_组名”okUseCrossSectionDeformationCapabilitiesyes，右侧需填写三个菜单：在“SectionandDimension”中选“Beam,ReinforcedconcreteSection”及定义好的“Sect_组名”在“BasicF-DRelationship”FUDXMoment-Curvature关系曲线中均可以读出；在“Deformation Capabilities”填写三水准下的曲率，读取----宋停云与您分享--------宋停云与您分享----Moment-Curvature关系曲线中对应的三个曲率值（第一水准一般取弹性1.521.5235）Compound中选“FrameMemberCompoundComponentnew菜单，输入“Com_组名”ok，右侧需填写两个菜单：在“BasicCompounds”CompoundTypeMomentHinge,CurvatureCompoundName中选已经定义好的“Hinge_组名“Tributarylengtproportion”0.02ADDCompoundTypeBeam,ReinforcedconcreteSectionCompoundName中选已经定义好的“Sect_组名”，在“Tributarylength，proportion”中输入0.5，点击ADD（该过程重复一次CompoundTypeMomentHinge,CurvatureCompoundName中选已经定义好的“Hinge_组名”，在“Tributarylength，proportion”0.98ADD添加。完成梁的塑性铰的定义。SelfWeightYL不需填该项。3）定义墙单元（关键是定义混凝土、钢筋和剪切材料的本构：Material中选钢筋材料“InelasticConcreteMaterial”new菜单，输入代“C60”ok，在StrainCapacitiesStrengthLoss中ShapeofRelationshipTrilinear。右侧需填写三个菜单：BasicRelationship中定义混凝土的弹性模量、本构关系图中对于的压应力和应变值，具体可参算例。可查规范确定各数值。StrengthLoss中定义本构关系中的压应变，具体可参算例。为保障计算的收敛可以不考虑混凝土的强度损失。StrainCapacities中定义三水准混凝土的应变，具体可参算例。Material中选混凝土材料“InelasticSteelMaterial,Non-Buckling”new并ok，在StrainCapacities和StrengthLossyesBasicRelationship中定义钢筋的弹性模量、极限应力和应变（与本构关系图对应，可查规范。----宋停云与您分享--------宋停云与您分享----StrengthLoss中定义本构关系中的拉应变，具体可参算例。可以暂不考虑。StrainCapacities中定义三水准钢筋的应变，具体可参算例。Material中选弹性剪切材料“ElasticShearMaterialforawall”new菜单，输入代表剪切材料的名称“ShearMaterial”ok，在右侧菜单只需要填（非线性剪切材料定义后费机时，不易收敛）CrossSects中选“ShearWallInelasticSectionnew菜单，输入“Sect_组名”AutoSizeok，在右侧菜单中直接选已经定义好的混凝土和钢筋材料作为剪力墙纤维束，同时填写墙厚度、纤维数量（4~8为宜，模33个纤维数，以简便考虑边缘约束构件的作用。Compound中选“ShearWallCompoundComponent”new菜单，输入“Com_组名”并ok，右侧填写基本信息即可。4）定义连梁单元（关键是定义一个剪切铰，但没有定型的剪切骨架曲线：Inelastic中选“ShearHinge,Displace