SARS_病毒传播的数学模型.doc

2010/2011第一学期数学建模课程设计作业 题 目：SARS 病毒传播的数学模型姓 名： 班 级： 专 业： 系 部： 理学系 2011年01月07摘 要 为研究传染病的传播规律、预测和控制传染病蔓延创造条件的重要性，我借助微分方程和曲线拟合等方面的知识，建立两个模型：确诊病例比例变化模型、疑似病例比例变化的模型。在模型求解中，我确定了几个重要参数的值：病人的日接触率、疑似病人的日接触率、日治愈率、日转化率和日死亡率。此外，我还用相轨线分析了病人及健康人群在总群体中的比例变化关系，并运用matlab软件拟合曲线直观地说明情况。由于模型中的参数都具有随机性，所以得出的结果还是会有误差的存在，但所建立的模型对控制疫情的发展有一定的借鉴意义。 关键词：、微分方程、matlab软件、曲线拟合问题阐述是21世纪第一个在世界范围内传播的传染病。SARS的爆发和蔓延给我国的经济发展和人民生活带来了很大影响，我们从中得到了许多重要的经验和教训，认识到定量地研究传染病的传播规律、为预测和控制传染病蔓延创造条件的重要性。为此，我们提供了相关数据供研究病毒的传播规律及其对经济某个方面的影响之用。模型假设(1) 山西省的总人口保持不变，且不考虑人口的迁入与迁出；(2) 人群分为：健康者、确诊病人和退出者（包括死亡者和治愈者）；(3) 在政府的强制措施下，人口基本不流动，故无病源的迁出，避免了交叉感染，降低了感染率,即:被隔离的人与外界的联系几乎全部中断，不具有传染性；(4) 由于对病原体的研究不够深入，加之属于突发型，认为人类在暴发前对无免疫性；(5) 治愈或死亡者均为确诊病人，疑似病人一旦被排除为病人则于第二天进入健康者之列。符号定义及说明：时刻健康者在总人群中所占的比例；：时刻病人在总人群中所占的比例；：时刻疑似病人在总人群中所占的比例；：时刻退出者在总人群中所占的比例；：病人对人群中健康者的日接触率，即每个病人每天有效接触（足以使之致病）的健康者平均人数；：日治愈率，即每天被治愈的病人人数占当天病人总数的比例；：日转化率，即每天疑似病人被确诊为病人人数占当天疑似病人总数的比例；：日死亡率，即每天病人死亡的数量占当天病人总数的比例；：疑似日感染率，即每天每个疑似病人使人群中的健康者感染为疑似病人的平均人数；:整个传染期内每个病人有效接触的平均人数,称为接触数。 模型分析与建立 由分析知，和的变化趋势是预测的传播规律的关键。因此，就它们建立相应的两个模型，健康者的比例分别用和表示。模型一 :病人的变化情况 由假设及符号定义说明，每个病人每天可使个健康者感染为病人，而病人总数为，因此，每天将有个健康者感染为病人。于是，就是病人数的增加率，又每天将有个病人被治愈，个病人死亡，从而得到病人的变化满足： （1）且有 （2）退出者变化情况满足 （3）且 （4）此外，令初始时刻，；将得到病人变化情况的模型的微分方程形式： （5）模型二 :疑似病人的变化情况类似于模型一的分析，得到模型二的微分方程形式： （6）模型求解 所建立的模型一和模型二都是常微分方程模型，所以将视参数、为待定常数。另外，由于和未知，无法求出和的解析解。但可以作数值计算，利用处理已知数据获得、和的合理的最佳估计值；再借助编程得到和的合理的最佳估计值。(一)参数的分析和确定 1. 、和的确定由假设得，用处理已知数据得到： ，在用处理已知数据时，发现了一些反常数据,现假设，在暴发前期，由于人类知之甚少，导致病人被误诊，而在后期，随着人类对的深入探究，误诊最终被更正。而实际的情形是，反常数据出现在后期，这就说明了假设的合理性。(处理过程见附件1)2.和的确定(1) 确定 根据附录1的处理情况,能够作出(或)的实际图像,而通过作数值解,可获得的图像。显然，它们是相似的，这就给确定指明了一种有效的途径：即：使的图像尽可能地吻合的图像，而最佳的吻合效果所对应的即为所求。的图像如下图1所示（程序见附录2）：图1以下将用龙格-可库塔法作的数值图像，(程序见附件2)。根据实际情况，为了方便计算，令山西2010年的总人口为万，取初值，结果如下图2所示：图2此时,。图1和图2的高峰期到平稳期十分吻合，但高峰期前的情况却有较大差别，这主要是缺乏潜伏期数据造成的。若获得潜伏期的实际数据，图1和图2仍将是很相似的。(2)用类似于（1）确定的方法，得到的最佳值为。模拟结果如图3和图4所示： 图3图4(二)相轨线分析和 平面称为相平面,相轨线在相平面上的定义域,为 (7)在方程（5）中消去并注意到的定义，可得 （8）为方便处理，令 （9）从而 （10）容易求出方程（10）的解为 （11）令，于是（11）化为 （12）在定义域内，（12）式表示的曲线即为相轨线，如图5所示。 图5 从图5中可以得到一些初步的信息，即随着时间的增加先增大后减少，但并未能肯定最终相交，即是否必定消失，这是由于把所用各参数均定为常数的缘故。实际上，随着疫情的发展，各参数都应是变化的，集中反映为的变化。而显然是递增的。图6即是取不同的值时所反映的情况比较。 图6其中红线、绿线、蓝线分别表示随着取值的增大（），和的发展情况。显然从图6中，可以肯定必定消失，即。由（5）式，令，此时，将达到其最大值，而将等于。 接下来，讨论：设最终未被感染的健康者的比例为，在（11）式中令得到: （13）而正是（13）式在内的根。在图形上是相轨线与轴在内交点的横坐标。 由图5和图6得到以下推断：（1） 若，则先增加，当时，达到最大值 （14）然后减小且趋于零，则单调减小至。（2） 若，则单调减小至零，单调减小至。 容易理解，只有当递增的阶段才被认为是传染病在蔓延，那么就相当于一个阈值，当时，传染病就会蔓延，而当时，传染病不会蔓延。就而言，在传播初期，由于政府和大众缺乏重视，此时的日接触率偏高，而日治愈率偏低导致较小，不可避免的引发的蔓延。于是，得到了有关群体免疫和预防的看法：（1） 提高医疗卫生水平和群众的警惕意识，使阈值变大，从而阻止等传染病的蔓延；（2） 通过有效的预防接种使群体免疫，降低。 如果忽略病人比例的初始值，则有。于是传染病不会蔓延的条件可表示为 （15）也就是说，只要通过群体免疫使初始时刻的移出者比例（即死亡者和治愈者比例之和）满足（15）式，就可以制止等传染病的蔓延。 但是途径（2）却是难以做到的，它要求免疫者要均匀分布在全体人口中。另外，面对诸如等突然暴发的新传染病而言，由于没有任何医学记录，是无法进行群体免疫的。所以，在任何时候，提高医疗卫生水平及人们的警惕意识都是极为必要的。但无论如何，针对通常的传染病，通过接种免疫来降低传染病蔓延的风险是行之有效的。(三)数值验证与估量 根据上面的分析，制止等传染病的蔓延有两种途径：一是提高医疗卫生水平及大众的警惕意识，即降低日接触率，提高日治愈率（而假设日死亡率不变）；二是群体免疫，即提高移出者的比例初值（相当于降低健康者比例的初值）。下面通过数值计算来验证并估量这两种途径（假设人类已经能够进行免疫接种）的效果。不妨用最终未感染的健康者的比例和病人比例的最大值作为等传染病蔓延程度的度量指标。我们将给定不同的，用（13）式计算，用（14）式计算。结果如下表1所示：途径一1.00.30.30.980.020.03990.34490.60.30.50.980.020.19650.16350.50.51.00.980.020.81220.02020.40.51.250.980.020.91720.0542途径二1.00.30.30.70.020.08400.16850.60.30.50.70.020.30560.05180.50.51.00.70.020.65780.07670.40.51.250.70.020.67550.1948由表1，分析得到：（1） 对于一定的，降低，提高，会使变大，变小；（2） 对于一定的，降低（即提高），也会使变大，变小；（1）和（2）是当时的情况，它们表明提高医疗卫生水平和人们的警惕性及进行群体免疫确实对控制等传染病的蔓延起到了积极有效的作用。（3） 当时，从的结果看，无论怎样，都相应变小了；且对于一定的，降低，提高，和对于一定的，降低，都会缩短等传染病的传染期（即使得的值变小）缩短；从而得到的结论和（1）、（2）相同。在（13）式中略去很小的，即有 （16）当等传染病再次到来时，如果估计没有多大变化，那么就可以用上面得到的分析这次等传染病的蔓延过程。(四)模型验证由（2）、（3）、（5）式解得 （17） （18）当时，取（18）式右端展开式的前3项，在初值下的解为 (19)其中.将（17）、（19）式代入（2）式得 （20）代入前面确定的各参数的值，就可以预估疫情的发展趋势了。对卫生部的措施的评估 仍就模型一进行讨论，对于卫生部采取的隔离措施，首先影响的是日接触率的值。显然，若隔离措施提前实施，则变小；如若延迟隔离，则增大。为了直观地表现不同隔离措施的效用，分别取，并作了图7（程序见附录3）：图7从图7知，若提前采取隔离措施，将延缓传播高峰期的到来，从而为疫情控制争取到了更多的时间；若延后隔离措施，将加速的传播，从而增加了疫情控制的难度，所以建议提前采取隔离措施。至于卫生部采取的提前或延后的时间定为天，有一定的道理。因为由前面确定的参数知天，即平均传染期为天，而高峰期在第天左右。显然提前天采取隔离措施是合理的。模型评价本模型中采用微分方程的模型，对疫情传播做出合理假设，并对其得过拟合，进一步研究和控制疫情的发展，但模型中的参数都具有随机性，所以得出的结果还是会有误差的存在，不能准确预报每次传染病高潮的到来的准确时间，只能限定在一定时间内。不过，所建立的模型对控制疫情的发展有一定的借鉴意义。 参考文献1.数学模型（第三版）姜启源 谢金星 叶俊 编 高等教育出版社 2003年版2.应用教程张圣勤 编 机械工业出版社 2009年版3.传播的数学模型 来源：4.传染病传播及预防的数学模型 来源：附录1的数据处理过程日 期已确诊病例累计现有疑似病例死亡累计治愈出院累计当天病人总数当天死亡总数当天确诊病人数当天治愈出院总数日治愈率日转化率日死亡率4月20日3394021833-4月21日48261025434317143100.0232020.2344260.0162414月22日5886662846520310630.0057690.1591590.0057694月23日6937823555619710590.014540.1342710.0113094月24日774863396468448190.0131580.0938590.0058484月25日8779544273774310390.0116280.1079660.0038764月26日98810934876873611130.0034360.1015550.0068734月27日111412555678990812620.002020.1003980.0080814月28日1199127559781065385000.0666670.0028174月29日1347135866831210714850.0041320.1089840.0057854月30日144014087590129199370.0054220.0660510.0069715月1日155314158210013887113100.0072050.0798590.0050435月2日1636146891109145498390.006190.056540.006195月3日17411493961151541510560.0038940.0703280.0032455月4日18031537100118159246230.0018840.0403380.0025135月5日18971510103121167939430.0017870.0622520.0017875月6日196015231071341736463130.0074880.0413660.0023045月7日20491514110141180838970.0038720.0587850.0016595月8日213614861121521885287110.0058360.0585460.0010615月9日217714251141681913241160.0083640.0287720.0010455月10日22271397116175194525070.0035990.0357910.0010285月11日226514111201861974438110.0055720.0269310.0020265月12日230413781292081998939220.0110110.0283020.0045055月13日234713381342442010543360.017910.0321380.0024885月14日23701308139252199252380.0040160.0175840.002515月15日23881317140257199711850.0025040.0136670.0005015月16日240512651412732008117160.0079680.0134390.0004985月17日242012501453072006415340.0169490.0120.0019945月18日243412501473321982214250.0126140.01120.0010095月19日24371249150349195833170.0086820.0024020.0015325月20日24441225154395194547460.023650.0057140.0020575月21日24441221156447189520520.02744100.0010555月22日245612051585281853212810.0437130.0099590.0010795月23日24651179160582177929540.0303540.0076340.0011245月24日249011341636671748325850.0486270.0220460.0017165月25日24991105167704166949370.0221690.0081450.0023975月26日25041069168747163315430.0263320.0046770.0006125月27日25121005172828159748810.050720.007960.0025055月28日2514941175866151432380.0250990.0021250.0019825月29日2517803176928147613620.0420050.0037360.0006785月30日25207601771006141613780.0550850.0039470.0007065月31日25217471811087133841810.0605380.0013390.002996月1日25227391811124125401370.0295060.00135306月2日25227341811157121700330.027116006月3日25227241811189118400320.027027006月4日25227181811263115200740.064236006月5日25227161811321107800580.053803006月6日25227131831403102020820.08039200.0019616月7日2523668183144693701430.0458910.00149706月8日252255018415438931-1970.108623-0.001820.001126月9日25224511841653795001100.138365006月10日2522351186174768520940.13722600.002926月11日2523257186182159001740.1254240.00389106月12日2523155187187651610550.10658900.0019386月13日25227118719444590-1680.148148-0.0140806月14日25224189199439120500.12787700.0051156月15日25223189201533900210.061947006月16日2521319020533171-1380.119874-0.333330.0031556月17日25215190212027800670.241007006月18日25214191215421110340.16113700.0047396月19日25213191217117600170.096591006月20日25213191218915900180.113208006月21日25212191223114100420.297872006月22日2521219122579900260.262626006月23日2521219122777300200.27397300均值：0.0550760.0249740.002443附录2模型一中图1的程序：x=0:63;y=431,520,619,684,774,873,990,1065,1210,1291,1388,1454,1541,1592,1679,1736,1808,1885,1913,1945,1974,1998,2010,1992,1997,2008,2006,1982,1958,1945,1895,1853,1779,1748,1669,1633,1597,1514,1476,1416,1338,1254,1217,1184,1152,1078,1020,937,893,795,685,590,516,459,391,339,317,278,211,176,159,141,99,73;p s=polyfit(x,y,7)y1=polyval(p,x)plot(x,y,or,x,y1,-k)grid onxlabel(t)ylabel(i)title(每天病人总数的变化情况)legend(实际数据,-拟合曲线)图2的程序：function y=ill2(t,x)a=0.56;b=0.055076;c=0.002443;y=a.*x(1).*x(2)-b.*x(1)-c.*x(1),-a.*x(1).*x(2);ts=0:0.001:63;x0=431/14000000,1-431/14000000;t,x=ode45(ill2,ts,x0);t,x;plot(t,x(:,1)grid onxlabel(t)ylabel(i)title(数值解病人的变化情况)模型二中图3的程序：x=0:63;y=402,610,666,782,863,954,1093,1255,1275,1358,1408,1415,1468,1493,1537,1510,1523,1514,1486,1425,1397,1411,1378,1338,1308,1317,1265,1250,1250,1249,1225,1221,1205,1179,1134,1105,1069,1005,941,803,760,747,739,734,724,718,716,713,668,550,451,351,257,155,71,4,3,3,5,4,3,3,2,2,2;p s=polyfit(x,y,5)y1=polyval(p,x)plot(x,y,or,x,y1,-k)grid onxlabel(t)ylabel(i)title(疑似病人总数的变化情况)legend(实际数据,-拟合曲线)图4的程序：function y=ill2(t,x)a=0.56;b=0.024974;y=a.*x(1).*x(2)-b.*x(1),-a.*x(1).*x(2);ts=0:0.001:63;x0=431/34000000,1-431/34000000;t,x=ode45(ill2,ts,x0);t,x;plot(t,x(:,1)grid onxlabel(t)ylabel(i)title(疑似病人数值解的变化情况)附录3（1） is相轨线程