本题可用
x=fminbnd('-0.1*x*(1-x)',0,1)
y=0.1*x*(1-x)
4)指示最大值坐标
用线性绘图函数plot,调用格式如下:
plot(x1,y1,’颜色线型数据点图标’,x2,y2,’颜色线型数据点图标’,…)
说明参见《数学实验》p225
本题可用
holdon;%在上面的同一张图上画线(同坐标系)
plot([0,x],[y,y],':
',[x,x],[0,y],':
');
3)图形的标注
使用文本标注函数text,调用格式如下:
格式1
text(x,y,文本标识内容,’HorizontalAlignment’,’字符串1’)
x,y给定标注文本在图中添加的位置。
’HorizontalAlignment’为水平控制属性,控制文本标识起点位于点(x,y)同一水平线
上。
’字符串1’为水平控制属性值,取三个值之一:
‘left’,点(x,y)位于文本标识的左边。
‘center’,点(x,y)位于文本标识的中心点。
‘right’,点(x,y)位于文本标识的右边。
格式2
text(x,y,文本标识内容,’VerticalAlignment’,’字符串2’)
x,y给定标注文本在图中添加的位置。
’VerticalAlignment’为垂直控制属性,控制文本标识起点位于点(x,y)同一垂直线上。
’字符串1’为垂直控制属性值,取四个值之一:
‘middle’,’top’,’cap’,’baseline’,’bottom’。
(对应位置可在命令窗口应用确
定)
本题可用
text(0,y,'(di/dt)m','VerticalAlignment','bottom');
text(x,-0.001,num2str(x),'HorizontalAlignment','center');
4)坐标轴标注
调用函数xlabel,ylabel和title
本题可用
title('SI模型di/dt~i曲线');
xlabel('i');
ylabel('di/dt');
2.实验7-2传染病模型2(SI模型)——画i~t曲线图
(参考教材p137-138)
传染病模型2(SI模型):
di/dt=ki(1-i),i(0)=i0;
其中,
i(t)是第t天病人在总人数中所占的比例。
k是每个病人每天有效接触的平均人数(日接触率)。
i0是初始时刻(t=0)病人的比例
求出微分方程的解析解i(t),画出如下所示的i~t曲线(i(0)=0.15,k=0.2,
t=0~30)。
试编写一个m文件来实现。
(在图形窗口菜单选择Edit/CopyFigure,复制
图形)
[提示]
1)求解微分方程
常微分方程符号解用函数dsolve,调用格式如下:
dsolve(‘equ1’,’equ2’,…,’变量名’)
以代表微分方程及初始条件的符号方程为输入参数,多个方程或初始条件可在一个输入
变量内联立输入,且以逗号分隔。
默认的独立变量为t,也可把t变为其他的符号变量。
字符D代表对独立变量的微分,通常指d/dt。
本题可用
x=dsolve(‘Dx=k*x*(1-x)’,’x(0)=x0’)
2)画出i~t曲线(i(0)=0.15,λ=0.2,t=0~30)
用for循环,函数length,eval,plot,axis,title,xlabel,ylabel
3.实验7-3传染病模型3(SIS模型)——画di/dt~i曲线图
(参考教材p138-139)
已知传染病模型3(SIS模型):
di/dt=-λi[i-(1-1/σ)],i(0)=i0
其中,
i(t)是第t天病人在总人数中所占的比例。
λ是每个病人每天有效接触的平均人数(日接触率)。
i0是初始时刻(t=0)病人的比例。
σ是整个传染期内每个病人有效接触的平均人数(接触数)。
取λ=0.1,σ=1.5,画出如下所示的di/dt~i曲线图。
试编写一个m文件来实现。
(在图
形窗口菜单选择Edit/CopyFigure,复制图形)
[提示]
用fplot函数画出di/dt~i曲线图;
在上图上用plot函数画一条过原点的水平
用title,xlabel,ylabel标注。
4.实验7-4传染病模型3(SIS模型)——画i~t曲线图
(参考教材p138-139)
已知传染病模型3(SIS模型):
di/dt=-λi[i-(1-1/σ)],i(0)=i0
其中,
i(t)是第t天病人在总人数中所占的比例。
λ是每个病人每天有效接触的平均人数(日接触率)。
i0是初始时刻(t=0)病人的比例。
σ是整个传染期内每个病人有效接触的平均人数(接触数)。
实验要求:
求出微分方程的解析解i(t)。
取λ=0.2,σ=3,t=0~40,画出如下所示的图形。
试编写
一个m文件来实现。
其中
蓝色实线为i(0)=0.2时的i~t曲线(第1条);
黑色虚点线为过点(0,1-1/σ)的水平线(第2条);
红色虚线为i(0)=0.9时的i~t曲线(第3条)。
[提示]
图例标注可用
legend('i(0)=0.2','1-1/¦σ','i(0)=0.9');
5.实验7-5传染病模型4(SIR模型)
(参考教材p140-141)
SIR模型的方程:
di/dt=λsi-μii(0)=i0
ds/dt=-λsis(0)=s0
实验要求:
1.设λ=1,μ=0.3,i(0)=0.02,s(0)=0.98。
输入p139的程序,并修改程序中
的[t,x],使得输出的数据格式如下(提示:
取4位小数,使用四舍五入取整函数round,
矩阵剪裁和拼接):
ans=
Columns1through6
012345
0.020.0390.07320.12850.20330.2795
0.980.95250.90190.81690.69270.5438
Columns7through12
67891015
0.33120.34440.32470.28630.24180.0787
0.39950.28390.20270.14930.11450.0543
Columns13through18
202530354045
0.02230.00610.00170.00050.00010
0.04340.04080.0401
0.03990.03990.0398
2.运行结果与教材p140的内容比较。
[提示]
1)求解微分方程的数值解函数ode45,格式如下:
[t,x]=ode45('fun',ts,x0)
fun是由一个或多个待解方程写成的函数式m文件;
ts=[t0,tf]表示此微分方程的积分限是从t0到tf,也可以是一些离散的点,形式
为ts=[t0,t1,…,tf];
x0为初值条件。
2)等待用户反应命令pause:
程序执行到该命令时暂停,直到用户按任意键后继续(处
在命令窗口有效)。
三、实验内容
1.实验7-1传染病模型2(SI模型)——画di/dt~i曲线图
在matlab中建立M文件fun1.m
代码如下:
functiony=fun(x)
k=0.1;
y=k*x*[1-x];
Fun2.m
代码如下:
functiony=fun(x)
k=0.1;
y=-k*x*[1-x];
在命令行输入以下代码:
fplot('fun1',[01.100.03]);
x=fminbnd('fun2',0,1);
y=0.1*x*(1-x);
holdon;
plot([0,x],[y,y],'-',[x,x],[0,y],'-');
text(0,y,'(di/dt)m','VerticalAlignment','bottom');
text(x,-0.001,num2str(x),'HorizontalAlignment','center');
title('SI模型di/dt~i曲线');
xlabel('i');
ylabel('di/dt');
holdoff
2.实验7-2传染病模型2(SI模型)——画i~t曲线图
在matlab中建立M文件fun22.m
代码如下:
k=0.2;
x0=0.15;
x=dsolve('Dx=k*x*(1-x)','x(0)=x0');
tt=linspace(0,31,1001);
fori=1:
1001
t=tt(i);
xx(i)=eval(x);
end
plot(tt,xx)
axis([0,31,0,1.1]);
title('图1SI模型i~t曲线');
xlabel('t(天)');
ylabel('i(病人所占比例)');
在命令行输入以下代码:
fun22;
3.实验7-3传染病模型3(SIS模型)——画di/dt~i曲线图
在matlab中建立M文件fun3.m
代码如下:
functiony=fun(x)
a=0.1;
b=1.5;
y=-a*x*[x-(1-1/b)];
在命令行输入以下代码:
fplot('fun3',[00.4-0.00050.003]);
x=fminbnd('fun3',0,1);
title('SIS模型di/dt~i曲线');
xlabel('i');
ylabel('di/dt');
>>holdon
>>plot([0,0.4],[0,0])
4.实验7-4传染病模型3(SIS模型)——画i~t曲线图
在matlab中建立M文件fun4.m
代码如下:
functiony=fun(x)
x=dsolve('Dx=-0.2*x*(x-(1-1/3))','x(0)=0.2');
tt=linspace(0,41,1001);
fori=1:
1001
t=tt(i);
xx(i)=eval(x);
end
plot(tt,xx);
holdon;
plot([0,40],[1-1/3,1-1/3],'-k');
x=dsolve('Dx=-0.2*x*(x-(1-1/3))','x(0)=0.9');
tt=linspace(0,41,1001);
fori=1:
1001
t=tt(i);
xx(i)=eval(x);
end
plot(tt,xx,'-r');
axis([0,40,0,1]);
title('图1SI模型i~t曲线(λ=0.2,σ=3)');
xlabel('t(天)');
ylabel('i(病人所占比例)');
legend('i(0)=0.2','1-1/σ','i(0)=0.9');
在命令行输入以下代码:
fun4;
5.实验7-5传染病模型4(SIR模型)
在matlab中建立M文件fun5.m
代码如下:
functiony=fun(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('fun5',ts,x0);
>>plot(t,x(:
1),t,x(:
2)),grid,pause
>>plot(x(:
2),x(:
1)),grid,
四、实验结果及其分析
1.实验7-1传染病模型2(SI模型)——画di/dt~i曲线图
分析:
当i=1/2时di/dt达到最大值(di/dt)m,这时病人增加得在最快,可以认为是医院的门诊量最大的一天,预示着传染病高潮的到来,是医疗卫生部门关注的时刻。
当t趋近于无穷时i趋近于1,即所有人终将被传染,全部变成病人,着显然不符合实际。
原因是模型中没有考虑到病人可以治愈,人群中的健康者只能变成病人,病人不会再变成健康者。
2.实验7-2传染病模型2(SI模型)——画i~t曲线图
分析:
当i=1/2时di/dt达到最大值(di/dt)m,这时病人增加得在最快,可以认为是医院的门诊量最大的一天,预示着传染病高潮的到来,是医疗卫生部门关注的时刻。
当t趋近于无穷时i趋近于1,即所有人终将被传染,全部变成病人,着显然不符合实际。
原因是模型中没有考虑到病人可以治愈,人群中的健康者只能变成病人,病人不会再变成健康者。
3.实验7-3传染病模型3(SIS模型)——画di/dt~i曲线图
分析:
是一个阈值,当
>1时,i(t)的增减性取决于i0的大小,但其极限值i(无穷)=1-1/
随
的增加而增加(试从
的含义给予解释);当
<=1时,病人的比例i(t)越来越小。
最终趋近于0,这是由于传染期内经有效接触从而使健康者变成病人数不超过原来的病人数的缘故。
4.实验7-4传染病模型3(SIS模型)——画i~t曲线图
分析:
是一个阈值,当
>1时,i(t)的增减性取决于i0的大小,但其极限值i(无穷)=1-1/
随
的增加而增加(试从
的含义给予解释);当
<=1时,病人的比例i(t)越来越小。
最终趋近于0,这是由于传染期内经有效接触从而使健康者变成病人数不超过原来的病人数的缘故。
5.实验7-5传染病模型4(SIR模型)
ans=
00.02000.9800
1.00000.03900.9525
2.00000.07320.9019
3.00000.12850.8169
4.00000.20330.6927
5.00000.27950.5438
6.00000.33120.3995
7.00000.34440.2839
8.00000.32470.2027
9.00000.28630.1493
10.00000.24180.1145
11.00000.19860.0917
12.00000.15990.0767
13.00000.12720.0665
14.00000.10040.0593
15.00000.07870.0543
16.00000.06140.0507
17.00000.04780.0480
18.00000.03710.0460
19.00000.02870.0445
20.00000.02230.0434
21.00000.01720.0426
22.00000.01330.0419
23.00000.01030.0415
24.00000.00790.0411
25.00000.00610.0408
26.00000.00470.0406
27.00000.00360.0404
28.00000.00280.0403
29.00000.00220.0402
30.00000.00170.0401
31.00000.00130.0400
32.00000.00100.0400
33.00000.00080.0400
34.00000.00060.0399
35.00000.00050.0399
36.00000.00040.0399
37.00000.00030.0399
38.00000.00020.0399
39.00000.00020.0399
40.00000.00010.0399
41.00000.00010.0399
42.00000.00010.0399
43.00000.00010.0399
44.00000.00000.0398
45.00000.00000.0398
46.00000.00000.0398
47.00000.00000.0398
48.00000.00000.0398
49.00000.00000.0398
50.00000.00000.0398
分析:
如果仅当病人比例i(t)有一段增长的时期才认为传染病在蔓延,那么1/σ是一个阈值,当s0>1/σ时传染病就会蔓延,而减小传染期接触数σ,即提高阈值1/σ,是的s0<=1/σ,传染病就不会蔓延。
σ减小时,s∞增加,im降低,也控制了蔓延的程度。
在σ=λ/μ中,人们的卫生水平提高,日接触率λ越小;医疗水平越高,日治愈率μ越大,于是σ越小,所以提高卫生水平和医疗水平有助于控制传染病的蔓延。