1、关于SARS模型的建立与相关的预测分析摘要：本文先根据材料提供的模型与数据较为扼要地分析了附件1的模型的优缺点，全面地评价了该模型的合理性与实用性。而后在对问题进行较为全面评价的基础上引入更为全面合理的假设和建立系统分析模型。运用联立微分方程组体现疫情发展过程中各类人的内在因果联系，并在此基础上运用经典的龙格库塔微分方程求解算法结合MATLAB编程(程序在附件一)拟合出与实际较为符合的曲线并进行了疫情预测。同时运用双线性函数模型对卫生部的措施进行了评价并给出建议。而后运用差分方程（程序在附件二）就SARS对经济（主要是旅游业）的影响进行了较为准确的分析，进而通过模型算出的理论预测数值与实际数值
2、进行对比，以数值上的显著差异直观地表现了SARS对经济（旅游）的影响，并对接下来的几个月进行了较为合理的预测。本文的最后，通过本次建模过程中的切身体会，以一篇短文评述去说明建立如SARS预测模型之类的传染病预测模型的重要意义。关键词：微分方程 龙格库塔算法 SARS 双线性函数模型 差分方程数学模型(一) 问题的重述SARS（Severe Acute Respiratory Syndrome，严重急性呼吸道综合症, 俗称：非典型肺炎）是21世纪第一个在世界范围内传播的传染病。SARS的爆发和蔓延给我国的经济发展和人民生活带来了很大影响，我们从中得到了许多重要的经验和教训，认识到定量地研究传染病
3、的传播规律、为预测和控制传染病蔓延创造条件的重要性。请你们对SARS 的传播建立数学模型，具体要求如下：（1）对附件1所提供的一个早期的模型，评价其合理性和实用性。（2）建立你们自己的模型，说明为什么优于附件1中的模型；特别要说明怎样才能建立一个真正能够预测以及能为预防和控制提供可靠、足够的信息的模型，这样做的困难在哪里？对于卫生部门所采取的措施做出评论，如：提前或延后5天采取严格的隔离措施，对疫情传播所造成的影响做出估计。附件2提供的数据供参考。（3）收集SARS对经济某个方面影响的数据，建立相应的数学模型并进行预测。附件3提供的数据供参考。（4）给当地报刊写一篇通俗短文，说明建立传染病数学
4、模型的重要性。（二）对附件1所提供的模型的评价该模型的合理性首先体现在模型假设上：“假定初始时刻的病例数为N0，平均每病人每天可传染K个人（K一般为小数），平均每个病人可以直接感染他人的时间为L天。”其一，一般来说每病人每天可传染的人数与当时的健康人数有关1，但由于北京的人数基数较大，SARS病人数相对较少并且SARA持续时间不是很长，所以这样假设也是可以的。其二，每个病人可以直接感染他人的时间是有限的，该模型考虑到了这一点，也是很合理的。该模型的合理性还在于用数理统计的方法估计相关参数。该模型的实用性是较好地模拟与预测了北京的SARA数据与发展。在传染病发病初期对疫情的预测结果还是较为理想的
5、，这主要得益于发病初期，由于病情来得突然，有关部门没有来得及采取措施加以控制，使病情得以蔓延迅速，而且发病初期在治疗方法上不是特别有效，治愈所需的时间长，所以使用作为模型进行估计以及参数的假设均较为合理，基本上是可行的。但是到了疫情发展中后期，由于政府部门采取强硬措施加强防治工作以及人民群众的防范意识与警觉程度上的普遍提高，加之治疗措施的改进，使得每天被传染的人数下降，并且治愈的人数在不断增加，治愈时间也在不断缩短，每天的病人数应在上一天的基础上减去治愈和死亡的人数，并且由于采取强硬措施,“L”的取值会大大的减小，“K”取值也会是个变量，而不是常数。大多数疑是病人往往在早期就会被隔离，所以，基
6、本上很少能转化成自由非典病人而去接触并传染别人。如果此时还是选取这样的单调递增函数作为预测模型，就会有较大的误差。该模型的另一个不足是没有考虑SARS的潜伏期，也没有对人群进行合理的分类（如易感染人群、病人、治愈人群等等），所以必须建立更为合理的假设与模型。（三）定义与符号说明S：表示易感染人群即健康者在人群中的比例。E：处于潜伏期人群在人群众的比率。，这种人暂时未发病，但他们最终将发病。I：已受感染者即病人在人群中的比例。R：移出者(包括“出院者”和“死亡者”)在人群中的比例。M：未被隔离的带菌者。X：疑似病人。a:每个病人每天有效接触并使之感染的平均人数（常数）。b:退出率，即SARS患者
7、的每日死亡率和每日治愈率之和。c：潜伏期的病人的日发病率。d: 每个未被隔离带菌者被隔离前平均每天感染有效人数。x1: 疑似中每日被排除的人数占疑似人数的比例。x2：疑似者中每日确诊的人数占疑似人数的比例。j：每个未被隔离的带菌者转化为病人的日转化率。k：被未被隔离的带菌者有效传染的人中可以控制的比例。0:SARS 潜伏期(天数) 。（四）模型的假设：总体假设：材料提供的疫情统计数据真实有效。基本假设：1、潜伏期一般是2到11天，我们假设为=5天。2、据医学权威表明，传染途径主要是SARS患者，不包括处于潜伏期的。3、SARS患者被治愈后不具有传染性，也不会再被传染。4、北京市的总人数可视为常
8、数，即不考虑流入人口与流出人口的影响。也不考虑这段时间的人口出生率与自然死亡率。把由SARS引起的死亡人数视为“移出者”（R）。5、与SARS病人接触后都会被感染。（五）模型的分析： 根据附件三（我们自己从中国网收集得到的数据，比赛题给的数据要详细）中的数据表明，由于疫情初期政府控制力度不够，大众的对SARS的防范意识不强，造成病情迅速蔓延。而当政府采取有力措施，人们的防患意识增强，疫情则趋于缓和，病患者人数迅速下降。所以SARS传播大体上可分为两个阶段：1控制前期：即认为病毒传播方式是自然传播。2控制后期：政府强力介入之后的病毒传播模型。我们以附件三提供的北京市疫情统计数据为基础建立模型。（
9、六）模型的建立：1、控前模型的建立：将人群分为易感染人群、已受感染者和移出者三类。时间为3月1日到4月19日。记S：表示易感染人群即健康者在人群中的比例。I：已受感染者即病人在人群中的比例。R：移出者(包括“出院者”和“死亡者”)在人群中的比例。a:每个病人每天有效接触并使之感染成为潜伏病人（在天后发病）的平均人数。b:移出率，即SARS患者的每日死亡率和每日治愈率之和。0:SARS 潜伏期(天数) 。则有：初始值：S(0), S(1), S(2), S(3), S(4), I(0) I(1) I(2) I(3) I(4),R(0)该微分方程组带有时滞因素，没有解析解。因此我们考虑用差分方程的
10、办法，因为潜伏期为5天，所以I（0），I（1），I（2），I（3），I（4）分别表示疫情开始前5天的病人数，从第5天开始有：I（5）=I（4）+aS（0）I（0）bI（5）即：I（5）=I（4）+aS（0）I（0）/1+b同理，I（6）=I（5）+aS（1）I（1）/1+b：I（n）=I(n-1)+aS(n-5)I(n-5)/ 1+b (5)对S（t）,亦有S（5）S(4)-aI(0)S(0) S（6）S(5)-aI(1)S(1)： S（n）S(n-1)-aI(n-5)S(n-5) (6)对R（t）,则有R（n）=R(n-1)+bI(n) (n=1,2,) (7)由此可得SARS的控前差分方程
11、模型，即：如果初始值给定，并将参数a，b确定（a：由有关数据推导得出。b：由医疗水平和有关数据分析得出，取其平均值）就可计算出任一天的易感染人群、已受感染者和移出者的数目， 但可惜的是，由于疫情初期政府控制力度不够，没有提供给我们真实有效数据，如北京首例SARS病人出现在3月1日，但只有4月19日以后的数据，所以我们只能进行模型建立和分析，而不能求解模型。这也是建立真正有效的能预测的模型的困难之一。困难之二是这个微分方程组的求解极其困难。困难之三是我们不知道政府在何时干预及力度如何。2、控后模型的建立：将人群分为易感染人群、已受感染者、移出者、疑似病人和未被隔离的带菌者五类。设控制开始时间为4
12、月21日。记S：表示易感染人群即健康者在人群中的比例。I：已受感染者即病人在人群中的比例。M：未被隔离的带菌者。X：疑似病人。R：移出者(包括“出院者”和“死亡者”)在人群中的比例。a:每个病人每天有效接触并使之感染的平均人数（常数）。b:移出率，即SARS患者的每日死亡率和每日治愈率之和。d: 每个未被隔离带菌者被隔离前平均每天感染有效人数。x1: 疑似者中每日被排除的人数占疑似人数的比例。x2：疑似者中每日确诊的人数占疑似人数的比例。j：未被隔离的带菌者转化为病人的日转化率。k：被未被隔离的带菌者有效传染的人中可以控制的比例。则有： S(t)I(t)R(t)X(t)M(t)1 (13) S
13、(0),I(0),R(0),X(0),M(0)为初始值 参数的确定：我们以材料提供的北京市疫情统计数据来说明参数的分析方法。（见附件三）以下全部图的坐标0均表示4月19日。（1）x1：x1=（每天新增的疑似排除人数）/ （当天疑似病人累计人数当天移出累计人数）首先我们先直观的观察一下X1的变化趋势。根据材料提供的数据，用MATLAB来出x1，并画图，如图1所示：图1接着用曲线拟合图1，如图2所示：图2从上图可看出，图2大概有两个峰值。第一个高峰可能是疑似者中非感染者较高；第二个峰值则是因大部分真正带病的疑似者已转化为确诊后，未带菌者相对比例增大造成的。由此4阶拟合得出的曲线误差很大，为此，我们
14、去掉几个偏差太大的点后，易看出，x1集中分布在0到0.05之间。从图中，可以发现，最集中的数据为0.035，这样我们就以0.035为x1的估计值。（2）x2：x2=（每天新增的疑似转化为确诊的人数） / （当天疑似累计人数当天累计移出者）首先观察x2的变化趋势，如图3所示：图3用5阶曲线拟合，如图4所示：图4从图4可见，x2在疫情得到重视后一直下降。由图还可以看出x2的值主要分布在0 .0005和0.015之间。我们根据SARS具有潜伏期的情况，估计x2分为两个阶段值：0.0223和0.006。从对x1与x2的数据处理来看，我们可以将控制后期的这段时期分为两个阶段：过渡期和平稳期；这两个阶段的
15、产生是与非典自身的特性分不开的。由于非典具有潜伏期，所以在控制初期，由于前一段时间对非典的控制力度不够，造成较多的人处于SARS潜伏期，这一部分人最终将转化为SARS病人；且因为他们为未被隔离带菌者，在进医院前会传染较多的人；加之各项措施从颁布到实行总会有一段反应时间，所以上述原因直接导致了过渡期的形成。（3）b：b=（每天新增的出院者和死亡者的人数）/（当天病人累计人数当天累计移出者人数）首先观察b的变化趋势。如图5所示：图5用3阶曲线拟合，如图6所示：图6由图还可以看出B的值主要分布在0 .005和0.09之间。b也可以分为三个阶段值：0.0085，0.036和0.085。（4）j：从材料
16、提供的数据可以估计出其值在0.12到0.25之间。现在我们令j=0.23。（5）k：根据各地的人们的意识和习惯等因素反映出来，比如在控制期对人口流动的控制严格程度，减少聚会等措施。由此我们估计k在0.7到0.9之间。（七）模型的解法与求解： 经过尝试，我们可以知道无法从建立的微分方程直接求出S，I，R，X，M的函数表达式。我们尝试编写龙格库塔微分求解算法求它们的数值解。 龙格库塔算法（求解微分方程组）： K1=f(t,Zn) K2=f(t,Zn+h/2*K1) K3=f(t,Zn+h/2*K2) K4=f(t,Zn+h*K3) Zn+1=Zn+h/2(K1+2K2+2K3+K4)步长h取1Zn
17、为矩阵（S，I，R，X，M）f为5个方程的形式K（K1K4）为（4X5）的矩阵 1由Zn的初值可顺序算出K（K1，K2，K3，K4）的值，进而推出Zn+1的值。据此递推可得任意的Zn的值。（matlab程序附于附件一中）我们通过求出的数值解跟根据材料算出的每一天的S、I、R、X、M的值作对比，并画出图像。根据图像可以看出实际数值和理论计算数值存在一定的差异。因此，我们可以通过调整由估计得来的数据来使两个图像趋于一致。通过分析估计有：x1=0.035，x2=0.02230.006，b=0.00850.1913， d=0.555，j=0.22，k=0.71。这组数据使图像与实际最为接近。图7是材料
18、提供的北京市累计病人比率图： 图7 以下是由龙格库塔方法求得的I的数值解画出的图：（蓝圈表示实际的数据，红线是I的数值解）图8注：0点表示为4月19日，从此开始往后计算天数。从上图曲线的分析可得：1 病人数在5月13号(4月19日往后25天)前后达到高峰期。2 病人数在5月18日左右出现最大值，且Imax=0.00015，此后开始下降。3 病人数在6月10号前后趋于缓解，并不断递减，在6月10号至7月10初为缓解阶段。这是由于政府采取了大力控制措施，使大多数带菌者被隔离，从而被传染的人数越来越少。4 病人数将在130天后即8月底左右趋于0，即疫情的最终控制期。（八）模型的优缺点一模型的优点：(
19、1)首先此模型进行了较为详细的分类，使得微分方程的联立显得较为紧凑，相关性较强。因为材料提供的数据齐全，各参数的设置合理，所以参数的设计在详细的数据的支持下和实际情况较为接近，为图形的拟合和对最终控制期的预测奠定了基础。(2)该模型适用范围较广，只要数据足够地精确详细，则求出估计参数便可求解。(3)该模型虽然是控后模型，但只要人的警惕性增强，把控制期提前到疫情初期，那么只要做适当的参数修改便可以作为疫情发展的全过程的预测模型。二模型的缺点：(1)首先是没有考虑年龄结构层次对疫情的影响。因为根据医学研究表明，儿童与老人极易感染非典病毒，而青壮年由于精力旺盛使的活动积极者，且由于控制后期的前期人们
20、的活动几率仍然较高，所以接触的几率较大。而在控后期的后期由于人们的活动水平降低，使得接触的几率降低，接触几率的不断降低这一点在模型中并没有得到很好的体现。(2)随着医疗水平的提高与人们的意识水平、政府的严厉的控制措施，退出率实际上是不断提高的，而在我们的模型中却认为这是一个常数。（九）控后模型的改进：(1)在模型中引入接触率与移出率应随着病人的减少而变化，可能会是趋于p的变化，随着时间推移有所调整，以附和预测的需要。(2)在模型中引入接触感染率的概念，即体现接触不一定感染，只不过是感染率较高。(3)对控后模型加入潜伏期对病毒的传播造成的影响。(十)对于卫生部门采取的措施的评价对于卫生部门提前或
21、延后5天采取严格的隔离措施的影响，我们可以建立下面的模型进行辅助分析估计：1模型参数定义：S(t)t时刻易感人群总数I(t) t时刻出现的新增患者 患者从患病起经过时间，仍为患者的概率患者距发病时间，具有传染性的概率患者与易感人群接触率近断时间的医学研究表明，从正式发病到治愈一般需714天或更长时间，假定平均治愈时间为12天。2基本条件假设：新患者出现的数量与现有患者的数量成正比，也与现有易感者的数量成正比，即发病率是患者人数和易感者人数的双线性函数。由基本假设条件可得：S（t+1）=S(t)-I(t+1) （1）I(t+1)=S(t)（2）经整理后得：S(t+1)=S(t)-（3）S(t+1
22、)=S(t)(1-（4）S(t+1)/S(t)=(1-（5） 虽然不能具体知道 的数值，那么我可以根据较为理想的均匀平均递减概率参数，可得下表：0123456789101112111/1210/129/128/127/126/125/124/123/122/121/120如果病人发病后5天才开始隔离，并且的值在疫情初期又较大的话，那么由上表可知病人已经分别以11/12、10/12、9/12、8/12、7/12的大概率在社会上与易感人群接触和传染。由（5）式得S(t+1)/S(t)=(1-1 （6）也就是S(t+1) 7/12。当处于潜伏期时，传染性几乎为0，因而同理我们有理由相信：S(t+1)
23、/S(t)=(1- （7）即S(t+1)近似于S(t)。所以，如果在病人发病前提前5天隔离的话，新增病人数将变得很小。（十一）SARS对经济的影响从附件四中2003年度4，5，6，7，8月的旅游人数与往年同期月份的数据作对比，可以发现Sars对经济（旅游业）的影响还是很显著的。下面建立模型进行预测和分析对比，用数值上的差异来说明问题。我们主要用一阶常系数差分方程作为预测模型，他的形式为:x(t+1)+px(t)=q (1)在应用一阶常系数差分方程的过程中，系数p与q是不易确定的，因为旅游的人数多少只对下一年度产生较大的影响，而对后年度的人数不产生较大的影响，所以对p与q的估计使用加权平均的方法
24、来近似计算，即使用1997到2002的旅游人数X(i)进行构造：-P1=X(1998)/X(1997) -P2=X(1999)/X(1998).-P5=X(2002)/X(2001)-P=(1/15)(-P1-2P2-3P3-4P4-5P5) (2)由-p的表达式反映了x(t+1)与x(t)间的加权比例，q则表示x(t+1)与x(t)之间的加权平均值，权数分配与-p的计算方法相同。即：-q=(1/15)(-q1-2(q2)-3(q3)-4(q4)-5(q5) (3)求出p与q具体值之后，可推导出表达式： x(t)=A(-P)+q/(1+p) (p-1) (4) 5其中求A值方法为代入t=0，x
25、(t)为最后一年即2002年时的数值后求出的对应数值，求出A值之后，将t=1代入，便可求出下一年的预测数值。下图9就是我们用所建立的模型的预测数值（红色曲线）与材料（附件四）提供的实际数值（绿色曲线）的对比图：图9从预测出的4,5,6,7,8月的数值与真实数值作对比，可明显看出他们之间的显著差异。所以sars的确对我国经济方面产生不良影响，如旅游，交通，餐饮业等。由于sars疫情逐渐趋于平缓，预计从9月开始旅游业恢复正常，为此我们也给了一些相关预测数据，仅供参考。下图10分别是用所建模型对2003年4月到12月预测数据（红色曲线）与根据材料提供的前6年后8个月的数据用MATLAB拟合的2003
26、年4月到12月预测数据（绿色曲线）的对比图图10月份456789101112模型预测31.567631.452729.809927.744334.437333.659734.895432.845724.8856拟合32.150032.830031.600029.330036.400033.140032.850026.850027.7900实际数11.61.782.618.816.2减少%63.394.391.268.252.9表2从表2可以知，从4月份起，SARS对北京的海外旅游人数打击很大，4月份就骤减了63%左右，当4月底5月初SARS达到高峰时，对人们的心理影响也最大，5月份几乎减少了9
27、5%，6月份因直到月底北京的SARS疫情才解禁，所以该月仍然减少了90%以上，7月与8月份仍心有余悸，分别减少近70%和53%，但影响已逐渐减少。随着SARS的彻底结束和时间的推移，相信从2003年9月至12月北京市的旅游人数将逐渐恢复正常。（十二）通俗短文20世纪末期，全世界医学科学突飞猛进，人类征服疾病的能力提高到了一个崭新的水平。随着人体器官的克隆成功，基因疗法临床应用的突破，为人类“永葆青春”、“延年益寿”展现出美好前景但今春以来，在世界范围有两大热点为人们所关注，这两大热点就是“两萨”，伊拉克战争和非典型肺炎。全球协力“剿萨”。前者指的是萨达姆（Saddam），后者指的则是令人发指的
28、萨斯（SARS）。作为21世纪的瘟神，SARS已经侵袭了世界东西方的33个国家和地区，堪称“东邪西毒”。它横行逞凶，荼毒生灵；变异的冠状病毒兴妖作怪，大展魔法。它袭击我们的呼吸系统，窜入肺腑为非作歹，使人们不能自由呼吸。据5月15日世界卫生组织统计公告，全球感染“非典”病的人数已逾8000，虽各个国家都采取了防范的应对措施，但由于SARS是突发的传染病的变种，使得人们没有及时有效的预测工具作为武器，仍未能在较短的时间内遏制蔓延的势头。“非典”对人类的摧残和毁灭，虽比不上一场大的战争，但波及面之广，影响之大，给人类心灵上造成的负面效应和伤害，远比一场战争深重得多。面对SARS，人类并不是无能为力
29、的。建立非典传染病的数学模型，可以描述传染病的传播过程，分析受感染人数的变化规律，预报传染病高潮的到来等等，这就为人们正确认识SARS，采取必要的预防和隔离措施提供了理论依据。事实上新的传染病疫情发展初期，人们不可能去做传染病传播的试验以获取数据，从医疗卫生部门得到的资料也是不完全与不充分的，所以通常主要是依据以往的机理分析的方法建立模型。造成传染病流行的主要有三个环节，一是传染源，二是传染途径，三是易感人群。但不同类型的传染病的传播过程有各自不同的特点，弄清这些特点需要相当多的病理知识，这里不可能从医学角度一一分析各种传染病的传播，而只是按照以往的传播可能机理建立几种相关模型。正确的判断并提
30、出这些“可能”，只能来自于对突发事件的现场初期的调查研究结果进行预测，因此，现场流行病学调查分析是必不可少甚至是首要的应急措施之一。这样就可以及时的在原有相似模型的基础上建立相应的新模型，对疫情发展作出相关的预测，并可在建立模型的过程中及时发现对疫情发展起关键作用的环节，做到擒贼先擒王，控制疫情发展中枢要害环节，采取强制措施加以控制，并做好应对相关突发情况的准备。而传染病数学模型无疑会给以上所论述的种种带来极大的帮助或起辅助作用。就像此次SARS疫情，如果我们能及时建立可能的预测模型，及时发现隔离疑是病人，SARS病人关键的关键，那么疫情就会的到及时有效的控制。在人类与各种疾病作斗争的时候，别
31、忘了数学模型也是一个强有力的武器！参考文献1姜启源，数学模型（第二版），北京：高等教育出版社，1993年2李庆扬、王能、易大义, 数值分析，武汉：华中理工大学出版社，1986年3寿纪麟，数学建模-方法与范例，西安：西安交大出版社，1993年4 SARS预测模型，2003年8月5 常备不懈，磨砺以须，预防和控制传染病，中华预防医学杂志，37卷第4、5期6 张恩详，来华境外游客的动态预测方法，数学的实践与认识，30卷第1期附录附件一：RungeKutta.my0=482588693 7748779881114119913471440155316361741180318971960204921362
32、1772227226523042347237023882405242024342437244424442456246524902499250425122514251725202521252225222522252225222522252325222522252225232523252225222522252125212521252125212521y1=43465564737678788390100109115118121134141152168175186203244252257273307332349395447528582667704747828866928100610871124115
33、71189126313211403144615431653174718211876194419942015205321202154217121892231y2=25283539424856596675829196100103107110112114116120129134139140141145147150154156158160163167168172175176177181181181181181181183183184184186186187187189189190190191191191191yy=y0-y1-y2; yy=yy/14000000;b=zeros(5,200);k=ze
34、ros(4,5);h=1;b(:,1)=13997526; 414; 68; 610; 600;%初值b(:,1)=b(:,1)/14000000;j=0.48;e=0.23;a=0.71;yy0=0.035;yy1=0.0223;q=0.0085;x=zeros(1,200);Min=1;for i=1:200 x(i)=i;endfor t=1:199 if t=25 yy1=0.006;q=0.036; if t=45 q=0.085; end ends=b(1,t);i=b(2,t);r=b(3,t);y=b(4,t);m=b(5,t);k(1,1)=yy0*y-j*m*s;k(1,2
35、)=e*m-q*i+yy1*y;k(1,3)=q*i;k(1,4)=j*m*s*a-yy0*y-yy1*y;k(1,5)=j*m*s*(1-a)-e*m;for p=1:2 k(p+1,1)=yy0*(y+h/2*k(p,4)-j*(m+h/2*k(p,5)*(s+h/2*k(p,1); k(p+1,2)=e*(m+h/2*k(p,1)-q*(i+h/2*k(p,2)+yy1*(y+h/2*k(p,4); k(p+1,3)=q*(i+h/2*k(p,2); k(p+1,4)=j*(m+h/2*k(p,5)*(s+h/2*k(p,1)*a-(yy0+yy1)*(y+h/2*k(p,4); k(p
36、+1,5)=j*(m+h/2*k(p,5)*(s+h/2*k(p,1)*(1-a)-e*(m+h/2*k(p,5);endk(4,1)=yy0*(y+h*k(3,4)-j*(m+h*k(3,5)*(s+h*k(3,1);k(4,2)=e*(m+h*k(3,5)-q*(i+h*k(3,2)+yy1*(y+h*k(3,4);k(4,3)=q*(i+h*k(3,2);k(4,4)=j*(m+h*k(3,5)*(s+h*k(3,1)*a-(yy0+yy1)*(y+h*k(3,4);k(4,5)=j*(m+h*k(3,5)*(s+h*k(3,1)*(1-a)-e*(m+h*k(3,5);b(:,t+1)
37、=b(:,t)+h/6*(k(1,:)+2*k(2,:)+2*k(3,:)+k(4,:);end plot(x(1:62),yy,o,x,b(2,:),+,x,b(2,:);附件二：economy1.mx=9.4 11.3 16.8 19.8 20.3 18.8 20.9 24.9 24.7 24.3 19.4 18.6;9.6 11.7 15.8 19.9 19.5 17.8 17.8 23.3 21.4 24.5 20.1 15.9;10.1 12.9 17.7 21.0 21.0 20.4 21.9 25.8 29.3 29.8 23.6 16.5;11.4 26.0 19.6 25.9
38、 27.6 24.3 23.0 27.8 27.3 28.5 32.8 18.5;11.5 26.4 20.4 26.1 28.9 28.0 25.2 30.8 28.7 28.1 22.2 20.7;13.7 29.7 23.1 28.9 29.0 27.4 26.0 32.2 31.4 32.6 29.2 22.9a=size(x)p=zeros(1,a(2);q=p;t=p;A=p;flag=0;for i=1:a(1)-1 p=p+i*x(i+1,:)./x(i,:);endp=p/15;for i=1:a(1)-1 q=q+i*(x(i+1,:)-p.*x(i,:);endq=q/1
39、5;A=x(a(1),:)-q./(1+p);for i=1:a(2) if p(i)=-1 flag=1; endendif flag=1t=A.*p+q./(1+p);%预测值endeconomy2.mx=9.4 11.3 16.8 19.8 20.3 18.8 20.9 24.9 24.7 24.3 19.4 18.6;9.6 11.7 15.8 19.9 19.5 17.8 17.8 23.3 21.4 24.5 20.1 15.9;10.1 12.9 17.7 21.0 21.0 20.4 21.9 25.8 29.3 29.8 23.6 16.5;11.4 26.0 19.6 25
40、.9 27.6 24.3 23.0 27.8 27.3 28.5 32.8 18.5;11.5 26.4 20.4 26.1 28.9 28.0 25.2 30.8 28.7 28.1 22.2 20.7;13.7 29.7 23.1 28.9 29.0 27.4 26.0 32.2 31.4 32.6 29.2 22.9;15.4 17.1 23.5 11.6 1.78 2.61 8.8 16.2 0 0 0 0xx=15.4 17.1 23.5 31.5676 31.4527 29.8099 27.7443 34.4373 33.6597 34.8954 32.8457 24.8856t=
41、1 2 3 4 5 6 7 8 9 10 11 12z=1 2 3 4 5 6 7 8fxy1=polyfit(z,x(7,1:8),4);y1=polyval(fxy1,t);plot(z,x(7,1:8),+,t,y1,t,xx);pausefor i=1:12 yy=x(1:6,i); fxy=polyfit(z(1:6),yy,2); f(i)=polyval(fxy,7); y=polyval(fxy,z(1:6); plot(z(1:7)+1996,x(:,i),+,z(1:6)+1996,y,2003,xx(i),ro); pauseendplot(t(4:12),f(4:12)
42、,g,t(4:12),xx(4:12),r);附件三：北京市有关“SARS” 疫情统计数据日期临床诊断病例数其中医务人员数出院人数死亡人数疑似病例数新增/(其中由疑似转为临床诊断数） 累计新增累计新增累计新增累计新增排除合计41969 42169435319610420744821078104342542110658821993 46328666422105693 271269557 3578242389774*17143964439112863424101861151589733421619544251069104162174345167971426113988*71673766481731093427126(50)*111420187278856162119142896(37)*11
