1、数学实验报告实验名称 数学建模与 MATLAB 学 院 材料学院 专业班级 材料 1014 姓 名 徐萌 孔德成 戴思雨 学 号 41071046 41030400 41030399 2012 年 6 月21、问题的提出。传染病是当今世界最严重的疾病之一，2009 年 4 月 26 日世界卫生组织以确认，美国和墨西哥发生了甲型 H1N1 流感，随后疫情迅速蔓延，截止 8 月中旬，全球感染人数约 5 万人。因此，运用传染病的数学模型来描述传染病甲型H1N1 流感的传播过程，分析受感染人数的变化规律，探索制止甲型 H1N1 蔓延的手段是值得关注的。2、模型的建立。考查中国内地疫情变化，在疾病传播期
2、间不考虑人口的出生率和死亡率，人口总数不变，为常量。中国的疫情研究发现易感染人数多为2050岁的青壮年，故保守估计在此传染病系统的人数N=50000人。甲型HINI流感的传播途径是与病源的直接接触，患者与健康者接触时，都使健康者感染病变故将人群分为3类：健康者(易感染者人群)、患者(已被感染人群)、治愈者(研究期间6月14日8月14日间中国内地感染病毒死亡人数为0，故此处不考虑死亡者)三者在总人数中的比例分别为 :s(t),i(t),r(t)且s(t)+i(t)+r(t)=1,io,So分别为患者人数，健康人数的比例初始值设每个患者每日感染健康者的平均人数为日感染率，记为j，则j=j 日新增病
3、例数/(j-1)日（累计确诊人数-累计出院人数） ；每日被治愈的患者人数占其总数的比例为日治愈率，记为 j，则j=j日被治愈的人数/j日累计确诊病人数;定义整个传染期内每个患者有效接触的平均人数为接触数,由s(t)+i(t)+r(t)=1可知，对于病愈免3疫的治愈者而言应有dr/dt=i,因此考虑SIR传染模型，该模型的方程为di/dt=si-i;ds/dt=-si （1）i(0)=io, s(0)=So;三、模型的求解1、数值运算由于在方程(1)中无法求出s(t)和i(t)的解析解，故先做数值运算据来自中国卫生部网站公布的2009年6月14日8月14日的疫情数据(见表1)包括日累计确诊病例、
4、日累计治愈病例等其中缺失的部分数据，将以通过给定的数据拟合得到. 表1疫情原始数据日期 新增病例 确诊病例累计 治愈累计 新增治愈数6月14日 20 185 736月15日 41 226 86 136月16日 11 237 97 116月17日 27 264 114 176月18日 33 297 135 216月19日 31 328 160 256月20日 28 356 185 256月21日 58 414 199 146月22日 27 441 227 286月23日 49 490 251 246月24日 38 528 275 246月25日 42 570 321 466月26日 48 618
5、 338 176月27日 60 678 373 356月28日 51 729 401 286月29日 37 766 445 446月30日 44 810 496 517月1日 56 866 554 587月2日 49 915 612 5847月3日 45 960 660 487月4日 40 1000 704 447月5日 40 1040 749 457月6日 57 1097 793 447月7日 54 1151 870 777月8日 36 1187 927 577月9日 36 1223 985 587月10日 40 1263 1035 507月11日 39 1302 1085 507月12日
6、26 1328 1110 257月13日 26 1354 1134 247月14日 45 1399 1166 327月15日 45 1444 1197 317月16日 41 1485 1230 337月17日 52 1537 1263 337月18日 44 1581 1293 307月19日 44 1625 1323 307月20日 43 1668 1355 327月21日 52 1720 1404 497月22日 52 1772 1454 507月23日 38 1810 1529 757月24日 42 1852 1604 757月25日 26 1878 1663 597月26日 26 190
7、4 1722 597月27日 26 1930 1781 597月28日 37 1967 1817 367月29日 36 2003 1853 367月30日 43 2046 1883 307月31日 44 2090 1912 298月1日 20 2110 1937 258月2日 21 2131 1962 258月3日 21 2152 1988 268月4日 29 2181 2031 438月5日 29 2210 2074 438月6日 27 2237 2098 248月7日 27 2264 2122 248月8日 28 2292 2137 158月9日 28 2320 2152 158月10日
8、28 2348 2167 158月11日 38 2386 2203 368月12日 39 2425 2240 378月13日 57 2482 2261 218月14日 55 2537 2283 225注：2009年疫情效据见文献8以6月15日为基日，当日累计确诊病例226例，累计出院者86例，故 s(0)=(50000-226+86)/50000=0.9972;I(0)=(226-86)/50000=0.0028;在研究期间，平均日感染率和平均日治愈率由每天相应数据平均求得设计程序为：新增病例A 确诊病例累计B 治愈累计C 新增治愈数DA=41 11 27 33 31 28 58 27 49
9、38 42 48 60 51 37 44 56 49 45 40 40 57 54 36 36 40 39 26 26 45 45 41 52 44 44 43 52 52 38 42 26 26 26 37 36 43 44 20 21 21 29 29 27 27 28 28 28 38 39 57 55B=226 237 264 297 328 356 414 441 490 528 570 618 678 729 766 810 866 915 960 1000 1040 1097 1151 1187 1223 1263 1302 1328 1354 1399 1444 1485 15
10、37 1581 1625 1668 1720 1772 1810 1852 1878 1904 1930 1967 2003 2046 2090 2110 2131 2152 2181 2210 2237 2264 2292 2320 2348 2386 2425 2482 2537C=86 97 114 135 160 185 199 227 251 275 321 338 373 401 445 496 554 612 660 704 749 793 870 927 985 1035 1085 1110 1134 1166 1197 1230 1263 1293 1323 1355 140
11、4 1454 1529 1604 1663 1722 1781 1817 1853 1883 1912 1937 1962 1988 2031 2074 2098 2122 2137 2152 2167 2203 2240 2261 2283D=13 11 17 21 25 25 14 28 24 24 46 17 35 28 44 51 58 58 48 44 45 44 77 57 58 50 50 25 24 32 31 33 33 30 30 32 49 50 75 75 59 59 59 36 36 30 29 25 25 26 43 43 24 24 15 15 15 36 37
12、21 22E=A./(B-C) %日感染率e=sum(E)/61 %平均日感染率F=D./(B-C) %日治愈率f=sum(F)/61 %平均日治愈率运行结果：A =Columns 1 through 1641 11 27 33 31 28 58 27 49 38 42 48 60 51 37 446Columns 17 through 3256 49 45 40 40 57 54 36 36 40 39 26 26 45 45 41Columns 33 through 4852 44 44 43 52 52 38 42 26 26 26 37 36 43 44 20Columns 49 th
13、rough 6121 21 29 29 27 27 28 28 28 38 39 57 55B =Columns 1 through 8226 237 264 297 328 356 414 441Columns 9 through 16490 528 570 618 678 729 766 810Columns 17 through 24866 915 960 1000 1040 1097 1151 1187Columns 25 through 321223 1263 1302 1328 1354 1399 1444 1485Columns 33 through 401537 1581 16
14、25 1668 1720 1772 1810 1852Columns 41 through 481878 1904 1930 1967 2003 2046 2090 2110Columns 49 through 562131 2152 2181 2210 2237 2264 2292 2320Columns 57 through 612348 2386 2425 2482 2537C =Columns 1 through 886 97 114 135 160 185 199 227Columns 9 through 16251 275 321 338 373 401 445 496Column
15、s 17 through 24554 612 660 704 749 793 870 927Columns 25 through 32985 1035 1085 1110 1134 1166 1197 1230Columns 33 through 401263 1293 1323 1355 1404 1454 1529 1604Columns 41 through 481663 1722 1781 1817 1853 1883 1912 1937Columns 49 through 561962 1988 2031 2074 2098 2122 2137 2152Columns 57 thro
16、ugh 612167 2203 2240 2261 2283D =7Columns 1 through 1613 11 17 21 25 25 14 28 24 24 46 17 35 28 44 51Columns 17 through 3258 58 48 44 45 44 77 57 58 50 50 25 24 32 31 33Columns 33 through 4833 30 30 32 49 50 75 75 59 59 59 36 36 30 29 25Columns 49 through 6125 26 43 43 24 24 15 15 15 36 37 21 22e =0
17、.173970451885603 %平均日感染率=0.173970451885603 f =0.164030384929960 %平均日治愈率=0.164030384929960 接触数：=/=0.173970451885603/0.164030384929960=1.060598936958463可得模型方程为：di/dt=0.173970451885603si-0.164030384929960i;ds/dt=-0.173970451885603si ； i(0)=0.9972, s(0)=0.0028; 然后用 Matlab 软件编程，解常微分方程做出患者人数比例 i(t)-时间 t/d
18、 关图，健康者比例 s(t)-时间 t/d 关系图 ，患者人数比例 i-健康者比例 s 图 。采用常微分方程数值求解的方法，预测之后的疫情，并绘制出曲线。设计程序为：M文件function wo=bean(t,x)wo=-0.17397045*x(2） 0;0.17397045*x(2) -0.16403038*x;命令：8t,x=ode45(bean,0,460,0.9972,0.0028)t=Columns 1 through 60 5.3145 10.6290 15.9435 21.2580 32.7580Columns 7 through 1244.2580 55.7580 67.25
19、80 78.7580 90.2580 101.7580Columns 13 through 18113.2580 124.7580 136.2580 147.7580 159.2580 170.7580Columns 19 through 24182.2580 193.7580 205.2580 216.7580 228.2580 239.7580Columns 25 through 30251.2580 262.7580 274.2580 285.7580 297.2580 308.7580Columns 31 through 36320.2580 331.7580 343.2580 354
20、.7580 366.2580 377.7580Columns 37 through 42389.2580 400.7580 412.2580 423.7580 435.2580 441.4435Columns 43 through 45447.6290 453.8145 460.0000x = Columns 1 through 30.9972, 0.0028 0.9946, 0.0020 0.9918, 0.0031Columns 4 through 60.9889, 0.0032 0.9859, 0.0034 0.9790, 0.0036Columns 7 through 90.9717,
21、 0.0039 0.9640, 0.0041 0.9560, 0.0042Columns 10 through 120.9479, 0.0043 0.9398, 0.0043 0.9318, 0.0042Columns 13 through 150.9240, 0.0041 0.9166, 0.0039 0.9096, 0.0037Columns 16 through 180.9031, 0.0034 0.8972, 0.0032 0.8918, 0.0029Columns 17 through 210.8869, 0.0026 0.8826, 0.0023 0.8788, 0.0020Col
22、umns 22 through 240.8755, 0.0018 0.8726, 0.0015 0.8700, 0.0013Columns 25 through 270.8679, 0.0012 0.8660, 0.0010 0.8644, 0.0009Columns 28 through 300.8631, 0.0007 0.8619, 0.0006 0.8609, 0.0005Columns 31 through 330.8601, 0.0004 0.8594, 0.0004 0.8588, 0.0003Columns 34 through 3690.8583, 0.0003 0.8578
23、, 0.0002 0.8575, 0.0002Columns 37 through 390.8572, 0.0002 0.8569, 0.0001 0.8567, 0.0001Columns 40 through 420.8565, 0.0001 0.8564, 0.0001 0.8563,0.0001Columns 43 through 450.8562, 0.0001 0.8562, 0.0001 0.8561, 0.0001 i=x(:,2); plot(t,i,g-) xlabel(时间t/d),ylabel(患者比例i) title(i-t关系图 ) s=x(:,1); plot(t
24、,s,k-) xlabel(时间t/d),ylabel(健康者比例s) title(s-t关系图) plot(s,i,r-) xlabel(健康者比例s),ylabel(患者比例i) title(i-s关系图) y=dsolve(Dy=1/x/1.06058937-1,y(0.9972)=0.0028,x)y=(235717995174701779256942769471657*log(x)/250000000000000000000000000000000-(235717995174701779256942769471657*log(2493/2500)/250000000000000000000000000000000 - x + 1即为i=1/*ln(s/s0)-s+s0+i0=1/1.06059894*ln(s/0.9972)-s+110图一 患者比例 i-时间 t关系图图二 健康者比例 s-时间 t关系图
