学院： 专 业： 姓 名：
学号：___ ____ 实验时间：__ ____ 实验地点：
一、实验项目：传染病模型求解
二、实验目的和要求
a.求解微分方程的解析解
b.求解微分方程的数值解
三、实验内容
问题的描述
各种传染病给人类带来的巨大的灾难，长期以来，建立传染病的数学模型来描述传染病的的传播过程，分析受感染人数的变化规律，探索制止传染病蔓延的手段等，一直是各国有关专家和关注的课题。
不同类型传染病有各自不同的特点，在此以一般的传播机理建立几种3模型。分别对3种建立成功的模型进行模型分析，便可以了解到该传染病在人类间传播的大概情况。
模型一（SI模型）：
（1）模型假设
1.在疾病传播期内所考察地区的总人数N不变，人群分为健康人和病人，时刻t这两类人在总人数中所占比例为s（t）和i（t）。
2.每个病人每天有效接触的平均人数是常数a，a成为日接触率，当病人与健康者有效接触时，可使其患病。
（2）建立模型
根据假设，每个病人每天可使as（t）个健康人变成病人，t时刻病人数为Ni（t），所以每天共有aNs（t）i（t）个健康者被感染，即病人的增加率为： Ndi/dt=aNsi
又因为s（t）+i（t）=1
再记时刻t=0时病人的比例为i0
则建立好的模型为：
i(0)=i0
（3）模型求解 （代码、计算结果或输出结果）
syms a i t i0 % a：日接触率，i：病人比例， s：健康人比例，i0：病人比例在t=0时的值
i=dsolve('Di=a*i*(1-i)','i(0)=i0','t');
y=subs(i,{a,i0},{0.3,0.02});
ezplot(y,[0,100])
figure
i=str2double(i);
i=0:0.01:1;
y=0.3*i.*(1-i);
plot(i,y)
SI模型的i~t曲线 SI 模型的di/dt~i 曲线
（4）结果分析
由上图可知，在i=0:1内，di/dt总是增大的，且在i=0.5时，取到最大值，即在t->inf时，所有人都将患病。
上述模型显然不符合实际，为修正上述结果，我们重新考虑模型假设，建立SIS 模型
模型二（SIS模型）
（1）模型假设
假设条件1.2与SI模型相同；
3.每天被治愈的病人数占病人总数的比例为常数u，成为日治愈率，病人治愈后成为仍可被感染的健康者。显然1/u是平均传染期。
（2）模型建立
病人的增加率：Ndi/dt=aNsi-uNi 且 i（t）+s(t)=1；
则有： di/dt=ai(1-i)-ui
在此定义k=a/b，可知k是整个传染传染期内每个病人有效接触的平均人数，成为接触数。
则建立好的模型为：
i(0)=i0;
（2）模型求解 （代码、计算结果或输出结果）
>> syms a i u t i0 % a：日接触率，i：病人比例，u：日治愈率，i0：病人比例在t=0时的值
>> dsolve('Di=a*i*(1-i)-u*i','i(0)=i0','t') % 求用u表示的i—t解析式
>> syms k % k：接触数
>> k=a/u;
>> i=dsolve('Di=-a*i*i+a*i*(1-1/k)','i(0)=i0','t') % 求用k表示的i—t解析式
% 给k、a、i0指定特殊值，作出相关图像
>> y=subs(i,{k,a,i0},{2,0.3,0.02}); %①k>1的情况，以k=2为例
>> ezplot(y,[0,100])
>>pause %作i—t图，分析随时间t的增加, i的变化
>> gtext('1/k')
>>legend('k>1 本例中k=2')
>>figure
>> i=str2double(i);
>> i=0:0.01:1;
>> y=-0.3*i.*[i-1/2];
>> plot(i,y) %作di/dt—i的图像
>> gtext('1-1/k,在此图中为0.5')
>> legend('k=2')
>> y=subs(i,{k,a,i0},{0.8,0.3,0.02}); %②k<1 的情况，以k=0.8为例
>> ezplot(y,[0,100]) %作i—t图，分析随时间t增加，i的变化
>> legend('k<1 本例中k=0.8')
>>figure
>> i=str2double(i);
>> i=0:0.01:1;
>> y=-0.3*i.*[i-(1-1/0.8)];
>> plot(i,y) %作di/dt—i 的图像
>> legend('k=0.8')
>> gtext('k<=1时的情况)
SIS 模型的di/dt—i曲线 （k>1） SIS模型的i—t曲线（k>1）
SIS 模型的di/dt—i曲线 （k<1） SIS模型的i—t曲线（k<1）
（4）结果分析
不难看出，接触数k=1是一个阈值，当k>1时，i（t）的增减性取决于i0的大小，但其极限值i(∞)=1-1/k随k的增加而增加；当k<=1时，病人比例i（t）越来越小，最终趋于0，这是由于传染期内经有效解除从而使健康者变为的病人数不超过原来病人数的缘故。
模型三．SIR模型
（1）模型假设
1.总人数N不变，人群分为健康者、病人和病愈免疫的移出者三类，称SIR模型。时刻三类人在总人数N中占得比例分别记作和。
2.病人的日接触率为，日治愈率为（与SI模型相同），传染期接触数为。
（2）模型建立
由假设1显然有
（1）
对于病愈免疫的移出者而言应有
（2）
再记初始时刻的健康者和病人的比例分别是s0(s0>0)和i0(i0>0)（不妨设移出者的初始值r0=0），则SIR模型的方程可以写作
（3）
(3)模型求解
我们无法求出解析解，先做数值计算：
设，用MATLAB软件编程：
function y=ill（t，x）
a=1;b=0.3;
y=[a*x(1)*x(2)-b*x(1), -a*x(1)*x(2)]';
ts=0:50;
x0=[0.02,0.98];
[t,x]=ode45('i11',ts,x0);[t,x]
plot(t,x(:,1),t,x(:,2)),grid,pause
plot(x(:,2),x(:,1))
表1 的数值计算结果
| t | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
| i(t) | 0.0200 | 0.0390 | 0.0732 | 0.1285 | 0.2033 | 0.2795 | 0.3312 | 0.3444 | 0.3247 |
| s(t) | 0.9800 | 0.9525 | 0.9019 | 0.8169 | 0.6027 | 0.5438 | 0.3995 | 0.2839 | 0.2027 |
| t | 9 | 10 | 15 | 20 | 25 | 30 | 35 | 40 | 45 |
| i(t) | 0.2863 | 0.2418 | 0.0787 | 0.0223 | 0.0061 | 0.0017 | 0.0005 | 0.0001 | 0 |
| s(t) | 0.1493 | 0.1145 | 0.0543 | 0.0434 | 0.0408 | 0.0401 | 0.0399 | 0.0399 | 0.0398 |
的图形 i—s图形（相轨线）
（4）结果分析
的图形见左图， 的图形见右图，称为相轨线，随着的增加，沿轨线自右向左运动。由上图结合表1可知，由初值增长至约时达到最大值，然后减少，则单调减少。
进行相轨线分析，可得：
平面称为相平面，相轨线在相平面上的定义域为
在方程（3）中消去，并注意到的定义，可得
， （4）
容易求出它的解为
（5）
在定义域D内，上式表示的曲线即为相轨线
1.不论初始条件如何，病人终将消失，即
（6）
其证明如下，首先，由（3），而故存在；由（2），，而，故存在，再由（1），对于充分大的有，这将导致，与存在相矛盾。
2.最终未被感染的健康者的比例是，在（5）式中令得到，是方程
（7）
在内的根。在图形上，是相轨线与轴在内交点的横坐标。
3.若，则先增加，当时，达到最大值
（8）
然后减小且趋近于0，则单调减小至。
4.若，则单调减少至0，单调减少至。如果仅当病人比例有一段增长的时期才认为传染病在蔓延，那么是一个阈值，当（即）时传染病就会蔓延。而减小传染期接触数，即提高阈值，使得（即），传染病就不会蔓延（健康者比例的初始值是一定的，通常可认为接近1）。
并且，即使，从（7），（8）式可以看出，减少时，增加（通过作图分析），降低，也控制了蔓延的程度，我们注意到，在中，人们的卫生水平越高，日接触率越小；医疗水平越高，日治愈率越大，于是越小，所以提高卫生水平和医疗水平有助于控制传染病的蔓延。
从另一方面看，是传染期内一个病人传染的健康者的平均数，称为交换数，其含义是一个病人被个健康者交换，所以当即时，必有，既然交换数不超过1，病人比例绝不会增加，传染病不会蔓延。
建模所得：
1.符号变量如何使用
2.如何求微分方程的解析解和数值解
3.对符号变量方程作图时，先将其中的符号变量赋值，再将其变成数值变量，这也是一种有效的解决方法。下载本文
