传染病模型(微分方程)

实用标准文案微分方程建模（传染病模型）的求解。1、模型 1： SI 模型。假设：（1） t 时刻人群分为易感者（占总人数比例的s(t) ）和已感染者（占总人数比例的y(t ) ）（2）每个病人每天有效接触的平均人数是常数， 称为日接触率，当健康者与病人接触时，健康者受感染成为病人。分析：根据假设，每个患者每天可以使s(t) 个健康者变为病人，因为病人数为Ny (t ) ，所以每天共有Ns(t ) y(t) 个健康者变为病人。即：N dyNsy ，且 s(t ) y(t ) 1，设初始时刻病人比例为b ，则：dtdyy(1 y) ，用 MATLAB解此微分方程：dty(0)b syms a b f=dsolve(Dy=a*y*(1-y),y(0)=b,t)f =1/(1-exp(-a*t)*(-1+b)/b)y(t)11%1 b 1 e t1 ( 1 1)e tbb当 b0.09,0.1 时 ， 分 别 在 坐 标 系 oty 中 作 出y(t) 的 图 像 ， 坐 标 系 oyy中 作 出yy(1y) 的图像，1/(1+91/9 exp(-1/10 t) a=0.1;10.9 b=0.09;0.80.7 h=dsolve(Dy=a*y*(1-y),y(0)=b,t)0.6h =0.50.41/(1-exp(-a*t)*(-1+b)/b)0.30.2 f=subs(h)0.1f =0t01020405060301/(1+91/9*exp(-1/10*t)y(t) 的图像 ezplot(f,0,60)0.025 grid on0.02 figure (2) fplot(0.1*y*(1-y),0,1)0.015 grid on0.010.005精彩文档00.10.20.30.40.50.60.70.80.910实用标准文案模型分析：（ 1）当（ 2）当yy(1y) 的图像y 1 时， dy 达到最大值，则此时病人增速最快。2dtt时， y(t)1，即所有的人被传染，全部变为病人，这显然是不符合实际的， 其原因是没有考虑到病人可以治愈，人群中的健康者只能变为病人，而病人不会变为健康者。2、模型 2： SIS 模型。假设：（1） t 时刻人群分为易感者（占总人数比例的s(t) ）和已感染者（占总人数比例的y(t ) ）（2）每个病人每天有效接触的平均人数是常数，称为日接触率，当健康者与病人接触时，健康者受感染成为病人。（ 3）病人每天被治愈的占病人总数的比例为，称为日治愈率，显然1 为这种传染病的平均传染期。则N dyNsyNy 。则建立微分方程模型为：dtdyy(1y)ydty(0)b用 MATLAB解此微分方程： h2=dsolve(Dy=a*y*(1-y)-c*y,y(0)=b,t)h2 = (a-c)/(a-exp(-(a-c)*t)*(-a+c+b*a)/b/(a-c)*a+exp(-(a-c)*t)*(-a+c+b*a)/b/(a-c)*c) pretty(h2)/ exp(-(a - c) t) (-a + c + b a) a (a - c)/|a - -b (a - c)exp(-(a - c) t) (-a + c + b a) c+ -|b (a - c)/化简：ac(a c )t .( ac ba)ae ( a c)t .( ac ba)ceab( ac)b(ac)精彩文档实用标准文案b( ac)2ab(ac)e (a c)t .(acba)e (a c )t .( a c ba)cb(ac)2ab(ac)(ca)e ( a c)t .(acba)ab(ac)(ca)e (a c) t .(ac1ba)b(ac)21a ( 1 a )e (a c)t a c b a c11即： y(t)() e ()t。b( 11当（ 1）时， y(t )e ()t；b（2）时， clear h2=dsolve(Dy=a*y*(1-y)-a*y,y(0)=b,t)h2 = 1/(a*t+1/b)1即： y(t )tb1。定义：一个传染期内每个病人有效接触的平均人数。1则： y( )1,(1)，用 MATLAB作图像：0(1)令0.01，0.05 ， b0.7 （0.2 1） clear a=0.01;b=0.7;c=0.05; h2=dsolve(Dy=a*y*(1-y)-c*y,y(0)=b,t); h22=subs(h2)h22 =-1/25/(1/100-47/700*exp(1/25*t)-1/25/(1/100-47/700 exp(1/25 t)0.60.50.40.30.20.10020406080100120t精彩文档实用标准文案 ezplot(h22,0,120) grid on令，0.15 ，b 0.7或 b 0.3分别作图（0.3 a=0.3;b=0.7;c=0.15; h2=dsolve(Dy=a*y*(1-y)-c*y,y(0)=b,t); h23=subs(h2)h23 =3/20/(3/10-3/35*exp(-3/20*t) subplot(2,1,1) ezplot(h23,0,25) grid on0.7 b=0.3;0.650.6 h24=subs(h2);0.55 subplot(2,1,2)0.505 ezplot(h24,0,25)grid on0.50.450.40.350.3050.21的图像21 ）3/20/(3/10-3/35 exp(-3/20 t)10152025t3/20/(3/10+1/5 exp(-3/20 t)10152025t21的图像（上面b0.7 ，下面 b0.3）模型分析：（ 1）1时，病人比例越来越少，最终趋于零，这是因为传染期内经有效接触从而使健康者变为病人数不超过原来病人数的缘故。（ 2）11时，病人比例 y(t) 增减性是由 b 来决定， 其极限值 y( ) 1随着的增加而增加。3、模型 3： SIR 模型。假设：（ 1）人群分为健康者，其比例s(t ) 、病人 i(t ) 、病愈免疫的移出者r (t ) 。（ 2）病人的日接触率为，日治愈率为，传染期接触数为。则 s(t)i (t )r (t)1，精彩文档实用标准文案dr对于病愈者而言，y ，dt设初始时刻的健康者和病人的比例为s0 和 y0 ，则建立微分方程模型为：dysyydtdsdtsyy(0)y0 , s(0)s0由于此微分方程组的解析解无法求出，则转为相平面sy上讨论解的性质。相轨线的定义域(s, y)D 应为：D( s, y) s0, y0, sy1 ，由方程组消去dt 并将得：dy1ds1sy s sy00用 matlb 求解： dsolve(Dy=1/cma/s-1,y(s0)=y0,s)ans =1/cma*log(s)-s-1/cma*log(s0)+s0+y0 pretty(ans)log(s)log(s0)- - s - - + s0 + y0cmacma即 y(s) ( s0y0 )s1 ln s（相轨线）s0定义域内，1时， ( y0 , s0 ) 分别取(0.3,0.65) ， (0.4,0.35) ， (0.5,0.45)， (0.7,0.25)在同一直角坐标系中作出其图像： cma=1;y0=0.3;s0=0.65; clear f=dsolve(Dy=1/cma/s-1,y(s0)=y0,s); cma=1;y0=0.3;s0=0.65; f1=subs(f); ezplot(f1,0,1) hold on精彩文档1-s10.90.80.70.60.50.40.30.20.1000.10.20.30.40.50.60.70.80.91s实用标准文案 y0=0.4;s0=0.35; f2=subs(f); ezplot(f2,0,1) hold on y0=0.5;s0=0.45; f3=subs(f); ezplot(f3,0,1) hold onSIR 模型的相轨线 y0=0.7;s0=0.25; f4=subs(f); ezplot(f4,0,1) hold on ezplot(1-s,0,1) grid on模型分析：（ 1）不论初始条件s0 ， y0 如何，病人比例越来越少，最终消失。（ 2）最终未被感染的健康者的比例是s ，在 y( s) ( s0 y0 )s1 ln s中。s0令 y(s ) 0 时， ( s0y0 ) s1 ln s0 的单根即为 s ：最终未被感染的健康者的比s0例。在图像上：相轨线与s 轴在 (0, 1 ) 内交点的横坐标。（ 3）当 s01时传染病不会蔓延，（如最左边的曲线， 随着 t 的增加，病人数 y(t)在减小）。所以提高医疗卫生水平（使日接触率减小或者使日治愈率增大），从而使1变大，也可降低s0 （设 y0 0 ，则 s0r0 1），则 s0 1 r01， r0 11，即使免疫者比例增大。这其实是比较困难的。如5 ， r0 80% 。精彩文档