微分方程建模（传染病模型）的求解。
1、模型1：SI模型。 假设：
（1）t时刻人群分为易感者（占总人数比例的s(t)）和已感染者（占总人数比例的y(t)） （2）每个病人每天有效接触的平均人数是常数?，?称为日接触率，当健康者与病人接触时，健康者受感染成为病人。
分析：根据假设，每个患者每天可以使?s(t)个健康者变为病人，因为病人数为Ny(t)，所以每天共有?Ns(t)y(t)个健康者变为病人。即：
Ndy??Nsy，且s(t)?y(t)?1，设初始时刻病人比例为b，则： dt?dy???y(1?y)，用MATLAB解此微分方程： ?dt??y(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 ?b?1??t11?e1?(?1)e??tbb当b?0.09,??0.1时，分别在坐标系oty中作出y(t)的图像，坐标系oyy?中作出
y???y(1?y)的图像，
>> a=0.1;
>> b=0.09;
>> h=dsolve('Dy=a*y*(1-y)','y(0)=b','t') h =
1/(1-exp(-a*t)*(-1+b)/b) >> f=subs(h) f =
1/(1+91/9*exp(-1/10*t)) >> ezplot(f,[0,60]) >> grid on
>> figure (2)
>> fplot('0.1*y*(1-y)',[0,1]) >> grid on
0.02510.90.80.70.60.50.40.30.20.100101/(1+91/9 exp(-1/10 t))2030t405060y(t)的图像
0.020.0150.010.005000.10.20.30.40.50.60.70.80.91
y???y(1?y)的图像
模型分析：（1）当y?1dy时，达到最大值，则此时病人增速最快。 2dt （2）当t??时，y(t)?1，即所有的人被传染，全部变为病人，这显然是不符合实际的，其原因是没有考虑到病人可以治愈，人群中的健康者只能变为病人，而病人不
会变为健康者。
2、模型2：SIS模型。 假设：
（1）t时刻人群分为易感者（占总人数比例的s(t)）和已感染者（占总人数比例的y(t)） （2）每个病人每天有效接触的平均人数是常数?，?称为日接触率，当健康者与病人接触时，健康者受感染成为病人。
（3）病人每天被治愈的占病人总数的比例为?，称为日治愈率，显然
1为这种传染病的?平均传染期。则Ndy??Nsy??Ny。则建立微分方程模型为： dt?dy???y(1?y)??y ?dt?y(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) / 化简：
a?c
e?(a?c)t.(?a?c?ba)ae?(a?c)t.(?a?c?ba)ca??b(a?c)b(a?c)b(a?c)2? ab(a?c)?e?(a?c)t.(?a?c?ba)?e?(a?c)t.(?a?c?ba)cb(a?c)2?ab(a?c)?(c?a)?e?(a?c)t.(?a?c?ba)?ab(a?c)?(c?a)?e?(a?c)t.(?a?c?ba)???? 2b(a?c)???11a?a????(?)e?(a?c)t? ?a?cba?c????1?即：y(t)???(?)e?(???)t?。
????b???????1?当（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?1?11??即：y(t)???t??。
b??定义???1?：一个传染期内每个病人有效接触的平均人数。 ??1?1?,(??1)则：y(?)???，用MATLAB作图像：
??0(??1)令??0.01，??0.05，b?0.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.100204060t80100120