数学建模实验答案_微分方程模型

1、实验07 微分方程模型（2学时）（第5章 微分方程模型）1.（验证）传染病模型2（SI模型）p136138传染病模型2（SI模型）：其中，i(t)是第t天病人在总人数中所占的比例。k是每个病人每天有效接触的平均人数（日接触率）。i0是初始时刻（t=0）病人的比例。1.1 画曲线图p136138取k=0.1，画出的曲线图，求i为何值时达到最大值，并在曲线图上标注。参考程序：%传染病模型2（SI模型）的di/dti曲线图%文件名：p137fig2.m%=0.1clear; clc;fplot('0.1*x*(1-x)',0 1.1 0 0.03);x=fminbnd('-0

2、.1*x*(1-x)',0,1)y=0.1*x*(1-x)hold onplot(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/dti曲线');xlabel('i');ylabel('di/d

3、t');hold off;提示：fplot, fminbnd, plot, text, title, xlabel1）画曲线图用fplot函数，调用格式如下：fplot(fun,lims)fun必须为一个M文件的函数名或对变量x的可执行字符串。若lims取xmin xmax，则x轴被限制在此区间上。若lims取xmin xmax ymin ymax，则y轴也被限制。本题可用fplot('0.1*x*(1-x)',0 1.1 0 0.03);2）求最大值用求解边界约束条件下的非线性最小化函数fminbnd，调用格式如下：x=fminbnd('fun',x1

4、,x2)fun必须为一个M文件的函数名或对变量x的可执行字符串。返回自变量x在区间x1<x<x2上函数取最小值时的x值。本题可用x=fminbnd('-0.1*x*(1-x)',0,1)y=0.1*x*(1-x)3）指示最大值坐标用线性绘图函数plot，调用格式如下：plot(x1,y1, '颜色 线型 数据点图标', x2,y2, '颜色 线型 数据点图标',)本题可用hold on; %在上面的同一张图上画线（同坐标系）plot(0,x,y,y,':',x,x,0,y,':');4）图形的标注使用文

5、本标注函数text，调用格式如下：格式1text(x,y,文本标识内容, 'HorizontalAlignment', '字符串1')x,y给定标注文本在图中添加的位置。'HorizontalAlignment'为水平控制属性，控制文本标识起点位于点(x,y)同一水平线上。'字符串1' 为水平控制属性值，取三个值之一：'left'，点(x,y)位于文本标识的左边。'center' ，点(x,y)位于文本标识的中心点。'right' ，点(x,y)位于文本标识的右边。格式2text(x

6、,y, 文本标识内容, 'VerticalAlignment', '字符串2')x,y给定标注文本在图中添加的位置。'VerticalAlignment'为垂直控制属性，控制文本标识起点位于点(x,y)同一垂直线上。'字符串1' 为垂直控制属性值，取四个值之一：'middle'，'top'，'cap'，'baseline'，'bottom'。（对应位置可在命令窗口应用确定）本题可用text(0,y,'(di/dt)m','Ver

7、ticalAlignment','bottom');text(x,-0.001,num2str(x),'HorizontalAlignment','center');5）坐标轴标注调用函数xlabel，ylabel和title本题可用title('SI模型di/dti曲线');xlabel('i');ylabel('di/dt'); 程序运行结果（比较138图2）：（在图形窗口菜单选择Edit/Copy Figure，复制图形）1.2 画it曲线图p136138求出微分方程的解析解i(t)，

8、画出it曲线（i(0)=0.15, k=0.2, t=030）（见138图1比较）。参考程序：% 5.1 传染病模型模型2% 文件名：p136fig1.m% di/dt=ki(1-i), i(0)=i0clear; clc;x=dsolve('Dx=k*x*(1-x)','x(0)=x0') %求微分方程的解析解，为符号表达式x0=0.15; k=0.2;%xi对应i，xi0对应i0，k对应tt=0:0.1:30;%时间单位为天for s=1:length(tt)%x的表达式中没有点运算，按标量运算取值xx t=tt(s); xx(s)=eval(x);%给出x

9、i0=0.2,k=0.2,t,求符号表达式xi的对应值end %xx为复数表示plot(tt,xx);axis(0 31 0 1.1);title('图1 SI模型的it曲线');xlabel('t （天）'); ylabel('i （病人所占比例）');提示：1) 求解微分方程dsolve，见提示；2) 画出it曲线（i(0)=0.15, =0.2, t=030）用for循环，函数length, eval, plot, axis, title, xlabel, ylabel。 程序运行结果（见138图1）：命令窗口中的结果：图形窗口中的结果（比

10、较138图1）：2.（编程）传染病模型3（SIS模型）已知传染病模型3（SIS模型）：其中，i(t)是第t天病人在总人数中所占的比例。是每个病人每天有效接触的平均人数（日接触率）。i0是初始时刻（t=0）病人的比例。是整个传染期内每个病人有效接触的平均人数（接触数）。2.1 画曲线图p138139取=0.1，=1.5，画出如下所示的曲线图。试编写一个m文件来实现。（在图形窗口菜单选择Edit/Copy Figure，复制图形）（注：p139图3）提示：用fplot函数画出的曲线图；在上图上用plot函数画一条过原点的水平线；用title, xlabel, ylabel标注。 编写的M文件和运行

11、结果（见139图3）：%传染病模型3（SIS模型）的di/dti曲线图%文件名：p138fig3.m%=0.1, =1.5clear; clc;fplot('-0.1*x*(x-(1-1/1.5)',0 0.4 -0.0005 0.003);hold on;plot(0,0.4,0,0,':');title('SIS模型的di/dti曲线');xlabel('i');ylabel('di/dt');2.2 画it曲线图p138139要求：求出微分方程的解析解i(t)。取=0.2, =3, t=040，画出如下所示的

12、图形。试编写一个m文件来实现。（注：p139图4）其中蓝色实线为i(0)=0.2时的it曲线（第1条）；黑色虚点线为过点（0, 1-1/）的水平线（第2条）；红色虚线为i(0)=0.9时的it曲线（第3条）。提示图例标注可用legend('i(0)=0.2','1-1/¦','i(0)=0.9'); 编写的M文件和运行结果（比较139图4）：解法一：程序：%传染病模型3（SIS模型）的it曲线图%文件名：p138fig4.mclear; clc;%=0.2, =3，i(0)=0.2, x代表ix=dsolve('Dx=-0.2*

13、x*(x-(1-1/3)','x(0)=0.2')%求微分方程的解析解，为符号表达式tt=0:0.1:40;%时间单位为天for i=1:length(tt) t=tt(i); xx(i)=eval(x);endplot(tt,xx);hold on;plot(0,41,1-1/3,1-1/3,'-.k');%=0.2, =3，i(0)=0.9x=dsolve('Dx=-0.2*x*(x-(1-1/3)','x(0)=0.9') tt=0:0.1:40;%时间单位为天for i=1:length(tt) t=tt(i);

14、xx(i)=eval(x);endplot(tt,xx,':r');legend('i(0)=0.2','1-1/','i(0)=0.9');axis(0 40 0 1);title('图1 SI模型的it曲线(=0.2,=3)');xlabel('t （天）'); ylabel('i （病人所占比例）');命令窗口的结果：图形窗口的结果：解法二：程序：%传染病模型3（SIS模型）的it曲线图%文件名：p138fig4.mclear; clc;%=0.2, =3，x代表ix=dsol

15、ve('Dx=-0.2*x*(x-(1-1/3)','x(0)=x0')%求微分方程的解析解，为符号表达式tt=0:0.1:40;%时间单位为天for x0=0.2,0.9%i(0)=0.2,0.9 for t=tt xx(2-(x0=0.2),round(t/0.1)+1)=eval(x); endendplot(tt,xx(1,:),'-b',0,41,1-1/3,1-1/3,'-.k',tt,xx(2,:),':r');legend('i(0)=0.2','1-1/','

16、;i(0)=0.9');axis(0 40 0 1);title('图1 SI模型的it曲线(=0.2,=3)');xlabel('t （天）'); ylabel('i （病人所占比例）');命令窗口的结果：图形窗口的结果：与解法一相同解法三：程序%传染病模型3（SIS模型）的it曲线图%文件名：p138fig4.mclear; clc;x=dsolve('Dx=-lam*x*(x-(1-1/si)','x(0)=x0')%求微分方程的解析解，为符号表达式tt=0:0.1:40;%时间单位为天lam=0.2

17、; si=3; %=0.2, =3，x代表ifor x0=0.2,0.9 %i(0)=0.2,0.9 for t=tt xx(2-(x0=0.2),round(t/0.1)+1)=eval(x); endendplot(tt,xx(1,:),'-b',0,41,1-1/3,1-1/3,'-.k',tt,xx(2,:),':r');legend('i(0)=0.2','1-1/','i(0)=0.9');axis(0 40 0 1);title('图1 SI模型的it曲线(=0.2,=3)&#

18、39;);xlabel('t （天）'); ylabel('i （病人所占比例）');命令窗口的结果：图形窗口的结果：与解法一相同3.（验证）传染病模型4（SIR模型）p140141SIR模型的方程：设 =1， =0.3，i(0)=0.02，s(0)=0.98。输入p140的程序并运行，结果与教材p141的图7和图8比较。ode45, pause的用法见提示。 2个M文件（见140）和运行结果（比较141图7、图8）：函数M文件：%5.1 传染病模型模型4（SIR模型）%文件名：ill.mfunction y=ill(t,x)a=1; b=0.3; % 用a表示

19、， 用b表示y=a*x(1)*x(2)-b*x(1), -a*x(1)*x(2)' %i用x(1)表示，s用x(2)表示命令M文件：%5.1 传染病模型模型4（SIR模型）%文件名：p140.mclear; clc;ts=0:50;x0=0.02, 0.98;t,x=ode45('ill',ts,x0); t,xplot(t,x(:,1),t,x(:,2), grid, pause%图7 i(t), s(t)图形plot(x(:,2),x(:,1); grid,%图8 is图形（相轨线）i(t), s(t)图形（比较141图7）： is图形（相轨线）（比较141图8）：

20、4.（验证）人口指数增长模型参数估计及结果分析（美国1790-2000年人口）p163164美国1790-2000年人口统计数据（以百万为单位）年17901800181018201830184018501860187018801890人口3.95.37.29.612.917.123.231.438.650.262.9年19001910192019301940195019601970198019902000人口76.092.0106.5123.2131.7150.7179.3204.0226.5251.4281.4人口指数增长模型：x(t) = x0 e r t(1) 用表中数据进行数据拟合求参

21、数r，x0。将x(t) = x0 e r t两边取对数，可得y = rt + a其中， y = ln x, a = ln x0，即x0 = exp(a)。采用线性最小二乘法进行数据拟合，用MATLAB中的函数polyfit计算。（r为每10年估计的人口增长率，x0为1790年估计的初始人口数）以下是M文件：clc; format compact;x=3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 . 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7 . 150.7 179.3 204.0 226.5 251.4 281.4;t=0:l

22、ength(x)-1; %t=0为1790年，t=1为1800年，y=log(x); %取1790-2000年的数据ra=polyfit(t,y,1);disp('用1790-2000年数据估计的参数为：') r=ra(1)x0=exp(ra(2)(1) 运行程序并给出结果（见167）：(2) 人口指数增长模型计算结果与实际数据比较（数据表）以下是M文件：clc; format compactx=3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 . 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7 . 150.7 179

23、.3 204.0 226.5 251.4 281.4;t=0:length(x)-1;y=log(x); %取1790-2000年的数据ra=polyfit(t,y,1);r=ra(1); x0=exp(ra(2);x2=x0*exp(r*t);disp('指数增长模型拟合美国人口数据的结果(x2)')format short g;disp(1790+10*t; x; round(10*x2)/10');(2) 运行程序并给出结果（见167表4中x2列）：(3) 人口指数增长模型计算结果与实际数据比较（拟合图形）以下是M文件：x=3.9 5.3 7.2 9.6 12.9

24、 17.1 23.2 31.4 . 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7 . 150.7 179.3 204.0 226.5 251.4 281.4;t=0:length(x)-1;y=log(x); %取1790-2000年的数据ra=polyfit(t,y,1);r=ra(1); x0=exp(ra(2);t2=linspace(0,length(x)-1,30);x2=x0*exp(r*t2);plot(t,x,'r+',t2,x2,'b');title('指数增长模型拟合图形');(3) 运

25、行程序并给出结果（见168图3(b)）：5.（验证，编程）估计阻滞增长模型的参数和绘制图形p165168美国1790-2000年人口统计数据（以百万为单位）年17901800181018201830184018501860187018801890人口3.95.37.29.612.917.123.231.438.650.262.9年19001910192019301940195019601970198019902000人口76.092.0106.5123.2131.7150.7179.3204.0226.5251.4281.4人口阻滞增长模型：或(1)（验证）用1860-1990年的数据拟合估计

26、参数r, xm用下面方程估计参数r, xm。程序如下：%用数值微分的三点公式计算美国人口增长率(/10年)clear; clc;x=3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 . 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7 . 150.7 179.3 204.0 226.5 251.4;%1790-1990年的人口xx=x(1,8:end);%取1860-1990年的数据len=length(xx); h=1; dx=-3*xx(1)+4*xx(2)-xx(3), xx(3:len)-xx(1:len-2), . %求数值微

27、分dx 3*xx(len)-4*xx(len-1)+xx(len-2)/(2*h);y=dx./xx; sr=polyfit(xx,y,1);r=sr(2)xm=-r/sr(1)(1) 运行程序并给出结果（见168差异大，书上结果数据处理过）：(2)（编程）计算阻滞增长模型拟合美国人口数据（1790-1990）取参数r=0.2557, xm=392.0886, x0=3.9，编写程序计算1790-1990年的人口，参考的输出如下（第1列为年，第2列为实际人口，第3列为计算人口）：(2) 给出程序及其运行结果（见167表4中x列）：clear; clc;format compact;x=3.9

28、5.3 7.2 9.6 12.9 17.1 23.2 31.4 . 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7 . 150.7 179.3 204.0 226.5 251.4; %1790-1990年的人口t=0:length(x)-1;x0=3.9; r=0.2557; xm=392.0886; %代入参数值xx=xm./(1+(xm/x0-1)*exp(-r*t); %计算人口format short g;1790+10*t;x;round(10*xx)/10'(3)（编程）绘制阻滞增长模型拟合图形（1790-1990）取参数r=0.2557, xm=392.0886, x0=3.9，编写程序绘制1790-1990年的人口数据拟合图形。参考图形如下，实线为模型曲线，+为实际数据值。(3) 给出程序及其运行结果（见168图4）：clear; clc; format compact;x=3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 . 38.6 50.2 62.9 76.0 92.0 106.5