【建模教程】-甲型 H1N1 流感傳播模型研究

----宋停云與您分享--------宋停云與您分享----H1N1摘要SIRh1n1流感病毒的傳播規(guī)律進(jìn)行了一、問題重述h1n1（又稱豬流感）正成為人們關(guān)注的h1n1二、問題分析h1n1SIRSIRSIRh1n1H1N1時間確診（包括死亡病例）死亡（累計(jì)）42350424804251104262004274004286404299104301091511411521601532261542791554031----宋停云與您分享--------宋停云與您分享----56642257896258163925922542510253235112600351230093513335245144298451547144三、建立模型（一、不考慮潛伏期的數(shù)學(xué)模型1、模型假設(shè)1H1N1NIR(括死ts(t)、i(t)、r(t)。2、i(t)關(guān)于時間的增長率與s(t)成正比，比例常數(shù)為 。病人的數(shù)量減少速度與當(dāng)時的病人總?cè)藬?shù)成正比，比例常數(shù)為。治愈的病人具有了免疫力，即治愈后不再會成為二次患者。3、s(t)、r(t)、i(t)1。2、模型構(gòu)成的易感者數(shù)為S(t，NI(t個發(fā)病者平均每天能使S(t)NI(t個易感者成為病毒潛伏者。所以有：dS(t)S(t)I(t) （1）dt單位時間內(nèi)退出者的變化等于發(fā)病人群的減少，即dR(t)I(t) （2）dt發(fā)病人群的變化等于易感人群轉(zhuǎn)入的數(shù)量，即dI(t)S(t)I(t)I(t) （3）dt----宋停云與您分享--------宋停云與您分享----SR0 3、模型求解

（R0

=0。----宋停云與您分享--------宋停云與您分享----方程組（1（2（3）無法求出解析解，我們定義一個新的變量 /，于是可以求出方程的解為：----宋停云與您分享--------宋停云與您分享----i(s0

i)s10

slns （4）0----宋停云與您分享--------宋停云與您分享----下面分析s(t)、i(t)、r(t)的變化情況：----宋停云與您分享--------宋停云與您分享----aSR0 0

如何，病人最終將消失，即i

0。----宋停云與您分享--------宋停云與您分享--------宋停云與您分享--------宋停云與您分享----b、最終未被感染者的健康者的比例是s

，是方程----宋停云與您分享--------宋停云與您分享----(sis1lns0在(0,1/內(nèi)的根。0 0 s0----宋停云與您分享--------宋停云與您分享----Cs0

1/i(ts0

1/時，i(t)達(dá)到最大值，然后i(t)----宋停云與您分享--------宋停云與您分享----減小且趨于零，s(t)則單調(diào)減小至s 。d、若s 1/，則i(t)單調(diào)減小至5，s(t)則單調(diào)減小至s 。0 我們發(fā)現(xiàn)人們的衛(wèi)生水平越高，日接觸率越小；醫(yī)療水平越高，日治愈率越高，于是越小，所以提高衛(wèi)生水平和醫(yī)療水平有利于傳染病的蔓延。結(jié)合的具體情況和假設(shè)條件進(jìn)行分析：根據(jù)所得的數(shù)據(jù)畫出患病人數(shù)變化曲線和治愈人數(shù)變化曲線：----宋停云與您分享--------宋停云與您分享--------宋停云與您分享--------宋停云與您分享----h1n1s0

1/----宋停云與您分享--------宋停云與您分享----0的情況，即/1/s。由假設(shè)條件可知的取值范圍在1.4~1.6之間?，F(xiàn)0----宋停云與您分享--------宋停云與您分享----在我們?nèi)?1.6，則表示/(1/s0

)1.6，即每天平均治愈的人數(shù)最多為----宋停云與您分享--------宋停云與您分享----1.61.642336.8人，這與實(shí)際的情況相差甚遠(yuǎn)。產(chǎn)生這個問題的原因有以下幾個方面：第二：疾病預(yù)防與控制中心所得到的數(shù)據(jù)具有滯后性。----宋停云與您分享--------宋停云與您分享----第三：在s0

0不一定成立?？梢园涯切┥眢w強(qiáng)壯的、注意自己個人衛(wèi)生的----宋停云與您分享--------宋停云與您分享----人排除在外。（二、考慮潛伏期的數(shù)學(xué)模型1、模型假設(shè)1H1N1NE,IR(括死亡者和治愈者)四類，時刻t內(nèi)這三類人在總?cè)藬?shù)中所占比例分別為setitr。2、每個病人每天有效接觸的平均人數(shù)為，稱為日接觸率，當(dāng)已感染者與易2、模型構(gòu)成觸的易感者數(shù)為(t)S(t，NI(t個發(fā)病者平均每天能使(t)S(t)NI(t個易感者成為病毒潛伏者。----宋停云與您分享--------宋停云與您分享----N

dS(t)dt

(t)S(t)NI(t)

化簡得：----宋停云與您分享--------宋停云與您分享----dS(t)(t)S(t)I(t)dt病毒潛伏人群的變化等于易感人群轉(zhuǎn)入數(shù)量減去轉(zhuǎn)化為發(fā)病人群的數(shù)量，即dE(t)(t)S(t)I(t)(t)E(t)dt其中(t)表示潛伏期日發(fā)病率，即每個潛伏者平均有效發(fā)病的人數(shù)。----宋停云與您分享--------宋停云與您分享----單位時間內(nèi)退出者的變化等于發(fā)病人群的減少，即dR(t)(t)I(t)dt其中(t)表示日退出率，即每個病人平均有效病情結(jié)束的人數(shù)。發(fā)病人群的變化等于潛伏人群轉(zhuǎn)入的數(shù)量，即dI(t)(t)E(t)(t)I(t)dts(t)i(t)r(t)1初始時刻易感染者，已感染者與病愈免疫者的比例分別是----宋停云與您分享--------宋停云與您分享----s(s0 0

0),i(i0 0

0),r00----宋停云與您分享--------宋停云與您分享----3、模型求解由于潛伏期的人群數(shù)量不能確定，所以可視為是易感人群的一部分，因此求解過程跟忽略潛伏期的一樣。四、模型的改進(jìn)就如何確定日接觸率：----宋停云與您分享--------宋停云與您分享----[1]姜啟源謝金星葉俊數(shù)學(xué)建模（第四版）高等教育出版社[2]數(shù)據(jù)來源：疾病預(yù)防控制中心----宋停云與您分享--------宋停云與您分享----Per3D技術(shù)文件Modelingphase中的工作流程和細(xì)節(jié)sapetabspkpmcad中導(dǎo)入幾何信息。（1）注意連梁一定要用梁單元模擬，不要用殼單元。從pkpm中導(dǎo)入的模型如對于深連梁的處理辦法：（2）所有幾何尺寸不同或配筋不同的構(gòu)件都有分組（前期較粗可以只按照幾何尺寸分組，以便后續(xù)分組導(dǎo)入。（3）用etabs分組導(dǎo)處perform文本時一定要計(jì)算一下恒載，然后導(dǎo)出文本，并注意修改文本前面不要有空行。（4）perform下建模菜單中的“Elements”new，定義需要導(dǎo)入組中構(gòu)（可與etabs中的分組名一致（5）perform下建模菜單中的“Import/Exportstructuredata”選單元類型，按（vistC盤或其他盤的根目錄下以便尋找skip行填實(shí)際文本的跳過的數(shù)據(jù)行，形成幾何模型。perform下建模菜單中的“Nodes”完成以下工作：（1）在supports中完成邊界約束。（2）massetabsD+0.5L每層數(shù)值按節(jié)點(diǎn)數(shù)分配到（N-m-k。（3）在slaving中定義剛性隔壁：注意要每層一個剛性隔壁名。perform下建模菜單中的“Loadpatterns”完成以下工作：1elementsloadetabs中的梁荷載。定了Dead和Live兩種工況。（2）自重按照單元組來定義其具有自重（componentproperties中可以定義各----宋停云與您分享--------宋停云與您分享----組單元的密度。注意延梁不需要定義自重。（3）其他荷載工況是也在此處定義，如靜力風(fēng)、雪等，在后續(xù)分析結(jié)果中也可以進(jìn)行組合。perform下建模菜單中的“Addordeleteframe”完成以下工作：perform下建模菜單中的“componentproperties”完成以下工作：（1）定義柱單元：CISSectionBuildingP-MMoment-Curvature關(guān)系曲線（kN-m的單位）和基本幾何特性。CrossSects中選“ColumnReinforcedconcreteSection”new“Sect_組名”并ok，右側(cè)出現(xiàn)三個菜單，需填寫前兩個：在“StiffnessDimension”中填寫柱的幾何尺寸和幾何特性（型鋼柱需從Sn”中抄寫過來，一般柱可讓其自己計(jì)算；在“InelasticStrength”中只填寫“Axial/Bending”CISSectionBuildingP-MPTPCPB2軸（0度）3軸（90度）PB0軸力對應(yīng)的彎矩。plotP-M關(guān)系曲線。其他選項(xiàng)可不定義。Inelastic中選“P-M2-M3HingeConcreteCurvaturenew菜單，輸入“Hinge_組名”okUseCrossSectionDeformationCapabilities點(diǎn)選yes，右側(cè)需填寫三個菜單：在“SectionandDimension”中選“Column,ReinforcedconcreteSection”及定義好的“Sect_組名”在“BasicF-DRelationship”PB2（0度3（90度）Moment-CurvaturePBDX；在“Deformation Capabilities”填寫三水準(zhǔn)下的曲率，讀取Moment-Curvature關(guān)系曲線中對應(yīng)的三個曲率值。----宋停云與您分享--------宋停云與您分享----Compound中選“FrameMemberCompoundComponentnew菜單，輸入“Com_組名”ok，右側(cè)需填寫兩個菜單：在“BasicCompounds”CompoundTypeP-M2-M3HingeConcreteCurvatureCompoundName中選已經(jīng)定義好的“Hinge_組名”，在“Tributarylength，proportion”0.02ADD添加；然后再在CompoundTypeColumn,ReinforcedconcreteSectionCompoundName中選已經(jīng)定義好的“Sect_組名”，在“Tributarylength，proportion”ADD（該過程重復(fù)一次CompoundTypeP-M2-M3HingeConcreteCurvatureTypeCompoundName中選已經(jīng)定義好的_組名yproportion0.98，ADD添加。完成柱的塑性鉸的定義。SelfWeight中定義單位長度的密度，注意長度單位。（2）定義梁單元：CISSectionBuildingMoment-Curvature（一kN-m的單位）和基本幾何特性。CrossSects中選“BeamReinforcedconcreteSection”new菜單，輸入“Sect_組名”并ok，右側(cè)出現(xiàn)三個菜單，需填寫第一個：在“StiffnessDimension”中填寫梁的幾何尺寸和幾何特性（型鋼梁需從Sn中抄寫過來，一般柱可讓其自己計(jì)算；同時要主要連梁、邊梁和中梁的放大系數(shù)。Inelastic中選“MomentHinge,Curvaturenew“Hinge_組名”okUseCrossSectionDeformationCapabilitiesyes，右側(cè)需填寫三個菜單：在“SectionandDimension”中選“Beam,ReinforcedconcreteSection”及定義好的“Sect_組名”在“BasicF-DRelationship”FUDXMoment-Curvature關(guān)系曲線中均可以讀出；在“Deformation Capabilities”填寫三水準(zhǔn)下的曲率，讀取----宋停云與您分享--------宋停云與您分享----Moment-Curvature關(guān)系曲線中對應(yīng)的三個曲率值（第一水準(zhǔn)一般取彈性1.521.5235）Compound中選“FrameMemberCompoundComponentnew菜單，輸入“Com_組名”ok，右側(cè)需填寫兩個菜單：在“BasicCompounds”CompoundTypeMomentHinge,CurvatureCompoundName中選已經(jīng)定義好的“Hinge_組名“Tributarylengtproportion”0.02ADDCompoundTypeBeam,ReinforcedconcreteSectionCompoundName中選已經(jīng)定義好的“Sect_組名”，在“Tributarylength，proportion”中輸入0.5，點(diǎn)擊ADD（該過程重復(fù)一次CompoundTypeMomentHinge,CurvatureCompoundName中選已經(jīng)定義好的“Hinge_組名”，在“Tributarylength，proportion”0.98ADD添加。完成梁的塑性鉸的定義。SelfWeightYL不需填該項(xiàng)。3）定義墻單元（關(guān)鍵是定義混凝土、鋼筋和剪切材料的本構(gòu)：Material中選鋼筋材料“InelasticConcreteMaterial”new菜單，輸入代“C60”ok，在StrainCapacitiesStrengthLoss中ShapeofRelationshipTrilinear。右側(cè)需填寫三個菜單：BasicRelationship中定義混凝土的彈性模量、本構(gòu)關(guān)系圖中對于的壓應(yīng)力和應(yīng)變值，具體可參算例?？刹橐?guī)范確定各數(shù)值。StrengthLoss中定義本構(gòu)關(guān)系中的壓應(yīng)變，具體可參算例。為保障計(jì)算的收斂可以不考慮混凝土的強(qiáng)度損失。StrainCapacities中定義三水準(zhǔn)混凝土的應(yīng)變，具體可參算例。Material中選混凝土材料“InelasticSteelMaterial,Non-Buckling”new并ok，在StrainCapacities和StrengthLossyesBasicRelationship中定義鋼筋的彈性模量、極限應(yīng)力和應(yīng)變（與本構(gòu)關(guān)系圖對應(yīng)，可查規(guī)范。----宋停云與您分享--------宋停云與您分享----StrengthLoss中定義本構(gòu)關(guān)系中的拉應(yīng)變，具體可參算例?？梢詴翰豢紤]。StrainCapacities中定義三水準(zhǔn)鋼筋的應(yīng)變，具體可參算例。Material中選彈性剪切材料“ElasticShearMaterialforawall”new菜單，輸入代表剪切材料的名稱“ShearMaterial”ok，在右側(cè)菜單只需要填（非線性剪切材料定義后費(fèi)機(jī)時，不易收斂）CrossSects中選“ShearWallInelasticSectionnew菜單，輸入“Sect_組名”AutoSizeok，在右側(cè)菜單中直接選已經(jīng)定義好的混凝土和鋼筋材料作為剪力墻纖維束，同時填寫墻厚度、纖維數(shù)量（4~8為宜，模33個纖維數(shù)，以簡便考慮邊緣約束構(gòu)件的作用。Compound中選“ShearWallCompoundComponent”new菜單，輸入“Com_組名”并ok，右側(cè)填寫基本信息即可。4）定義連梁單元（關(guān)鍵是定義一個剪切鉸，但沒有定型的剪切骨架曲線：Inelastic中選“ShearHinge,DisplacementnewShearHingeokUseCrossSectioDeformationCapabilitiesStrengthLossShapeofRelationshipE-P-P。右側(cè)需填寫三個菜單，具體可參算例。Compound中選“FrameMemberCompoundComponentnew菜單，輸入“Com_組名”ok，右側(cè)填寫基本信息即可，特別需要在中間加一個剪力鉸。（5）定義彈性梁：方法同上，只需要定義其彈性特性。perform下建模菜單中的“Elements”將裝配好的特性賦予給構(gòu)件組使其形成單元。先選組，然后在Properties中點(diǎn)擊“AssignComponent”按鈕，選取裝配類型ChoosecomponenttypeAssign按鈕完成特性賦值。將各個構(gòu)件定義局部坐標(biāo)系，否則后續(xù)不給計(jì)算。perform下建模菜單中的“Driftsanddeflections”訂制層間位移和變形。（1）訂制層間變形及頂部到底部的監(jiān)測變形；----宋停