实验07讲评参考答案微分方程模型2学时.docx
- 文档编号:9037165
- 上传时间:2023-02-02
- 格式:DOCX
- 页数:48
- 大小:2.16MB
实验07讲评参考答案微分方程模型2学时.docx
《实验07讲评参考答案微分方程模型2学时.docx》由会员分享,可在线阅读,更多相关《实验07讲评参考答案微分方程模型2学时.docx(48页珍藏版)》请在冰豆网上搜索。
实验07讲评参考答案微分方程模型2学时
实验07讲评、参考答案
讲评
未按时交的同学
数学:
01边清水,14黄浦,28陆杭涛,34谭世韬,50钟鑫
信科:
19施磊
批改情况:
附参考答案:
实验07微分方程模型(2学时)
(第5章微分方程模型)
1.(验证)传染病模型2(SI模型)p136~138
传染病模型2(SI模型):
其中,
i(t)是第t天病人在总人数中所占的比例。
k是每个病人每天有效接触的平均人数(日接触率)。
i0是初始时刻(t=0)病人的比例。
1.1画
曲线图p136~138
取k=0.1,画出
的曲线图,求i为何值时
达到最大值,并在曲线图上标注。
参考程序:
%传染病模型2(SI模型)的di/dt~i曲线图
%文件名:
p137fig2.m
%λ=0.1
clear;clc;
fplot('0.1*x*(1-x)',[01.100.03]);
x=fminbnd('-0.1*x*(1-x)',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;
提示:
fplot,fminbnd,plot,text,title,xlabel
1)画曲线图
用fplot函数,调用格式如下:
fplot(fun,lims)
fun必须为一个M文件的函数名或对变量x的可执行字符串。
若lims取[xminxmax],则x轴被限制在此区间上。
若lims取[xminxmaxyminymax],则y轴也被限制。
本题可用
fplot('0.1*x*(1-x)',[01.100.03]);
2)求最大值
用求解边界约束条件下的非线性最小化函数fminbnd,调用格式如下:
x=fminbnd('fun',x1,x2)
fun必须为一个M文件的函数名或对变量x的可执行字符串。
返回自变量x在区间x1 本题可用 x=fminbnd('-0.1*x*(1-x)',0,1) y=0.1*x*(1-x) 3)指示最大值坐标 用线性绘图函数plot,调用格式如下: plot(x1,y1,'颜色线型数据点图标',x2,y2,'颜色线型数据点图标',…) 本题可用 holdon;%在上面的同一张图上画线(同坐标系) plot([0,x],[y,y],': ',[x,x],[0,y],': '); 4)图形的标注 使用文本标注函数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'); 5)坐标轴标注 调用函数xlabel,ylabel和title 本题可用 title('SI模型di/dt~i曲线'); xlabel('i'); ylabel('di/dt'); ☆程序运行结果(比较[138]图2): (在图形窗口菜单选择Edit/CopyFigure,复制图形) 1.2画i~t曲线图p136~138 求出微分方程的解析解i(t),画出i~t曲线(i(0)=0.15,k=0.2,t=0~30)(见[138]图1比较)。 参考程序: %5.1传染病模型——模型2 %文件名: p136fig1.m %di/dt=ki(1-i),i(0)=i0 clear;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;%时间单位为天 fors=1: length(tt)%x的表达式中没有点运算,按标量运算取值xx t=tt(s); xx(s)=eval(x);%给出xi0=0.2,k=0.2,t,求符号表达式xi的对应值 end%xx为复数表示 plot(tt,xx); axis([03101.1]); title('图1SI模型的i~t曲线'); xlabel('t(天)');ylabel('i(病人所占比例)'); 提示: 1)求解微分方程dsolve],见提示; 2)画出i~t曲线(i(0)=0.15,λ=0.2,t=0~30) 用for循环,函数length,eval,plot,axis,title,xlabel,ylabel。 ☆程序运行结果(见[138]图1): 命令窗口中的结果: 图形窗口中的结果(比较[138]图1): 2.(编程)传染病模型3(SIS模型) 已知传染病模型3(SIS模型): 其中, i(t)是第t天病人在总人数中所占的比例。 λ是每个病人每天有效接触的平均人数(日接触率)。 i0是初始时刻(t=0)病人的比例。 σ是整个传染期内每个病人有效接触的平均人数(接触数)。 2.1画 曲线图p138~139 取λ=0.1,σ=1.5,画出如下所示的 曲线图。 试编写一个m文件来实现。 (在图形窗口菜单选择Edit/CopyFigure,复制图形) (注: p139图3) 提示: 用fplot函数画出 的曲线图; 在上图上用plot函数画一条过原点的水平线; 用title,xlabel,ylabel标注。 ★编写的M文件和运行结果(见[139]图3): %传染病模型3(SIS模型)的di/dt~i曲线图 %文件名: p138fig3.m %λ=0.1,σ=1.5 clear;clc; fplot('-0.1*x*(x-(1-1/1.5))',[00.4-0.00050.003]); holdon; plot([0,0.4],[0,0],': '); title('SIS模型的di/dt~i曲线'); xlabel('i'); ylabel('di/dt'); 2.2画i~t曲线图p138~139 要求: 求出微分方程的解析解i(t)。 取λ=0.2,σ=3,t=0~40,画出如下所示的图形。 试编写一个m文件来实现。 (注: p139图4) 其中 蓝色实线为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'); ★编写的M文件和运行结果(比较[139]图4): 解法一: 程序: %传染病模型3(SIS模型)的i~t曲线图 %文件名: p138fig4.m clear;clc; %λ=0.2,σ=3,i(0)=0.2,x代表i x=dsolve('Dx=-0.2*x*(x-(1-1/3))','x(0)=0.2')%求微分方程的解析解,为符号表达式 tt=0: 0.1: 40;%时间单位为天 fori=1: length(tt) t=tt(i); xx(i)=eval(x); end plot(tt,xx); holdon; plot([0,41],[1-1/3,1-1/3],'-.k'); %λ=0.2,σ=3,i(0)=0.9 x=dsolve('Dx=-0.2*x*(x-(1-1/3))','x(0)=0.9') tt=0: 0.1: 40;%时间单位为天 fori=1: length(tt) t=tt(i); xx(i)=eval(x); end plot(tt,xx,': r'); legend('i(0)=0.2','1-1/σ','i(0)=0.9'); axis([04001]); title('图1SI模型的i~t曲线(λ=0.2,σ=3)'); xlabel('t(天)');ylabel('i(病人所占比例)'); 命令窗口的结果: 图形窗口的结果: 解法二: 程序: %传染病模型3(SIS模型)的i~t曲线图 %文件名: p138fig4.m clear;clc; %λ=0.2,σ=3,x代表i x=dsolve('Dx=-0.2*x*(x-(1-1/3))','x(0)=x0')%求微分方程的解析解,为符号表达式 tt=0: 0.1: 40;%时间单位为天 forx0=[0.2,0.9]%i(0)=0.2,0.9 fort=tt xx(2-(x0==0.2),round(t/0.1)+1)=eval(x); end end plot(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([04001]); title('图1SI模型的i~t曲线(λ=0.2,σ=3)'); xlabel('t(天)');ylabel('i(病人所占比例)'); 命令窗口的结果: 图形窗口的结果: 与解法一相同 解法三: 程序 %传染病模型3(SIS模型)的i~t曲线图 %文件名: p138fig4.m clear;clc; x=dsolve('Dx=-lam*x*(x-(1-1/si))','x(0)=x0')%求微分方程的解析解,为符号表达式 tt=0: 0.1: 40;%时间单位为天 lam=0.2;si=3;%λ=0.2,σ=3,x代表i forx0=[0.2,0.9]%i(0)=0.2,0.9 fort=tt xx(2-(x0==0.2),round(t/0.1)+1)=eval(x); end end plot(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([04001]); title('图1SI模型的i~t曲线(λ=0.2,σ=3)'); xlabel('t(天)');ylabel('i(病人所占比例)'); 命令窗口的结果: 图形窗口的结果: 与解法一相同 3.(验证)传染病模型4(SIR模型)p140~141 SIR模型的方程: 设λ=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.m functiony=ill(t,x) a=1;b=0.3;%λ用a表示,μ用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.m clear;clc; ts=0: 50; x0=[0.02,0.98]; [t,x]=ode45('ill',ts,x0);[t,x] plot(t,x(: 1),t,x(: 2)),grid,pause%图7i(t),s(t)图形 plot(x(: 2),x(: 1));grid,%图8i~s图形(相轨线) i(t),s(t)图形(比较[141]图7): i~s图形(相轨线)(比较[141]图8): 4.(验证)人口指数增长模型参数估计及结果分析(美国1790-2000年人口)p163~164 美国1790-2000年人口统计数据(以百万为单位) 年 1790 1800 1810 1820 1830 1840 1850 1860 1870 1880 1890 人口 3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 年 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 人口 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 281.4 人口指数增长模型: x(t)=x0ert (1)用表中数据进行数据拟合求参数r,x0。 将x(t)=x0ert两边取对数,可得 y=rt+a 其中,y=lnx,a=lnx0,即x0=exp(a)。 采用线性最小二乘法进行数据拟合,用MATLAB中的函数polyfit计算。 (r为每10年估计的人口增长率,x0为1790年估计的初始人口数) 以下是M文件: clc;formatcompact; x=[3.95.37.29.612.917.123.231.4... 38.650.262.976.092.0106.5123.2131.7... 150.7179.3204.0226.5251.4281.4]; t=0: length(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;formatcompact x=[3.95.37.29.612.917.123.231.4... 38.650.262.976.092.0106.5123.2131.7... 150.7179.3204.0226.5251.4281.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)') formatshortg; disp([1790+10*t;x;round(10*x2)/10]'); ☆ (2)运行程序并给出结果(见[167]表4中x2列): (3)人口指数增长模型计算结果与实际数据比较(拟合图形) 以下是M文件: x=[3.95.37.29.612.917.123.231.4... 38.650.262.976.092.0106.5123.2131.7... 150.7179.3204.0226.5251.4281.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)运行程序并给出结果(见[168]图3(b)): 5.(验证,编程)估计阻滞增长模型的参数和绘制图形p165~168 美国1790-2000年人口统计数据(以百万为单位) 年 1790 1800 1810 1820 1830 1840 1850 1860 1870 1880 1890 人口 3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 年 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 人口 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 281.4 人口阻滞增长模型: 或 (1)(验证)用1860-1990年的数据拟合估计参数r,xm 用下面方程 估计参数r,xm。 程序如下: %用数值微分的三点公式计算美国人口增长率(/10年) clear;clc; x=[3.95.37.29.612.917.123.231.4... 38.650.262.976.092.0106.5123.2131.7... 150.7179.3204.0226.5251.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),...%求数值微分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;formatcompact; x=[3.95.37.29.612.917.123.231.4... 38.650.262.976.092.0106.5123.2131.7... 150.7179.3204.0226.5251.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));%计算人口 formatshortg; [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;formatcompact; x=[3.95.37.29.612.917.123.231.4... 38.650.262.976.092.0106.5123.2131.7... 150.7179.3204.0226.5251.4];%1790-1990年的人口 t=0: length(x)-1; x0=3.9;r=0.2557;xm=392.0886; t3=linspace(0,length(x)-1,30); x3=xm./(1+(xm/x0-1)*exp(-r*t3)); plot(t,x,'r+',t3,x3,'b-'); 附1: 实验提示 第1.2题 提示: 求解微分方程dsolve 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题 ode45,pause的用法 1)求解微分方程的数值解函数ode45,格式如下: [t,x]=ode45('fun',ts,x0) fun是由一个或多个待解方程写成的函数式m文件; ts=[t0,tf]表示此微分方程的积分限是从t0到tf,也可以是一些离散的点,形式为ts=[t0,t1,…,tf]; x0为初值条件。 2)等待用户反应命令pause: 程序执行到该命令时暂停,直到用户按任意键后继续(处在命令窗口有效)。 附2: 第5章微分方程模型 [136]5.1传染病模型 [138]题1.1、1.2答案 [139]题2.1
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 实验 07 讲评 参考答案 微分方程 模型 学时
![提示](https://static.bdocx.com/images/bang_tan.gif)