實(shí)驗(yàn)07講評、參考答案_微分方程模型(2學(xué)時(shí))

1、實(shí)驗(yàn)07講評、參考答案講 評未按時(shí)交的同學(xué)數(shù)學(xué)：01邊清水，14黃浦，28陸杭濤，34譚世韜，50鐘鑫信科：19施磊批改情況：附參考答案：數(shù)學(xué)建模實(shí)驗(yàn)王平實(shí)驗(yàn)07 微分方程模型（2學(xué)時(shí)）（第5章 微分方程模型）1.（驗(yàn)證）傳染病模型2（SI模型）p136138傳染病模型2（SI模型）：其中，i(t)是第t天病人在總?cè)藬?shù)中所占的比例。k是每個病人每天有效接觸的平均人數(shù)（日接觸率）。i0是初始時(shí)刻（t=0）病人的比例。1.1 畫曲線圖p136138取k=0.1，畫出的曲線圖，求i為何值時(shí)達(dá)到最大值，并在曲線圖上標(biāo)注。參考程序：%傳染病模型2（SI模型）的di/dti曲線圖%文件名：p137fig2

2、.m%=0.1clear; clc;fplot('0.1*x*(1-x)',0 1.1 0 0.03);x=fminbnd('-0.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','cen

3、ter');title('SI模型的di/dti曲線');xlabel('i');ylabel('di/dt');hold off;提示：fplot, fminbnd, plot, text, title, xlabel1）畫曲線圖用fplot函數(shù)，調(diào)用格式如下：fplot(fun,lims)fun必須為一個M文件的函數(shù)名或?qū)ψ兞縳的可執(zhí)行字符串。若lims取xmin xmax，則x軸被限制在此區(qū)間上。若lims取xmin xmax ymin ymax，則y軸也被限制。本題可用fplot('0.1*x*(1-x)',0 1

4、.1 0 0.03);2）求最大值用求解邊界約束條件下的非線性最小化函數(shù)fminbnd，調(diào)用格式如下：x=fminbnd('fun',x1,x2)fun必須為一個M文件的函數(shù)名或?qū)ψ兞縳的可執(zhí)行字符串。返回自變量x在區(qū)間x1<x<x2上函數(shù)取最小值時(shí)的x值。本題可用x=fminbnd('-0.1*x*(1-x)',0,1)y=0.1*x*(1-x)3）指示最大值坐標(biāo)用線性繪圖函數(shù)plot，調(diào)用格式如下：plot(x1,y1, '顏色 線型 數(shù)據(jù)點(diǎn)圖標(biāo)', x2,y2, '顏色 線型 數(shù)據(jù)點(diǎn)圖標(biāo)',)本題可用hold o

5、n; %在上面的同一張圖上畫線（同坐標(biāo)系）plot(0,x,y,y,':',x,x,0,y,':');4）圖形的標(biāo)注使用文本標(biāo)注函數(shù)text，調(diào)用格式如下：格式1text(x,y,文本標(biāo)識內(nèi)容, 'HorizontalAlignment', '字符串1')x,y給定標(biāo)注文本在圖中添加的位置。'HorizontalAlignment'為水平控制屬性，控制文本標(biāo)識起點(diǎn)位于點(diǎn)(x,y)同一水平線上。'字符串1' 為水平控制屬性值，取三個值之一：'left'，點(diǎn)(x,y)位于文本標(biāo)識的左邊。

6、'center' ，點(diǎn)(x,y)位于文本標(biāo)識的中心點(diǎn)。'right' ，點(diǎn)(x,y)位于文本標(biāo)識的右邊。格式2text(x,y, 文本標(biāo)識內(nèi)容, 'VerticalAlignment', '字符串2')x,y給定標(biāo)注文本在圖中添加的位置。'VerticalAlignment'為垂直控制屬性，控制文本標(biāo)識起點(diǎn)位于點(diǎn)(x,y)同一垂直線上。'字符串1' 為垂直控制屬性值，取四個值之一：'middle'，'top'，'cap'，'baseline&#

7、39;，'bottom'。（對應(yīng)位置可在命令窗口應(yīng)用確定）本題可用text(0,y,'(di/dt)m','VerticalAlignment','bottom');text(x,-0.001,num2str(x),'HorizontalAlignment','center');5）坐標(biāo)軸標(biāo)注調(diào)用函數(shù)xlabel，ylabel和title本題可用title('SI模型di/dti曲線');xlabel('i');ylabel('di/dt'); 程序運(yùn)

8、行結(jié)果（比較138圖2）：（在圖形窗口菜單選擇Edit/Copy Figure，復(fù)制圖形）1.2 畫it曲線圖p136138求出微分方程的解析解i(t)，畫出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') %求微分方程的解析解，為符號表達(dá)式x0=0.15; k=0.2;%xi對應(yīng)i，xi0對應(yīng)i0，k對應(yīng)tt=0:0.1:3

9、0;%時(shí)間單位為天for s=1:length(tt)%x的表達(dá)式中沒有點(diǎn)運(yùn)算，按標(biāo)量運(yùn)算取值xx t=tt(s); xx(s)=eval(x);%給出xi0=0.2,k=0.2,t,求符號表達(dá)式xi的對應(yīng)值end %xx為復(fù)數(shù)表示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循環(huán)，函數(shù)length

10、, eval, plot, axis, title, xlabel, ylabel。 程序運(yùn)行結(jié)果（見138圖1）：命令窗口中的結(jié)果：圖形窗口中的結(jié)果（比較138圖1）：2.（編程）傳染病模型3（SIS模型）已知傳染病模型3（SIS模型）：其中，i(t)是第t天病人在總?cè)藬?shù)中所占的比例。是每個病人每天有效接觸的平均人數(shù)（日接觸率）。i0是初始時(shí)刻（t=0）病人的比例。是整個傳染期內(nèi)每個病人有效接觸的平均人數(shù)（接觸數(shù)）。2.1 畫曲線圖p138139取=0.1，=1.5，畫出如下所示的曲線圖。試編寫一個m文件來實(shí)現(xiàn)。（在圖形窗口菜單選擇Edit/Copy Figure，復(fù)制圖形）（注：p139圖

11、3）提示：用fplot函數(shù)畫出的曲線圖；在上圖上用plot函數(shù)畫一條過原點(diǎn)的水平線；用title, xlabel, ylabel標(biāo)注。 編寫的M文件和運(yùn)行結(jié)果（見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

12、('di/dt');2.2 畫it曲線圖p138139要求：求出微分方程的解析解i(t)。取=0.2, =3, t=040，畫出如下所示的圖形。試編寫一個m文件來實(shí)現(xiàn)。（注：p139圖4）其中藍(lán)色實(shí)線為i(0)=0.2時(shí)的it曲線（第1條）；黑色虛點(diǎn)線為過點(diǎn)（0, 1-1/）的水平線（第2條）；紅色虛線為i(0)=0.9時(shí)的it曲線（第3條）。提示圖例標(biāo)注可用legend('i(0)=0.2','1-1/¦','i(0)=0.9'); 編寫的M文件和運(yùn)行結(jié)果（比較139圖4）：解法一：程序：%傳染病模型3（SIS模型）的

13、it曲線圖%文件名：p138fig4.mclear; clc;%=0.2, =3，i(0)=0.2, x代表ix=dsolve('Dx=-0.2*x*(x-(1-1/3)','x(0)=0.2')%求微分方程的解析解，為符號表達(dá)式tt=0:0.1:40;%時(shí)間單位為天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-

14、1/3)','x(0)=0.9') tt=0:0.1:40;%時(shí)間單位為天for i=1:length(tt) t=tt(i); 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 （病人所占比例）');命令窗口的結(jié)果：圖形窗口

15、的結(jié)果：解法二：程序：%傳染病模型3（SIS模型）的it曲線圖%文件名：p138fig4.mclear; clc;%=0.2, =3，x代表ix=dsolve('Dx=-0.2*x*(x-(1-1/3)','x(0)=x0')%求微分方程的解析解，為符號表達(dá)式tt=0:0.1:40;%時(shí)間單位為天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,'-.

16、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)');xlabel('t （天）'); ylabel('i （病人所占比例）');命令窗口的結(jié)果：圖形窗口的結(jié)果：與解法一相同解法三：程序%傳染病模型3（SIS模型）的it曲線圖%文件名：p138fig4.mclear; clc;x=dsolve('Dx=-lam*x*(

17、x-(1-1/si)','x(0)=x0')%求微分方程的解析解，為符號表達(dá)式tt=0:0.1:40;%時(shí)間單位為天lam=0.2; 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-

18、1/','i(0)=0.9');axis(0 40 0 1);title('圖1 SI模型的it曲線(=0.2,=3)');xlabel('t （天）'); ylabel('i （病人所占比例）');命令窗口的結(jié)果：圖形窗口的結(jié)果：與解法一相同3.（驗(yàn)證）傳染病模型4（SIR模型）p140141SIR模型的方程：設(shè) =1， =0.3，i(0)=0.02，s(0)=0.98。輸入p140的程序并運(yùn)行，結(jié)果與教材p141的圖7和圖8比較。ode45, pause的用法見提示。 2個M文件（見140）和運(yùn)行結(jié)果（比較141圖7、

19、圖8）：函數(shù)M文件：%5.1 傳染病模型模型4（SIR模型）%文件名：ill.mfunction y=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.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)圖形plo

20、t(x(:,2),x(:,1); grid,%圖8 is圖形（相軌線）i(t), s(t)圖形（比較141圖7）： is圖形（相軌線）（比較141圖8）：4.（驗(yàn)證）人口指數(shù)增長模型參數(shù)估計(jì)及結(jié)果分析（美國1790-2000年人口）p163164美國1790-2000年人口統(tǒng)計(jì)數(shù)據(jù)（以百萬為單位）年17901800181018201830184018501860187018801890人口3.95.37.29.612.917.123.231.438.650.262.9年19001910192019301940195019601970198019902000人口76.092.0106.5123.

21、2131.7150.7179.3204.0226.5251.4281.4人口指數(shù)增長模型：x(t) = x0 e r t(1) 用表中數(shù)據(jù)進(jìn)行數(shù)據(jù)擬合求參數(shù)r，x0。將x(t) = x0 e r t兩邊取對數(shù)，可得y = rt + a其中， y = ln x, a = ln x0，即x0 = exp(a)。采用線性最小二乘法進(jìn)行數(shù)據(jù)擬合，用MATLAB中的函數(shù)polyfit計(jì)算。（r為每10年估計(jì)的人口增長率，x0為1790年估計(jì)的初始人口數(shù)）以下是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

22、 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; %t=0為1790年，t=1為1800年，y=log(x); %取1790-2000年的數(shù)據(jù)ra=polyfit(t,y,1);disp('用1790-2000年數(shù)據(jù)估計(jì)的參數(shù)為：') r=ra(1)x0=exp(ra(2)(1) 運(yùn)行程序并給出結(jié)果（見167）：(2) 人口指數(shù)增長模型計(jì)算結(jié)果與實(shí)際數(shù)據(jù)比較（數(shù)據(jù)表）以下是M文件：clc; format compactx=3.9 5.3 7.2 9.6

23、 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:length(x)-1;y=log(x); %取1790-2000年的數(shù)據(jù)ra=polyfit(t,y,1);r=ra(1); x0=exp(ra(2);x2=x0*exp(r*t);disp('指數(shù)增長模型擬合美國人口數(shù)據(jù)的結(jié)果(x2)')format short g;disp(1790+10*t; x; round(10*x2)/10');(2) 運(yùn)行程

24、序并給出結(jié)果（見167表4中x2列）：(3) 人口指數(shù)增長模型計(jì)算結(jié)果與實(shí)際數(shù)據(jù)比較（擬合圖形）以下是M文件：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:length(x)-1;y=log(x); %取1790-2000年的數(shù)據(jù)ra=polyfit(t,y,1);r=ra(1); x0=exp(ra(2);t2=linspace(0,length(x)-1,30);x2=x0*exp(r*t

25、2);plot(t,x,'r+',t2,x2,'b');title('指數(shù)增長模型擬合圖形');(3) 運(yùn)行程序并給出結(jié)果（見168圖3(b)）：5.（驗(yàn)證，編程）估計(jì)阻滯增長模型的參數(shù)和繪制圖形p165168美國1790-2000年人口統(tǒng)計(jì)數(shù)據(jù)（以百萬為單位）年17901800181018201830184018501860187018801890人口3.95.37.29.612.917.123.231.438.650.262.9年19001910192019301940195019601970198019902000人口76.092.0106

26、.5123.2131.7150.7179.3204.0226.5251.4281.4人口阻滯增長模型：或(1)（驗(yàn)證）用1860-1990年的數(shù)據(jù)擬合估計(jì)參數(shù)r, xm用下面方程估計(jì)參數(shù)r, xm。程序如下：%用數(shù)值微分的三點(diǎn)公式計(jì)算美國人口增長率(/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年的數(shù)據(jù)l

27、en=length(xx); h=1; dx=-3*xx(1)+4*xx(2)-xx(3), xx(3:len)-xx(1:len-2), . %求數(shù)值微分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) 運(yùn)行程序并給出結(jié)果（見168差異大，書上結(jié)果數(shù)據(jù)處理過）：(2)（編程）計(jì)算阻滯增長模型擬合美國人口數(shù)據(jù)（1790-1990）取參數(shù)r=0.2557, xm=392.0886, x0=3.9，編寫程序計(jì)算1790-1990年的人口，參考的輸出如下（第1列為年，

28、第2列為實(shí)際人口，第3列為計(jì)算人口）：(2) 給出程序及其運(yùn)行結(jié)果（見167表4中x列）：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 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; %代入?yún)?shù)值xx=xm./(1+(xm/x0-1)*exp(-r*t); %計(jì)算人口format short g;1790+10*t;x;round(10*xx)/10'(3)（編程）繪制阻滯增長模型擬合圖形（1790-1990）取參數(shù)r=0.2557, xm=392.0886, x0=3.9，編寫程序繪制1790-1990年的人口數(shù)據(jù)擬合圖形。參考圖形如下，實(shí)線為模型曲線，+為實(shí)際數(shù)據(jù)值。(3) 給出程序及其運(yùn)行結(jié)果（見168圖4）：clear; clc; format compact;x=3.9 5.3 7.2 9.6 12.9 17.1