1、實(shí)驗(yàn)07講評(píng)、參考答案未按時(shí)交的同學(xué)數(shù)學(xué)：01邊清水，14黃浦，28陸杭濤，34譚世韜，50鐘鑫 信科：19施磊批改情況:附參考答案:海南大學(xué)信息科學(xué)技術(shù)學(xué)院實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)課程:數(shù)學(xué)建模實(shí)驗(yàn)成績(jī)指導(dǎo)教師王平學(xué)號(hào):姓名:班級(jí):同組成員:完成日期：20_年月 日實(shí)驗(yàn)07微分方程模型（2學(xué)時(shí)）（第5章微分方程模型）1.（驗(yàn)證）傳染病模型2 （SI模型）pl36138傳染病模型2 （SI模型）：I J =i（O） = ioat其中，是'第f天病人在總?cè)藬?shù)中所占的比例。*是每個(gè)病人每天有效接觸的平均人數(shù)（日接觸率）。 巾是初始時(shí)刻（eo）病人的比例。di .1.1畫不7曲線圖p136*138取Jt

2、=0.1.畫出弓的曲線圖，求r為何值時(shí)專達(dá)到最大值，并在曲線圖atat上標(biāo)注。參考程序: %傳染病模型2 （SI模型）的di/dti曲線圖%文件名 J P 137ng2.m%k=O.ldear; cic;fplot( 0.1*x*(l-x)0 1.1 0 0.03); x=fminbn(l('-0,l*x*(l-x)0,l) y=0.1*x*(l-x)hold onplot(0,x,y,y/:x,x,0,y,':'); text(O,y/(dl/dt)iiV/VertlcalAIIgninenf/bottom*); text(Xr0.001,nuin2str(x)/Ho

3、rizontalAhginnent','center*); tltleCSI 模型的 dl/clt-1 曲線J;xlabel('r); ylabel('dl/(lt'); hold off;提示：fplot, fininbnd, plot，text, title, xlabel1) 畫曲線圖用fplot函數(shù)，調(diào)用格式如下：fplotffunjims)fun必須為一個(gè)M文件的函數(shù)名或?qū)ψ兞縓的可執(zhí)行字符串。 若lims取xmhi xniax,則x軸被限制在此區(qū)間上。若 lims 取xmhi xmax vinin ymax,則 y 軸也被限制。本題可用&qu

4、ot;*fplotC().l*x*(l-x)0 1.10 0.03);2) 求最大值用求解邊界約束條件下的非線性最小化函數(shù)fininbnd.調(diào)用格式如下: x=finlnbnd(lunxl,x2)fun必須為一個(gè)M文#的函數(shù)名或?qū)ψ兞縓的可執(zhí)行字符串。返回自變董X在區(qū)間xl<x<x2上函數(shù)取最小值時(shí)的X值。本題可用x=finlnbiid('-0.1*x*(l-x)0,l)y=0.1*x*(l-x)3) 指示*大值坐標(biāo)用線性繪圖函數(shù)plot,調(diào)用格式如下：plot(xl，yV顏色 線型 數(shù)據(jù)點(diǎn)圖標(biāo)：x2,y2J顏色 線型 數(shù)據(jù)點(diǎn)圖標(biāo)，) 本題可用hold on; %在上面的同

5、一張圖上畫線(同坐標(biāo)系) plot(0,x,y,y,':x,xJ0,y/:');4園冨帥娠&使用文本標(biāo)注函數(shù)text,調(diào)用格式如下：格式1text(x,y,文本標(biāo)識(shí)內(nèi)容,*HorizoiitalAIIgninent','字符串 1J x,y給定標(biāo)注文本在圖中添加的位置。HorlzontalAllgnmenf為水平控制屬性，控制文本標(biāo)識(shí)起點(diǎn)位于點(diǎn)(x,y)同一 水平線上。，字符串1-為水平控制屬性值，取三個(gè)值之一：'left*,點(diǎn)(x,y)位于文本標(biāo)識(shí)的左邊。'center' , ；(x,y)位于文本標(biāo)識(shí)的中心點(diǎn)。Tight，,點(diǎn)(

6、x，y)位于文本標(biāo)識(shí)的右邊。格式2text(x,y,文本標(biāo)識(shí)內(nèi)容,'VerticaIAIIginnent '字符串 2")x,y給定標(biāo)注文本在圖中添加的位置。*Vertica!Ahgninenf為垂直控制屬性.控制文本標(biāo)識(shí)起點(diǎn)位于點(diǎn)(x,y)同一垂 直線上。*字將串r為垂直控制屬性值，取四個(gè)值之一：*inlddle Sop打'cap', 'baseline', 'bottom'(對(duì)應(yīng)位置可在命令窗 口應(yīng)用 確定)本題可用text(0,y,'(di/(it)in','VerticaLAngninen

7、t'/bottoin');text(x,-0-001,num2str(x)/HorizontalAllgninent'/center');5)坐標(biāo)軸標(biāo)注調(diào)用函數(shù) xiabeb vlabel 和 title本題可用tItleCSI 模型 di/dt-i 曲線J;xlabel(T);ylabel('dl/(lf);程序運(yùn)行結(jié)果(比較138圖2)：(在圖形窗口菜單選擇Edit/Copy Figure.復(fù)制圖形)S槓型的di/dH曲線0,030.025di/dt)m0.020,010.005§ 0.015P0.20.40.50.60.81.2 畫 it

8、 曲線圖 p 136-138求出微分方程的解析解，畫出H曲線(Z(0)=0.15,=0.2,=030)(見(jiàn)138 圖1比校)。歩考程序：% 5.1傳染病模型模型2% 文件名 J P 136flgLin% dVdt=ki(l-i), l(0)=i0dear; de;x=dsolveCDx=k*x*(l-x)-；x(0)=x0*)%求微分方程的解析解，為符號(hào)表達(dá)式 x0=0.15; k=0.2;%xi 對(duì)應(yīng) i. xiO 對(duì)應(yīng) iO, k 對(duì)應(yīng) att=0:0.1:30;%時(shí)間單位為天for s=l:length(tt)%x的表達(dá)式中沒(méi)有點(diǎn)運(yùn)算，按標(biāo)量運(yùn)算取值xx t=tt(s);xx(s)=ev

9、al(x);%給出xl0=0.2,k=0.2,t.求符號(hào)表達(dá)式xl的對(duì)應(yīng)值end %XX為復(fù)數(shù)表示plot(仏 XX)；axls(0 31 0 L1);皿嗽圖1 SI模型的卜t曲線xlabeUt (天)*);ylabel(4 (病人所占比例)提示J1) 求解徹分方程dsolve,見(jiàn)提示；2) ® 出卜f 曲線(i(0)=05,啟02,e030)用 for 循環(huán)，函數(shù) length, eval, plo(, axis, title, xiabet ylabek程序運(yùn)行結(jié)果(見(jiàn)138圖1):命令窗口中的結(jié)果；-1 / (expdog (xO - 1) /xO) - k«i) -

10、 0Warning: Imaginary parts of comptex X and/or Y arguments ignored > In d74Jim a12A» i圖形窗口中的結(jié)果(比校138圖1)：2.（編程）傳染病模型3 （SIS模型）巳知傳染病模型3 （SIS模型）：J = _兄訂j_（l一一）, /（0） = Z,at<T其中，訛）是'第t天病人在總?cè)藬?shù)中所占的比例。2是每個(gè)病人每天有效接觸的平均人數(shù)（日接觸*）。巾是初始時(shí)刻（/=0）病人的比例。CT是整個(gè)傳染期內(nèi)每個(gè)病人有效接觸的平均人數(shù)（接觸數(shù)）。di 2.1畫不7曲線圖P138"1

11、39I 取z=0.b 0=1.5,畫出如下所示的4 曲線圖。試編寫一個(gè)m文件來(lái)實(shí)現(xiàn)。at（在圖形窗口菜單選擇Edit/Copy Figure,復(fù)制圖形）SIS模型的di/dt-i曲線252050異1十1.54Wg。0.05 OJ 015020.250.30.350.4（注：P139 圖 3） 提示：di用fplot函數(shù)畫出的曲線圖;at在上圖上用plot函數(shù)畫一條過(guò)原點(diǎn)的水平線;用 title, xlabeU ylabel 標(biāo)注。編寫的M文件和運(yùn)行結(jié)果（見(jiàn)139圖3）：%傳染病模型3 (SIS模型)的d"dt-i曲線圖%文件名 J pBSngXin%k=0d, 0=1.5 dear;

12、 de;fplot(-0.1*x*(x-(M/1.5)0 04 -0.0005 0.003); hold on;plot(0,0.4,0,0,*:*);tltle( SIS 模型的 d"dti 曲線 J; xiabeicr);ylabel('di/<lt');2.2 畫 i"曲線圖 pl38139要求：求出微分方程的解析解垃）。取z=0<2, 0=3, Z=040, B出如下所示的圖形。 試編寫一個(gè)m文件來(lái)實(shí)現(xiàn)。090807、 r - - 6 5 4 3 O O O O (凰丑血百Y昶)一二 一 二 O2 10 O O1(0)=0.21-1/0|

13、(0)=0.9(注：P139 圖 4)其中藍(lán)色實(shí)線為/(0)=0.2時(shí)的Z曲線(第1條)；黑色虛點(diǎn)線為過(guò)點(diǎn)(0,1-1/0)的水平線(第2條)；紅色虛線為/(0)=0.9時(shí)的Z曲線(第3條)。提示圖例標(biāo)注可用kgend(4(0)=0.24-l/!o'i(0)=0.93;編寫的M文件和運(yùn)行結(jié)果(比較1391圖4)：解法一：程序： %傳染病模型3 (SIS模型)的卜t曲線圖%文件名S P 138flg4mdear; de;% X =0.2, 0=3. i(0)=02,x 代表 ix=dsolve( Dx=-0.2*x*(x-(M/3)'；x(0)=0.2)%求微分方程的解析解，為符

14、號(hào)表達(dá) 式tt=0:0.1:40;%時(shí)間單位為天for l=l:length(tt)t=tt(l);xx(l)=eval(x);endplot(仏 XX)；hold on;p lot(0,41LM/34-l/31/-k")； % X =02 o =3, i(0)=09 x=dsolve(*Dx=-0.2*x*(x-(l-l/3)'/x(0)=0.9*) tt=O:O.l:4O;%0f間單位為天for l=l:length(tt)t=tt(l); xx(l)=eval(x); endplot(tt,xx/:r*); legendCI(0)=0-2*；l-l/o*/l(0)=0.

15、9); axls(0 40 0 1);tItleC 圖 1 SI 模型的 i-t 曲線(X=0.2,a=3)3; xIabelCt (天)lylabefi (病人所占比例)')； 命令窗口的結(jié)果:圖形窗口的結(jié)果:0908076 5 4 3 O O O O (區(qū)衛(wèi)血密Y昶二- - 二 O 2 10 O Oi(0)=0.21-1/0i(0)=0.9解法二: 程序：%傳染病模型3 (SIS模型)的卜t曲線圖%文件名S P13811典.mclear; de;% X =0.2, o =3, X 代表 ix=dsolve( Dx=-0.2*x*(x-(M/3)'/x(0)=x0*)%求微分

16、方程的解析解，為符號(hào)表達(dá) 式tt=0:0.1:40;%時(shí)間單位為天for x0=0.2,0.9 % i(0)=0.2,0.9for t=ttxx(2-(x0=0.2),round(t/0.1)+l)=eval(x);endend legendCI(0)=0.24-l/o'i(0)=0.9); axls(0 40 0 1);(MW圖1 SI模型的卜t曲線(1=0.2,a=3)*); xlabelCt (天)J;ylabelCi (病人所占比例)命令窗口的結(jié)果:File Edi t Debug Sesktop Vindow HelpX =-2/(3*(exp(logC(3*xO - 2)/

17、(3*x0) - (2*t)/15) - D)Warning： Imaginary parts of complex X and/or Y arguments ignored> In Untitled at 15A»圖形窗口的結(jié)果:與解法一相同解法三: 程序%傳染病模型3 (SIS模型)的卜t曲線圖%文件名 J pl38fig4.indear; de;x=dsolve(*Dx=-lain*x*(x-(M/sl)'x(0)=xO')%求微分方程的解析解，為符號(hào)表達(dá) 式tt=0:0.1:40;%B|間單位為天lain=0.2; si=3; % X =0,2, o =

18、3, x 代表 ifor x0=0.2,0.9 %i(0)=0.2,0.9for t=ttxx(2-(x0=0.2),round(t/0.1)+l)=eval(x);endendplot(仏 legendCI(0)=0.2*M/o*/i(0)=0.9*);axls(0 40 0 1);Mier 圖 1 SI 模型的 i-t 曲線=02,<y=3y); xlabelCt (天)J;yhtbelCi (病人所占比例)')； 命令窗口的結(jié)果:(tan(<(3i*i)/2 - i/2)*(- <lanr*t)/si + (2*atan<(2*3i*xC>*i)/(

19、si - I) - i)*L)Z<si - I)- i)*(si - Marnios： I co aginary parts of ccmplcx X and>or Y arguments ignored> In Oruitled Bt 15A»圖形窗口的結(jié)果:與解法一相同3.（驗(yàn)證）傳染病模型4 （SIR模型）P140-141SIR模型的方程J-J-/(O) = /oat=-Asi 5(0) = So dt0設(shè)z=b A =0.3, /(0)=0.02, s(0)=0.98o輸入pl40的程序并運(yùn)行，結(jié)果與教#pl41的圖7和圖8比較。ode45, pause的用

20、法見(jiàn)提示。 2個(gè)M文件(jar 1401)和運(yùn)行結(jié)果(比較141圖7、圖8):%51傳染病模型 %文件名；ill.ni function y=ill(t,x) a=l; b=0.3; %X用a表不，卩用b表不汕x(l)*x(2)b*x(l),曲*(1)*(2)??； i 用 x(l)表示，s 用 x(2)表示 命令M文件:'模型4 (SIR模型)函數(shù)M文件：模型4 (SIR模型)%51傳染病模型-%文件名；pl40.ni dear; clc;ts=0:50;x0=0.02,0.98;t,x=ode45CiHts,x0); t,x plot(t,x(:,l),t,x(:,2), grid,

21、 pause% 圖 7 l(t), s圖形 plot(x(:,2),x(:,l);gM,%ffl 8is 圖形(相軌線) i(t),s圖形(比較141圖7)：0.90807060504030.20101 111111、/khh h'V4J/051015202530354045501-S圖形（相軌線）（比較141圖8）：0.35030.25020.15OJ0.050 (1)010.20.304050.60.70.80.9-4.（驗(yàn)證）人口指數(shù)增長(zhǎng)模型參數(shù)估計(jì)及結(jié)果分析（美國(guó)1790-2000 年人口）pl63164美國(guó)1790-2000年人口統(tǒng)計(jì)數(shù)據(jù)（以百萬(wàn)為單位）年1790180018

22、1018201830184018501860187018801890人口3.95.37.29.612.917.123.231.43&650.262.9年19001910192019301940195019601970198019902000人口76.092.0106.5123.213L7150.71793204.0226.5251.42814人口指數(shù)增長(zhǎng)模型：xt） =xoe"(1)用表中數(shù)據(jù)進(jìn)行數(shù)據(jù)擬合求參數(shù) 將兩邊取對(duì)數(shù)，可得y = rt + a其中，J = In X,« = Inxo,即 xo = exp(“)©采用線性最小二乘法進(jìn)行數(shù)據(jù)擬合，用MA

23、TLAB中的函數(shù)polyfH計(jì)算。 "為每10年估計(jì)的人口增長(zhǎng)率，X0為1790年估計(jì)的初始人口數(shù)) 以下是M文件：de; format compact;x=3.95.37.29.612.917423.231.4 38.650.262.976.092.01063123.2131.7».150.7 1793 204.0 226.5 251.4 281.4;t=O:length(x)-l; %t=0 為 1790 年，t=l 為 1800 年， y=10g(x)； %取 1790-2000 年的數(shù)據(jù)ra=polynt(t,y4);dispC用1790-2000年數(shù)據(jù)估計(jì)的參數(shù)為

24、：Jr=ra(l)xO=ex p(ra(2)EE岡（1）運(yùn)行程序并給出結(jié)果（見(jiàn)畫）:Couand VindovA»y-Tt（2）人口指數(shù)增長(zhǎng)模型計(jì)算結(jié)果與實(shí)際數(shù)據(jù)比較（數(shù)據(jù)表）以下是M文件:cic; format coin pactX二3.9537.29.612.917.123.2 31.4 3&650.262.976.092.0106.5123-2 131.7 150.7 179.3 204.0 226.5 251.4 2814;t=O:length(x)-l;y=log(x)； %取 1790-2000 年的數(shù)據(jù) ra=polynt(t,y4);r=ra(l); xO=e

25、xp(ra(2);x2=x0*exp(r*t);dIspC指數(shù)增長(zhǎng)模型擬合美國(guó)人口數(shù)據(jù)的結(jié)果(x2)*) format short g;(llsp(1790+10*t; X； round(10*x2)/10');(2)運(yùn)行程序并給出結(jié)果(見(jiàn)167表4中X2列):Elie lilt Jctuc lesklcp 世indo V Kelp指數(shù)增長(zhǎng)模型擬合美國(guó)人口數(shù)據(jù)的結(jié)杲(x2)179039618005. S7.418107. 29. 118209. 6IL 1183013. 913.6184017. 11&6185023220.3I86031424.9187038. 630.518

26、8050. 237,3189062945.719007655.919109263.41920106 583.71930123. 2102.51940131. 7125.51950150. 7153.6i960179. 31881970204230. 11980226 5281.71990251 4344-82000281 4422. 1A»07R（3）人口指數(shù)增長(zhǎng)模型計(jì)算結(jié)果與實(shí)際數(shù)據(jù)比較（擬合圖形） 以下是M文件：X二3.9537.29.612.917.1 23.2 31.4 3&650.262.976.092.0106.5 12X2 131.7 150.7 179.3 2

27、04.0 226.5 251.4 281.4;t=0:length(x)-l;y=log(x)； %取 1790-2000 年的數(shù)據(jù) ra=polynt(t,y4);r=ra(l); x0=exp(ra(2); t2=llnspace(0,length(x)-130);x2=x0*exp(r*t2);plot(t,x,'r+t2,x2/b*);mbC指數(shù)增長(zhǎng)模型擬合圖形')；（3）運(yùn)行程序并給出結(jié)果（見(jiàn)1681圖3（b）：5-（驗(yàn)證，編程）估計(jì)阻滯增長(zhǎng)模型的參數(shù)和繪制圖形pl65168美國(guó)1790-2000年人口統(tǒng)計(jì)數(shù)據(jù)（以百萬(wàn)為單位）年17901800181018201830

28、184018501860187018801890人口3.9537.29.612.917.123.231.438.650.262.9年19001910192019301940195019601970198019902000人口76.092.0106.5123-213L7150.71793204.0226.5251.42814人口阻滯增長(zhǎng)模型：Xx(t) =/”、-rtl+d-lez丿(1)(驗(yàn)證)用1860-1990年的數(shù)據(jù)擬合估計(jì)參數(shù) 用下面方程dx=r- sx.估計(jì)拱數(shù)幾 程序如下：%用數(shù)值微分的三點(diǎn)公式計(jì)算美國(guó)人口增長(zhǎng)率(/10年)clear; clc;537.29.6X二3.9537.2

29、9.612.917.1 23.2 314 3&650.262.976.092.0106.5 12X2 131.7.«150.7 179.3 204.0 226.5 25L4；% 1790-1990 年的人口 xx=x(l,8:en(l)；%取 1860-1990 年的數(shù)據(jù) len=length(xx); h=l;dx=-3*xx(l)+4*xx(2)-xx(3), xx(3:len)-xx( 1 :len-2),求數(shù)值微分 dx3*xx(len)-4*xx(len-l)+xx(len-2)/(2*h);v=dx7xx;sr=p olyflt(xx,yj); r=sr(2) x

30、m=r/sr(l)(1)運(yùn)行程序并給出結(jié)果(見(jiàn)168差異大.書上結(jié)果數(shù)據(jù) 處理過(guò))：r1 丿 CoBsand Vindov口間區(qū)fEm EditDy kt op Wtnd<w HelpY =0. 22182xm =418. 2BA»阿(2)(編程)計(jì)算阻滯增長(zhǎng)模型擬合美國(guó)人口數(shù)據(jù)(1790.1990)取參數(shù)r=02557，切=392.0886,x(=39,編寫程序計(jì)算17901990年的人口, 參考的輸出如下(第1列為年，第2列為實(shí)際人口，第3列為計(jì)算人口)：F)CoBSXftd Vinduvim劉PH gdil"iadw州令ans =17903-93918005-3

31、518107-26518209683183012910.7164017-1IS. 7185023-217. 5186031-4223167038626.3188050-235. 8189062.94519007656. 219109269. 71920106.58551930123-2103. 91940131-7124. 51950150-7147.21960179-3171. 31970204196. 2198022&5221.2199025L4245. 3A» 1(2)給出程序及其運(yùn)行結(jié)果(見(jiàn)11167表4中X列兒clear; clc;foriiiat coin pac

32、t;x=3.9537.29.612.9174 23.2 314 38.650.262.976.092.0106.5123.2131.7».150.7 179.3 204.0 226.5 251.4； % 1790-1990 年的人口t=O:length(x)-l;x0=3.9; r=0.2557; xin=392.0886; %代入?yún)?shù)值 xx=xm7(l+(xin/x0-l)*exp(-r*t); %計(jì)算人口 format short g;1790+10*t;x;round(10*xx)/10'F)CoBSXftd Vinduvim劉PH gdil"iadw州令a

33、ns =17903-93918005-3518107-26518209683183012910.7164017-1IS. 7185023-217. 5186031-4223167038626.3188050-235. 8189062.94519007656. 219109269. 71920106.58551930123-2103. 91940131-7124. 51950150-7147.21960179-3171. 31970204196. 2198022&5221.2199025L4245. 3A» 1(3)(編程)繪制阻滯增長(zhǎng)模型擬合圖形(1790-1990)取參數(shù)r

34、=02557皿廬392.0886,x()=39,編寫程序繪制17901990年的人口數(shù) 據(jù)擬合圖形。參考圖形如下，實(shí)線為模型曲線，+為實(shí)際數(shù)據(jù)值。300250200150100+/+ Z異/50+68101214161820十0 2(3)給出程序及其運(yùn)行結(jié)果(見(jiàn)168圖4)：dear; de; format coin pact;x=3.9537.29.612.9174 23.2 3143&6 50.2 62.9 76.0 92.0106.5123.2131.7 150.71793 204.0 226.5 251.4； % 1790-1990 年的人口t=0: length(x)-l;x

35、0=3.9; r=0.2557; xin=392.0886; t3=llnspace(0,Iength(x)-130)； x3=xm./(l+(xin/x0-l)*ex p( r*t3);p Iot(t,x,'r+t3,x3,'b-*);OAA4yu250+200-+150-十七+100+/+/+50+-亠十+"AV<)246810121416 18 20附1:實(shí)驗(yàn)提示第12題提示：求解微分方程dsolve1)求解微分方程常微分方程符號(hào)解用函數(shù)dsoNc,調(diào)用格式如下：dsolvereqiil； 'equ2'9J變量名 J以代表微分方程及初始條件

36、的符號(hào)方程為輸入?yún)?shù)，多個(gè)方程或初始條件 可在一個(gè)輸入變量?jī)?nèi)聯(lián)立輸入，且以逗號(hào)分隔。默認(rèn)的獨(dú)立變*為f,也可把f變?yōu)槠渌姆?hào)變量。 字符D代表對(duì)獨(dú)立變*的微分，通常指d/dt。本題可用x=dsolve('Dx=k*x*( 1-X) *x(O)=xO') 2)畫出卜f 曲線(i(O)=OJ5M=O2,eO3O)用 for 循環(huán)，函數(shù) length, eval, plot, axis, title, xlabel, ylabeU第3題ode45, pause 的用法1) 求解微分方程的數(shù)值解函數(shù)ode45,格式如下：t,x=ode45(*funts,x0)fun是由一個(gè)或多個(gè)待解

37、方程寫成的函數(shù)式m文件：ts=to,tn表示此微分方程的積分限是從W到tf,也可以是一些離散的點(diǎn)，形 式為 ts=tO,tl,.,tf；XO為初值條件02) 等待用戶貝應(yīng)命令pause：程序執(zhí)行到該命令時(shí)暫停，直到用戶按任意 鍵后繼續(xù)(處在命令窗口有效)。微分方程模型當(dāng)我們描述實(shí)際對(duì)象的某些特性隨時(shí)間（或空間）而演變的過(guò)程、分析它的 變化規(guī)律、預(yù)測(cè)它的未來(lái)性態(tài)、研究它的控制手段時(shí)，通常要建立對(duì)象的動(dòng)態(tài)模 型建模時(shí)首先要根據(jù)建模目的和對(duì)問(wèn)題的具體分析作出簡(jiǎn)牝假設(shè)，然后按照對(duì) 象內(nèi)在的或可以類比的其他對(duì)象的規(guī)律列出微分方程，求出方程的解并將結(jié)果 翻譯回實(shí)際對(duì)象，就可以進(jìn)行描述、分析、預(yù)測(cè)或控制了.

38、我們?cè)?.5節(jié)建立的胃 腸道和血液中藥量變化的模型，大致就是這個(gè)過(guò)程.事實(shí)上，在微分方程課程中，解所謂應(yīng)用題時(shí)我們已經(jīng)遇到簡(jiǎn)單的建立動(dòng)態(tài) 模型問(wèn)題，例如“一質(zhì)量為m的物體自高血處自由下落，初速是0,設(shè)阻力與下落 速度的平方成正比，比例系數(shù)為筑求下落速度隨時(shí)間的變化規(guī)律”，又如“容器 內(nèi)有鹽水100 1,內(nèi)含鹽10燉，今以3 1/min的速度從一管放進(jìn)凈水，以2 1/min的 速度從另一管抽出鹽水，設(shè)容器內(nèi)鹽水濃度始終是均勻的，求容器內(nèi)含鹽量隨時(shí) 間變化的規(guī)律”本章討論的動(dòng)態(tài)模型與這些問(wèn)題的主要區(qū)別是，所謂微分方程 應(yīng)用題大多是物理或幾何方面的典型問(wèn)題，假設(shè)條件已經(jīng)給出，只需用數(shù)學(xué)符號(hào) 將已知規(guī)

39、律表示出來(lái)，即可列出方程，求解的結(jié)果就是問(wèn)題的答案，答案是惟一 的，已經(jīng)確定的.而本章的模型主要是非物理領(lǐng)域的實(shí)際問(wèn)題，要分析具體情況 或進(jìn)行類比才能給出假設(shè)條件作出不同的假設(shè)，就得到不同的方程，所以是事 先沒(méi)有答案的求解結(jié)果還要用來(lái)解釋實(shí)際現(xiàn)象并接受檢驗(yàn)-隨著衛(wèi)生設(shè)施的改善、醫(yī)療水平的提高以及人類文明的不斷發(fā)展，諸如耀 亂、天花等曾經(jīng)肆虐全球的傳染性疾病已經(jīng)得到有效的控制但是一些新的、不 斷變異著的傳染病毒卻悄悄向人類襲來(lái) 20世紀(jì)80年代，十分險(xiǎn)惡的艾滋病毒 開始肆慮全球，至今仍在蔓延；2003年春來(lái)歷不明的SARS病毒突襲人間，給人 們的生命財(cái)產(chǎn)帶來(lái)極大的危害長(zhǎng)期以來(lái)，建立傳染病的數(shù)學(xué)模

40、型來(lái)描述傳染病 的傳播過(guò)程,分析受感染人數(shù)的變化規(guī)律，探索制止傳染病蔓延的手段等，一宜 是各國(guó)有關(guān)專家和官員關(guān)注的課題.分析各種傳染病的傳播，而只是按不同類型傳染病的傳播過(guò)程有其各自不同的特點(diǎn)，弄清這些待點(diǎn)需要相當(dāng) 多的病理知識(shí)，這里不可能從醫(yī)學(xué)的角度 照一般的傳播機(jī)理建立幾種模型".模型1在這個(gè)最簡(jiǎn)單的模型中，設(shè)時(shí)刻£的病人人數(shù)兀0)是連續(xù)、可微函 數(shù),并且每天每個(gè)病人有效接觸(足以使人致病的接觸)的人數(shù)為常數(shù)入，考察上 到t +山病人人數(shù)的增加，就有%(/ + )x(t)再設(shè)1=0時(shí)有孔個(gè)病人，即得微分方程(1)誓=入兀,*(0) =3：0方稈(1)的解為兀(上)=欣&

41、quot;(2)結(jié)果表明，隨著十的增加，病人人數(shù)戈(/)無(wú)限增長(zhǎng)，這顯然是不符合實(shí)際的.建模失敗的原因在于：在病人有效接觸的人群中，有健康人也有病人， 而其中只有健康人才可以被傳染為病人，所以在改進(jìn)的模型中必須區(qū)別這 兩種人.模型2(SI模型)假設(shè)條件為1. 在疾病傳播期內(nèi)所考察地區(qū)的總?cè)藬?shù)N不變，既不考慮生死，也不考慮 遷移.人群分為易感染者(Susceptible)和已感染者(Infective)兩類(取兩個(gè)詞的 第1個(gè)字母，稱之為SI模型)，以下簡(jiǎn)稱健康者和病人.時(shí)刻：這兩類人在總?cè)?數(shù)中所占的比例分別記作s(f)和i(t).2. 每個(gè)病人每天有效接觸的平均人數(shù)是常數(shù)入，入稱為日接觸率.

42、當(dāng)病人與 健康者有效接觸時(shí),使健康者受感染變?yōu)椴∪?根據(jù)假設(shè)，每個(gè)病人每天可使入$(£)個(gè)健康者變?yōu)椴∪?，因?yàn)椴∪藬?shù)為 Ni(t)，所以每天共有入個(gè)健康者被感染，于是入Nd就是病人數(shù)M的 增加率，即有nJ"曲乂因?yàn)?再記初始時(shí)刻(20)病人的比例為d則 = Az(l - j),/(O) =/o方程(5)是logistic模型(詳細(xì)描述參見(jiàn)5. 6節(jié))它的解為137139 趙 2 1、2 2?累(6)iO)=:i(t) T和加i的圖形如圖1和圖2所示.139圖1 SI模型的iT曲線圖2 SI模型的半曲線at由(5),(6)式及圖1、圖2可知，第一，當(dāng)1 = 1/2時(shí)專達(dá)到最大

43、值(気)，這個(gè)時(shí)刻為(7)«« = A " 1"亍 1 %這時(shí)病人增加得最快，可以認(rèn)為是醫(yī)院的門診量最大的一夭，預(yù)示著傳染病高潮 的到來(lái)，是醫(yī)療衛(wèi)生部門關(guān)注的時(shí)刻-與入成反比，因?yàn)槿战佑|率入表示該地 區(qū)的衛(wèi)生水平，入越小衛(wèi)生水平越高.所以改善保健設(shè)施、提高衛(wèi)生水平可以推 遲傳染病高潮的到來(lái)第二，當(dāng)£-8時(shí)if 1,即所有人終將被傳染，全變?yōu)椴∪? 這顯然不符合實(shí)際情況.其原因是模型中沒(méi)有考慮到病人可以治愈，人群中的健 康者只能變成病人，病人不會(huì)再變成健康者為了修正上述結(jié)果，必須重新考慮模型的假設(shè)，下面兩個(gè)模型中我們討論病 人可以治愈的情況.橫型

44、3( SIS模型)有些傳染病如傷鳳、痢疾等愈后免疫力很低，可以假定 無(wú)免疫性，于是病人被治愈后變成健康者，健康者還可以被感染再變成病人，所 以這個(gè)模型稱SIS模型.SIS模型的假設(shè)條件1,2與SI模型相同，增加的條件為3. 每天被治愈的病人數(shù)占病人總數(shù)的比例為常數(shù)稱為日治敵率病人 治愈后成為仍可被感染的健康者顯然l/p是這種傳染病的平均傳染期不難看出，考慮到假設(shè)3,SI模型的(3)式應(yīng)修正為"半=入皿-屮(8)d£（4）式不變，于是（5）式應(yīng)改為 = Ai（l -£）-fjLi, i（0） =io（9）我們不去求解方程（9）（雖然它的解可以解析地表出），而是通過(guò)圖

45、形分析 i（t）的變化規(guī)律定義CT = A/£（10）注意到入和1么的含義，可知"是整個(gè)傳染期內(nèi)每個(gè)病人有效接觸的平均人 數(shù),稱為接觸數(shù).利用6方程（9）可以改寫作(11)由方程（11）容易先畫出尋i的圖形（圖3,圖5）,再畫出i 7的圖形（圖4,其中虛線是的情況 &圖5 SIS模型的字"曲線（crWl）0（圖6 SIS模型的i-t曲線（"Ml）不難看岀，接觸數(shù)7 = 1是一個(gè)網(wǎng)值.當(dāng)a>i時(shí)i（0的增減性取決于的大?。ㄒ?jiàn)圖4）,但其極限值i（8）=1 -丄隨b的增加而增加（試從的含義給<7110 趙 3仔予解釋)；當(dāng)時(shí)病人比例i(t

46、)越來(lái)越小，最終趨于0,這是由于傳染期內(nèi)經(jīng) 有效接觸從而使健康者變成的病人數(shù)不超過(guò)原來(lái)病人數(shù)的緣故-SI模型可視為本模型的特例，請(qǐng)讀者考慮它相當(dāng)于本模型中H或C取何值 的情況.模型4(SIR模型)大多數(shù)傳染病如天花、流感、肝炎、麻疹等治愈后均有 很強(qiáng)的免疫力，所以病愈的人既非健康者(易感染者)，也非病人(已感染者)，他 們已經(jīng)退岀傳染系統(tǒng).這種情況比較復(fù)雜，下面將詳細(xì)分析建模過(guò)程.橫型假設(shè)1. 總?cè)藬?shù)N不變?nèi)巳悍譃榻】嫡?、病人和病愈免疫的移出?Removed) 三類，稱SIR模型.時(shí)刻/三類人在總?cè)藬?shù)N中占的比例分別記作$ (f),i(C和 r(t).2. 病人的日接觸率為人，日治愈率為“(

47、與S1模型相同)，傳染期接觸數(shù)為 (T = A/i.樓型構(gòu)成由假設(shè)1顯然有(12)(13)$ (£)+ i( f) + r( t) =1根據(jù)條件2方程(8)仍成立.對(duì)于病愈免疫的移出者而育應(yīng)有N半=沏dz(14)再記初始時(shí)刻的健康者和病人的比例分別是旬( >0)和，0(4>0)(不妨設(shè)移 出者的初始值心=0),則由(8),(12),(13)式，SIR模型的方程可以寫作 r竽 m -皿，i(o)= io I 竽=-A5I, $(O) =Sa方程(14)無(wú)法求出sQ)和i(門的解析解，我們先作數(shù)值計(jì)算數(shù)值計(jì)算 在方程(14)中設(shè)入=0. 3,£(0) =0. 02,

48、$(0) =0.98,用 MATLAB軟件編程：function y = ill(t,x)a = 1 ；b = 0 <3 ；y = a x(l) * x(2 ) -b * x(l) , - a * x(l ) *x(2)'ts = 0 ：50 ;xO 二0 .02 ,0 .98;t ,x = ode45 Cill", ts ,xO ) ； t ,x#111 題3存案141plot(t,x( : ,1) ,t ,x( : ,2),grid,pauseplot(x( ： ,2 ) ,x( ： ,1 ) ) ,grid,輸出的簡(jiǎn)明計(jì)算結(jié)果列入表1，心)，$a)的圖形見(jiàn)圖7,圖

49、8是 S 的® 形，稱為相軌線，初值i(0) =0.02,$(0) =0.98相當(dāng)于圖8中的耳點(diǎn)，隨著/的 增加，($/)沿軌線自右向左運(yùn)動(dòng)由表1、圖7、圖8可以看出，認(rèn)門由初值增長(zhǎng)至 約« = 7時(shí)達(dá)到最大值，然后減少，£ oo,iTO；$a)則單調(diào)滅少， t80- 039 8.t012345678i(t)0. 020 00. 039 00. 073 20- 128 50. 203 30.279 50. 331 20. 344 40. 324 70. 980 00 952 50.901 90.816 90 692 70, 543 80 399 50. 283 9

50、0. 202 7t91015202530354045i(t)0 286 30.241 80 078 70. 022 30.006 10.001 70.000 50,000 I0$0)0. 149 30. 114 50. 054 30. 043 40. 040 80. 040 I0 039 90. 039 90. 039 8表1的數(shù)值計(jì)算結(jié)果0.4/、/、/r0-30.20-1圖8 iT圖形(相軌線)0 020.6 0.8圖7i(t)$(O圖形為了分析的一般變化規(guī)律，需要進(jìn)行相軌線分析-相軌線分析 我們?cè)跀?shù)值計(jì)算和圖形觀察的基礎(chǔ)上，利用相軌線討論解 i(t) ,$Q)的性質(zhì).s-i平面稱為相平面

51、，相軌線在相平面上的定義域(sM) eO為D = ) (s,i) l$m0,iM0,s + iWl 1在方程(14)中消去dr,并注意到a-的定義(10),可得di 1=ds as容易求出方程(16)的解為"10 = 0"(”)匚1叱(15)(16)(17)在定義域D內(nèi)，（17）式表示的曲線即為相軌線，如圖9所示，其中箭頭表示了隨 著時(shí)間t的增加$（£）和i（t）的變化趨向.下面根據(jù)（14）, （17）式和圖9分析s（r）, iQ）和rU）的變化情況（W«時(shí)它們的極限值 分別記作九，i和r.）.1. 不論初始條件如何，病人終將消 =0(18)失，即 其證明如下.圖9 SIR模型的相軌線首先，由（普wO,而故九存在；由而心）W1,故口存在;再由（12）知i存在其次，若is =£>0,則由（12）,對(duì)于充分大的£有竽 My,這將導(dǎo)致/= 8 ,與九存在相矛盾.從圖形上看，不論相軌線從P|或從巴點(diǎn)出發(fā)，它終將與$軸相交a充 分大）.2. 最終未被感染的健康者的比例是兒，在（17）式中令i=0得到，九是方程(19)1 1九 nSo + 切一$ + In =0