常微分方程模型

1、 微微分方程的基本概念分方程的基本概念 微分方程的解析解與數(shù)值解微分方程的解析解與數(shù)值解建立微分方程模型的基本方法建立微分方程模型的基本方法典型的建模案例簡(jiǎn)介典型的建模案例簡(jiǎn)介 在研究實(shí)際問題時(shí)，常常會(huì)聯(lián)系到某些變量的在研究實(shí)際問題時(shí)，常常會(huì)聯(lián)系到某些變量的或或，這樣所得到變量之間的關(guān)系式就是，這樣所得到變量之間的關(guān)系式就是微分方程模型。微分方程模型。 微分方程模型在自然科學(xué)中的應(yīng)用主要以物理、微分方程模型在自然科學(xué)中的應(yīng)用主要以物理、力學(xué)等客觀規(guī)律為基礎(chǔ)建立起來，而在經(jīng)濟(jì)學(xué)、人力學(xué)等客觀規(guī)律為基礎(chǔ)建立起來，而在經(jīng)濟(jì)學(xué)、人口預(yù)測(cè)等社會(huì)科學(xué)方面的應(yīng)用則是在類比、假設(shè)等口預(yù)測(cè)等社會(huì)科學(xué)方面的應(yīng)用則

2、是在類比、假設(shè)等措施下建立起來。措施下建立起來。 微分方程的基本概念微分方程的基本概念 微分方程的解析解與數(shù)值解微分方程的解析解與數(shù)值解建立微分方程模型的基本方法建立微分方程模型的基本方法典型的建模案例簡(jiǎn)介典型的建模案例簡(jiǎn)介 求微分方程（組）的求微分方程（組）的解析解解析解命令命令: :記號(hào)記號(hào): : 在表達(dá)微分方程時(shí)，用字母在表達(dá)微分方程時(shí)，用字母D表示求微分，表示求微分，D2、D3等表示求高階微分等表示求高階微分. .任何任何D后所跟的字母為后所跟的字母為因變量，自變量可以指定或由系統(tǒng)規(guī)則選定為確省。因變量，自變量可以指定或由系統(tǒng)規(guī)則選定為確省。例如，微分方程例如，微分方程 022 dxy

3、d應(yīng)表達(dá)為：應(yīng)表達(dá)為：dsolve(D2y=0)例例1 1 求求 21udtdu 的通解的通解. .解解 輸入命令：輸入命令：dsolve(Du=1+u2,t)例例2 2 求微分方程的特解。求微分方程的特解。 15)0( , 0)0(029422yyydxdydxyd 解解: : 輸入命令輸入命令: : y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x)結(jié)結(jié) 果果 為為 : y =3e-2xsin（5x） 例例3 3 求微分方程組的通解求微分方程組的通解. . zyxdtdzzyxdtdyzyxdtdx244354332解解 輸入命令輸入命令 ：x,y,z=

4、dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z,t); 結(jié)結(jié) 果果 為：為：x =C2*exp(-t)+C3*exp(2*t) y =C1*exp(-2*t)+C2*exp(-t)+C3*exp(2*t) z =C1*exp(-2*t)+C3*exp(2*t)微分方程的數(shù)值解微分方程的數(shù)值解（一）常微分方程數(shù)值解的定義（一）常微分方程數(shù)值解的定義 在生產(chǎn)和科研中所處理的微分方程往往很復(fù)在生產(chǎn)和科研中所處理的微分方程往往很復(fù)雜且大多得不出一般解。而在實(shí)際上對(duì)初值問題，雜且大多得不出一般解。而在實(shí)際上對(duì)初值問題，一般是要求得到解在若干個(gè)點(diǎn)上滿

5、足規(guī)定精確度一般是要求得到解在若干個(gè)點(diǎn)上滿足規(guī)定精確度的近似值，或者得到一個(gè)滿足精確度要求的便于的近似值，或者得到一個(gè)滿足精確度要求的便于計(jì)算的表達(dá)式。計(jì)算的表達(dá)式。用用Matlab軟件求常微分方程的數(shù)值解軟件求常微分方程的數(shù)值解t，x=solver（f,ts,x0,options）ode45 ode23 ode113ode15sode23s由待解由待解方程寫方程寫成的成的m-m-文件名文件名ts=t0，tf，t0、tf為自為自變量的初變量的初值和終值值和終值函數(shù)的函數(shù)的初值初值ode23：組合的組合的2/3階龍格階龍格- -庫(kù)塔算法庫(kù)塔算法ode45：運(yùn)用組合的運(yùn)用組合的4/5階龍格階龍格-

6、 -庫(kù)塔算法庫(kù)塔算法自變自變量值量值函數(shù)函數(shù)值值用于設(shè)定誤差限用于設(shè)定誤差限( (缺省時(shí)設(shè)定相對(duì)誤差缺省時(shí)設(shè)定相對(duì)誤差10-3, , 絕對(duì)誤差絕對(duì)誤差10-6),命令為：命令為：options=odeset（reltol,rt,abstol,at）, rt，at：分別為設(shè)定的相對(duì)誤差和絕對(duì)誤差分別為設(shè)定的相對(duì)誤差和絕對(duì)誤差. .在解在解n個(gè)未知函數(shù)的方程組時(shí)，個(gè)未知函數(shù)的方程組時(shí)，x0和和x均為均為n維向維向量，量，m-文件中的待解方程組應(yīng)以文件中的待解方程組應(yīng)以x的分量形式寫成。的分量形式寫成。使用使用Matlab軟件求數(shù)值解時(shí)，高階微分方程必軟件求數(shù)值解時(shí)，高階微分方程必須等價(jià)地變換成一階

7、微分方程組。須等價(jià)地變換成一階微分方程組。注意注意: :例例4 4)1(0)0( ; 2)0(0)1(1000222 xxxdtdxxdtxd解解: : 令令 y1=x，y2=y1 0)0(, 2)0()1(1000)1(211221221yyyyyyyy使用使用Matlab軟件求數(shù)值解時(shí)，高階微分方程必須等價(jià)軟件求數(shù)值解時(shí)，高階微分方程必須等價(jià)地變換成一階微分方程組地變換成一階微分方程組。1 1、建立、建立m-m-文件文件vdp1000.m如下：如下： function dy=vdp1000(t,y)%y(1)=y1,y(2)=y2 dy=zeros(2,1);%dy=(0,0)T dy(1

8、)=y(2); dy(2)=1000*(1-y(1)2)*y(2)-y(1); 2 2、取、取t0=0，tf=3000，輸入命令：，輸入命令： T,Y=ode15s(vdp1000,0 3000,2 0); plot(T,Y(:,1),-)1222121121000(1)(0)2,(0)0yyyyyyyy 例例5 5 求解方程求解方程()dxxyzdtdyyzdtdzxyyzdt 其中其中 ，初始條件為，初始條件為9,30,12,4(0)1, (0)0, (0)0.002.xyz1 1、建立、建立m-m-文件文件lorenz.m如下：如下： function dy=lorenz(t,y)%y(

9、1)=x,y(2)=y,y(3)=z beta=9/4;ruo=30;segma=12; dy=zeros(3,1); dy(1)=-beta*y(1)+y(2)*y(3); dy(2)=-segma*(y(2)-y(3); dy(3)=-y(1)*y(2)+ruo*y(2)-y(3);()dxxyzdtdyyzdtdzxyyzdt 其中其中 ，初始條件為，初始條件為9,30,12,4(0)1, (0)0, (0)0.002.xyz2 2、取、取t0=0，tf=100，輸入命令：，輸入命令： T,Y=ode15s(lorenz,0 100,1 0 0.002); plot(T,Y(:,1),

10、+) plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+)練習(xí)：練習(xí)：1：求下列方程的解析解：求下列方程的解析解：2(1)46ln ;t xtxxtt2(2), (0)1,(/ )0ya y yya 2：求下面方程的數(shù)值解：求下面方程的數(shù)值解2(3), (0)1,(/ )0,4ya y yyaa 123213312123(4)0.51(0)0,(0)1,(0)1yy yyy yyy yyyy 微分方程的基本概念微分方程的基本概念 微分方程的解析解與數(shù)值解微分方程的解析解與數(shù)值解建立微分方程模型的基本方法建立微分方程模型的基本方法典型的建模案例簡(jiǎn)介典型的建模案例簡(jiǎn)介1 1

11、. .按規(guī)律直接列方程按規(guī)律直接列方程在數(shù)學(xué)、力學(xué)、物理、化學(xué)等學(xué)科中許多自然現(xiàn)象所在數(shù)學(xué)、力學(xué)、物理、化學(xué)等學(xué)科中許多自然現(xiàn)象所滿足的規(guī)律已為人們所熟悉，并直接由微分方程所描滿足的規(guī)律已為人們所熟悉，并直接由微分方程所描述。如牛頓第二定律、放射性物質(zhì)的放射性規(guī)律等。述。如牛頓第二定律、放射性物質(zhì)的放射性規(guī)律等。我們常利用這些規(guī)律對(duì)某些實(shí)際問題列出微分方程我們常利用這些規(guī)律對(duì)某些實(shí)際問題列出微分方程。 一個(gè)較熱的物體置于室溫為一個(gè)較熱的物體置于室溫為180c的房間內(nèi)，該物體最的房間內(nèi)，該物體最初的溫度是初的溫度是600c，3分鐘以后降到分鐘以后降到500c 。想知道它的想知道它的溫度降到溫度降

12、到300c 需要多少時(shí)間需要多少時(shí)間？10分鐘以后它的溫度分鐘以后它的溫度是多少？是多少？牛頓冷卻（加熱）定律：牛頓冷卻（加熱）定律：將溫度為將溫度為T的物體放入處于的物體放入處于常溫常溫 m 的介質(zhì)中時(shí)，的介質(zhì)中時(shí)，T的變化速率正比于的變化速率正比于T與周圍與周圍介質(zhì)的溫度差。介質(zhì)的溫度差。物體散熱時(shí)間的確定物體散熱時(shí)間的確定 建模范例建模范例1 1分析：分析：假設(shè)房間足夠大，放入溫度較低或較高的物體時(shí)，假設(shè)房間足夠大，放入溫度較低或較高的物體時(shí)，室內(nèi)溫度基本不受影響，即室溫分布均衡室內(nèi)溫度基本不受影響，即室溫分布均衡, ,保持為保持為m。 建立模型：建立模型：設(shè)物體在冷卻過程中的溫度為設(shè)物

13、體在冷卻過程中的溫度為T(t)，t0， 轉(zhuǎn)化為?dTTdt 與與成成 正正 比比成正比成正比與與mTdtdT 牛頓冷卻（加熱）定律：牛頓冷卻（加熱）定律：將溫度為將溫度為T的物體放入處于的物體放入處于常溫常溫 m 的介質(zhì)中時(shí)，的介質(zhì)中時(shí)，T的變化速率正比于的變化速率正比于T與周圍與周圍介質(zhì)的溫度差。介質(zhì)的溫度差。()dTk Tmdt (),(0)60,(3)50.dTk TmdtTT 建立微分方程建立微分方程其中參數(shù)其中參數(shù)k 0，m=18. 求得一般解為求得一般解為18,0,ktTcet (T=dsolve(DT=-k*(T-m) , t)代入條件：代入條件：T(0)=60T(0)=60 求

14、得求得c=42(T=dsolve(DT=-k*(T-m) , T(0)=60, t)代入條件代入條件T(3)=50T(3)=50*50501842ke (solve(50=18+42*exp(-k*3) , k)116ln321k 結(jié)果結(jié)果 ：T(10)=18+42 = 34.9662 0，102116ln31 e T(t)=18+42 , t 0. te2116ln31問題問題: （1）10分鐘以后它的溫度是多少？分鐘以后它的溫度是多少？問題問題: （2）溫度降到溫度降到300c需要多少時(shí)間？需要多少時(shí)間？ 30=18+42 , t 0. te2116ln31該物體溫度降至該物體溫度降至30

15、0c 需要需要13.8224分鐘。分鐘。( solve(30=18+42*exp(1/3*log(16/21)*t), t) ) 作業(yè)作業(yè)2(刑事偵查中死亡時(shí)間的鑒定刑事偵查中死亡時(shí)間的鑒定) 警察于清晨警察于清晨7:107:10在一所住宅內(nèi)發(fā)現(xiàn)一具尸體，在一所住宅內(nèi)發(fā)現(xiàn)一具尸體，測(cè)得尸體的溫度是測(cè)得尸體的溫度是25250 0C C，當(dāng)時(shí)環(huán)境溫度是，當(dāng)時(shí)環(huán)境溫度是20200 0C C，一，一小時(shí)后再次測(cè)溫尸體溫度下降為小時(shí)后再次測(cè)溫尸體溫度下降為22220 0C C，若人的正常，若人的正常體溫是體溫是37370 0C C，請(qǐng)協(xié)助警察估計(jì)死者的死亡時(shí)間。，請(qǐng)協(xié)助警察估計(jì)死者的死亡時(shí)間。牛頓冷卻（

16、加熱）定律：牛頓冷卻（加熱）定律：將溫度為將溫度為T的物體放入處于常溫的物體放入處于常溫 m 的介質(zhì)中時(shí)，的介質(zhì)中時(shí)，T的變化速率正比于的變化速率正比于T與周圍介質(zhì)的溫度差。與周圍介質(zhì)的溫度差。設(shè)設(shè)t t時(shí)刻尸體的溫度為時(shí)刻尸體的溫度為T(t)，死亡死亡時(shí)間時(shí)間t =0t=t0發(fā)現(xiàn)尸發(fā)現(xiàn)尸體的時(shí)體的時(shí)間間t=t0+1solve()注求解方程組的命令：注求解方程組的命令：solve例：例： 2222119119xyxy syms x yS1=(x-1)2+(y-1)2-9;S2=(x+1)2+(y+1)2-9;x,y=solve(S1,S2)一場(chǎng)筆墨官司（放射性廢物的一場(chǎng)筆墨官司（放射性廢物的處

17、理問題）處理問題） 美國(guó)原子能委員會(huì)（現(xiàn)為核管理委員會(huì)）處理濃縮美國(guó)原子能委員會(huì)（現(xiàn)為核管理委員會(huì)）處理濃縮放射性廢物，是將廢物放入密封性能很好的圓桶中，然放射性廢物，是將廢物放入密封性能很好的圓桶中，然后扔到水深后扔到水深300300英尺的海里英尺的海里. .他們這種做法他們這種做法安全安全嗎？嗎？分析：分析：可從各個(gè)角度去分析造成危險(xiǎn)的因素，這里可從各個(gè)角度去分析造成危險(xiǎn)的因素，這里僅考慮圓桶泄露的可能。僅考慮圓桶泄露的可能。聯(lián)想：聯(lián)想：安全安全 、危險(xiǎn)、危險(xiǎn)建模范例建模范例2 21 1. .按規(guī)律直接列方程按規(guī)律直接列方程問題的關(guān)鍵問題的關(guān)鍵l圓桶至多能承受多大的圓桶至多能承受多大的沖撞

18、速度沖撞速度？(40(40英尺英尺/ /秒秒););l圓桶和海底碰撞時(shí)的速度有多大？圓桶和海底碰撞時(shí)的速度有多大？ 新問題：新問題：求這一種桶沉入求這一種桶沉入300300英尺的海底時(shí)英尺的海底時(shí)的末速度。（原問題是什么的末速度。（原問題是什么? ?）可利用的數(shù)據(jù)條件：可利用的數(shù)據(jù)條件： 圓桶的總重量圓桶的總重量 W=527.327（磅）（磅） 圓桶受到的浮力圓桶受到的浮力 B =470.327（磅）（磅） 圓桶下沉?xí)r受到的海水阻力圓桶下沉?xí)r受到的海水阻力 D=Cv，C=0.08 可利用牛頓第二定律，建立圓桶下沉位移滿足的可利用牛頓第二定律，建立圓桶下沉位移滿足的微分方程：微分方程： )1(2

19、2DBWdtydm vdtdyCvDgwm ,其中其中22d ymWBDdt ,WdymDCvvgdt 其其 中中W dvCvWBg dt 已知：已知：求求 ( )v t變形變形()dvCggvWBdtWW (0)0V (),(2)(0)0.dvcggvWBdtWWV 方程的解為方程的解為0),1()( teCBWtvtWCg計(jì)算碰撞速度，需確定圓桶和海底的碰撞時(shí)間計(jì)算碰撞速度，需確定圓桶和海底的碰撞時(shí)間t0 ?原問題得到解決了嗎原問題得到解決了嗎? ?解決思路：解決思路：避開求避開求t0的難點(diǎn)的難點(diǎn) 令令 v(t)=v(y(t), 其中其中 y=y(t) 是圓桶下沉深度。是圓桶下沉深度。 d

20、tdydydvdtdv. 將將代入代入22d ymWBDdt dvmWBDdt ,.CvBWdtdydydvm .,W dvvWBCvg dy vdvgWBCv dyW (0)0, (0)0.vy ,(0 )0 ,(0 )0 .vd vgWBC v d yWvy 兩邊積分得函數(shù)方程：兩邊積分得函數(shù)方程： ,ln2WgyBWCvBWCBWCv 若能求出函數(shù)若能求出函數(shù)v=v(y), ,就可求出碰撞速度就可求出碰撞速度v(300). .,ln2WgyBWCvBWCBWCv v=v(y) 是一個(gè)單調(diào)上升函數(shù)，而是一個(gè)單調(diào)上升函數(shù)，而v 增大增大, ,y 也增大也增大, ,其反函數(shù)其反函數(shù) y = y

21、(v)也是一個(gè)增函數(shù)。也是一個(gè)增函數(shù)。 ),ln(2BWCvBWCBWCWgWy 令令 v=40( (英尺英尺/ /秒秒) )，g=32.2( (英尺英尺/ /秒秒2 2),),算出算出y= 238.4 ( (英尺英尺) )300( (英尺英尺) )問題的實(shí)際解答：?jiǎn)栴}的實(shí)際解答： 美國(guó)原子能委員會(huì)處理放射性廢物的做美國(guó)原子能委員會(huì)處理放射性廢物的做法是極其危險(xiǎn)的，法是極其危險(xiǎn)的，必須改變必須改變。 在生物、經(jīng)濟(jì)等學(xué)科中，許多現(xiàn)象所滿足的規(guī)律在生物、經(jīng)濟(jì)等學(xué)科中，許多現(xiàn)象所滿足的規(guī)律并不很清楚而且相當(dāng)復(fù)雜，因而需要根據(jù)實(shí)際資料或并不很清楚而且相當(dāng)復(fù)雜，因而需要根據(jù)實(shí)際資料或大量的實(shí)驗(yàn)數(shù)據(jù)，提出

22、各種假設(shè)。在一定的假設(shè)下，大量的實(shí)驗(yàn)數(shù)據(jù)，提出各種假設(shè)。在一定的假設(shè)下，給出實(shí)際現(xiàn)象所滿足的規(guī)律，然后利用適當(dāng)?shù)臄?shù)學(xué)方給出實(shí)際現(xiàn)象所滿足的規(guī)律，然后利用適當(dāng)?shù)臄?shù)學(xué)方法列出微分方程。法列出微分方程。 ( (交通管理中的黃燈問題交通管理中的黃燈問題) ) 在十字路口的交通管理中，亮紅燈之前，要在十字路口的交通管理中，亮紅燈之前，要亮一段時(shí)間黃燈，這是為了讓那些正行駛在十字路亮一段時(shí)間黃燈，這是為了讓那些正行駛在十字路口上或距十字路口太近以致無法停下來的車輛通過口上或距十字路口太近以致無法停下來的車輛通過路口。路口。建模范例建模范例3 3試建立合理的數(shù)學(xué)模型以確定黃燈應(yīng)該亮多長(zhǎng)時(shí)間。試建立合理的數(shù)學(xué)

23、模型以確定黃燈應(yīng)該亮多長(zhǎng)時(shí)間。(1) (1) 問題的分析問題的分析 在十字路口行駛的車輛中，主要因素是機(jī)動(dòng)車輛。在十字路口行駛的車輛中，主要因素是機(jī)動(dòng)車輛。 駛駛近交叉口的駕駛員，近交叉口的駕駛員， 在看到黃色信號(hào)要作出決定在看到黃色信號(hào)要作出決定: : 是停車是停車還是要通過十字路口。如果它按法定速度還是要通過十字路口。如果它按法定速度( (或低于法定速度或低于法定速度) )行駛，當(dāng)決定停車時(shí)，行駛，當(dāng)決定停車時(shí)，它必須有足夠的停車距離它必須有足夠的停車距離。少于此距離時(shí)不能停車；少于此距離時(shí)不能停車；大于此距離時(shí)必須停車；大于此距離時(shí)必須停車；等于此距離時(shí)可以停車，也可以通過路口。等于此距

24、離時(shí)可以停車，也可以通過路口。等于此距離時(shí)可以停車，也可以通過路口。當(dāng)決定通過等于此距離時(shí)可以停車，也可以通過路口。當(dāng)決定通過路口時(shí)，他必須有足夠的時(shí)間使他完全通過路口，這包括路口時(shí)，他必須有足夠的時(shí)間使他完全通過路口，這包括作出決定的時(shí)間、通過停車所需要的最短距離的駕駛時(shí)間作出決定的時(shí)間、通過停車所需要的最短距離的駕駛時(shí)間以及通過十字路口的時(shí)間。以及通過十字路口的時(shí)間。 于是，黃燈的狀態(tài)應(yīng)該持續(xù)的時(shí)間于是，黃燈的狀態(tài)應(yīng)該持續(xù)的時(shí)間包括駕駛員的決定包括駕駛員的決定時(shí)間（反應(yīng)時(shí)間）、停車距離的駕駛時(shí)間和他通過十字路時(shí)間（反應(yīng)時(shí)間）、停車距離的駕駛時(shí)間和他通過十字路口的時(shí)間?？诘臅r(shí)間。(2) 建模

25、與求解建模與求解 記記 T1 -駕駛員反應(yīng)時(shí)間駕駛員反應(yīng)時(shí)間 T2 -汽車通過十字路口的時(shí)間汽車通過十字路口的時(shí)間 T3 -停車距離的駕駛時(shí)間停車距離的駕駛時(shí)間則則T= T1 + T2 + T3 為黃燈應(yīng)該亮的時(shí)間。下面計(jì)算為黃燈應(yīng)該亮的時(shí)間。下面計(jì)算T2 ,T3 設(shè)法定行駛速度為設(shè)法定行駛速度為 v0，十字路口的長(zhǎng)度為，十字路口的長(zhǎng)度為I，典型的，典型的車身長(zhǎng)度為車身長(zhǎng)度為 L，則汽車通過十字路口的時(shí)間為，則汽車通過十字路口的時(shí)間為 注意車的尾部必須通過路口。注意車的尾部必須通過路口。02vLIT 停車過程是通過駕駛員踩動(dòng)剎車踏板產(chǎn)生一種摩停車過程是通過駕駛員踩動(dòng)剎車踏板產(chǎn)生一種摩擦力，使汽車減速直到停止。擦力，使汽車減速直到停止。設(shè)設(shè)m為汽車質(zhì)量，為汽車質(zhì)量，f為剎車摩擦系數(shù)，為剎車摩擦系數(shù)，x(t)為行駛距為行駛距離，剎車制動(dòng)力為離，剎車制動(dòng)力為fmg ,（g 為重力加速度為重力加速度）。）。由牛由牛頓第二定律，剎車過程滿足下述方程：頓第二定律，剎車過程滿足下述方程：T3 ：停車距離停車距離的駕駛時(shí)間的駕駛時(shí)間 0022,0)0(vdtdxxfmgdtxdmttvfgttx0221)( 0022,0)0(vdtdxxfmgdtxdmttvfgttx0221)( 0vxfg 1002TfgvvLIT T3 ：停車