matlab教學第二章微分方程

1、第二章第二章 微分方程模型微分方程模型2021-12-131本章學習目的：本章學習目的：l學習微分方程模型的建立、求解方法、分析結學習微分方程模型的建立、求解方法、分析結果及解決實際問題的全過程。果及解決實際問題的全過程。l熟練掌握使用熟練掌握使用MATLAB軟件的函數求微分方軟件的函數求微分方程的解析解、數值解和圖形解程的解析解、數值解和圖形解 。2021-12-1322.1 引例引例 l對于圓柱形狀容器壁上的容積刻度，對于圓柱形狀容器壁上的容積刻度，可以利用圓柱體體積公可以利用圓柱體體積公式：式： ，其中容器的直徑，其中容器的直徑D為常數，體積為常數，體積V與相對于容器底部的與相對于容器底

2、部的任意高度任意高度H成正比，因此在容器壁上成正比，因此在容器壁上可以方便地標出容積刻度?？梢苑奖愕貥顺鋈莘e刻度。l而對于幾何形狀不規(guī)則的容器，比如而對于幾何形狀不規(guī)則的容器，比如“倒葫蘆形狀倒葫蘆形狀”的容器壁上如何標出的容器壁上如何標出容積刻度呢？容積刻度呢？x4/2HDV2021-12-133l建立坐標系，由微元法分析可知：建立坐標系，由微元法分析可知：l l其中其中x表示高度，直徑是高度的函數，記為表示高度，直徑是高度的函數，記為D(x)l可得微分方程：可得微分方程：l如果該方程中的函數如果該方程中的函數D(x)無解析表達式，只給出無解析表達式，只給出D(x)的部分測試數據，如何求解此

3、微分方程呢？的部分測試數據，如何求解此微分方程呢？dxxDdV2)(410)0()(412VxDdxdV2021-12-134lh=0.2;ld=0.04,0.11,0.26,0.56,1.04,1.17;lx(1)=0;v(1)=0;lfor k=1:5lx(k+1)=x(k)+h;lv(k+1)=v(k)+(h/2)*(pi/4)*(d(k)2+d(k+1)2);lendlx=x(1:6),v=v(1:6),lplot(x,v) lx =l Columns 1 through 5 l 0 0.2000 0.4000 0.6000 0.8000l Column 6 l 1.0000lv =l

4、 Columns 1 through 5 l 0 0.0011 0.0073 0.0373 0.1469l Column 6 l 0.339300.10.20.30.40.50.60.70.80.9100.050.10.150.20.250.30.352021-12-1352.2 微分方程模型的建立微分方程模型的建立 l在工程實際問題中，在工程實際問題中，“改變改變”、“變化變化”、“增加增加”、“減少減少”等關鍵詞提示我們注意什么量在變化，關鍵詞等關鍵詞提示我們注意什么量在變化，關鍵詞“速率速率”、“增長增長”、“衰變衰變”、“邊際的邊際的”等常涉及到等常涉及到導數。導數。l我們熟悉的速度公

5、式：我們熟悉的速度公式： 就是一個簡單的一階微分方就是一個簡單的一階微分方程。程。l微分方程是指含有導數或微分的等式。微分方程是指含有導數或微分的等式。l一般形式：一般形式：l常用的建立微分方程的方法有：運用已知物理定律；利用常用的建立微分方程的方法有：運用已知物理定律；利用平衡與增長式；運用微元法；應用分析法。平衡與增長式；運用微元法；應用分析法。vdtdy).,(0),()1()()(nnnyyyxfyyyyxF或：2021-12-1362.2.1 運用已知物理定律運用已知物理定律l例例2.1 一個較熱的物體置于室溫為一個較熱的物體置于室溫為180C的房間內，的房間內，該物體最初的溫度是該

6、物體最初的溫度是600C，3分鐘以后降到分鐘以后降到500C。想。想知道它的溫度降到知道它的溫度降到300C 需要多少時間？需要多少時間？10分鐘以后分鐘以后它的溫度是多少？它的溫度是多少？l牛頓冷卻（加熱）定律牛頓冷卻（加熱）定律：將溫度為：將溫度為T的物體放入處于的物體放入處于常溫常溫 m 的介質中時，的介質中時，T的變化速率正比于的變化速率正比于T與周圍介與周圍介質的溫度差。質的溫度差。l分析：假設房間足夠大，放入溫度較低或較高的物體分析：假設房間足夠大，放入溫度較低或較高的物體時，室內溫度基本不受影響，即室溫分布均衡，保持時，室內溫度基本不受影響，即室溫分布均衡，保持為為m，采用牛頓冷

7、卻定律是一個相當好的近似。，采用牛頓冷卻定律是一個相當好的近似。2021-12-137l建立模型：設物體在冷卻過程中的溫度為建立模型：設物體在冷卻過程中的溫度為T(t)，t0， l根據牛頓加熱（冷卻）定律：根據牛頓加熱（冷卻）定律： ，建立微分方程建立微分方程 l其中參數其中參數k 0，m=18。）成正比與（mTdtdT60)0()(TmTkdtdT2021-12-1382.2.2 利用平衡與增長式利用平衡與增長式 l許多研究對象在數量上常常表現出某種不變的特性，許多研究對象在數量上常常表現出某種不變的特性，如封閉區(qū)域內的能量、貨幣量等。利用變量間的平衡如封閉區(qū)域內的能量、貨幣量等。利用變量間

8、的平衡與增長特性，可分析和建立有關變量間的相互關系。與增長特性，可分析和建立有關變量間的相互關系。l此類建模方法的關鍵是分析并正確描述基本模型的右此類建模方法的關鍵是分析并正確描述基本模型的右端，使平衡式成立。端，使平衡式成立。l例例2.2 戰(zhàn)斗模型：兩方軍隊交戰(zhàn)，希望為這場戰(zhàn)斗建戰(zhàn)斗模型：兩方軍隊交戰(zhàn)，希望為這場戰(zhàn)斗建立一個數學模型，應用這個模型達到如下目的：立一個數學模型，應用這個模型達到如下目的： l1. 預測哪一方將獲勝？預測哪一方將獲勝？ l2. 估計獲勝的一方最后剩下多少士兵？估計獲勝的一方最后剩下多少士兵？ l3. 計算失敗的一方開始時必須投入多少士兵才能贏得計算失敗的一方開始時

9、必須投入多少士兵才能贏得這場戰(zhàn)斗？這場戰(zhàn)斗？ 2021-12-139l解：模型建立：解：模型建立：l設設 x(t)： t 時刻時刻X方存活的士兵數方存活的士兵數l y(t)： t 時刻時刻Y方存活的士兵數方存活的士兵數l假設：假設：l1）雙方所有士兵不是戰(zhàn)死就是活著參加戰(zhàn)斗，）雙方所有士兵不是戰(zhàn)死就是活著參加戰(zhàn)斗，x(t)與與y(t)都是連續(xù)變量；都是連續(xù)變量； l2）Y方軍隊的一個士兵在單位時間內殺死方軍隊的一個士兵在單位時間內殺死X 方軍隊方軍隊 a 名名士兵；士兵； l3）X 方軍隊的一個士兵在單位時間內殺死方軍隊的一個士兵在單位時間內殺死Y方軍隊方軍隊 b 名名士兵；士兵；lt 時間內

10、時間內X軍隊減少的士兵數軍隊減少的士兵數 = t 時間內時間內Y軍隊消滅軍隊消滅對方的士兵數對方的士兵數 l即有即有 x =aytl同理同理 y =bxt 令0t， 得得 )0()0(bbxdtdyaaydtdx2021-12-13102.2.3 微元法微元法 l基本思想：通過分析研究對象的有關變量在一個基本思想：通過分析研究對象的有關變量在一個很短時間內的變化情況建立微分方程。很短時間內的變化情況建立微分方程。l本章引例本章引例2021-12-13112.2.4 分析法分析法l基本思想：根據對現實對象特性的認識，分析其因果基本思想：根據對現實對象特性的認識，分析其因果關系，找出反映內部機理的

11、規(guī)律。關系，找出反映內部機理的規(guī)律。l例例2.4 獨家廣告模型獨家廣告模型 廣告是調整商品銷售的強有力廣告是調整商品銷售的強有力的手段，廣告與銷售量之間有什么內在聯系？如何評的手段，廣告與銷售量之間有什么內在聯系？如何評價不同時期的廣告效果？價不同時期的廣告效果？l解：解：1分析廣告的效果，可做如下的條件假設：分析廣告的效果，可做如下的條件假設：l商品的銷售速度會因廣告而增大，當商品在市場上商品的銷售速度會因廣告而增大，當商品在市場上趨于飽和時，銷售速度將趨于一個極限值；趨于飽和時，銷售速度將趨于一個極限值；l商品銷售率（銷售加速度）隨商品銷售速度的增高商品銷售率（銷售加速度）隨商品銷售速度的

12、增高而降低。而降低。 2021-12-1312l2符號說明符號說明lA(t) t 時刻的廣告費用時刻的廣告費用lS(t) t 時刻商品的銷售速度；時刻商品的銷售速度； lM 銷售飽和水平，即銷售速度的上限；銷售飽和水平，即銷售速度的上限； l 衰減因子，廣告作用隨時間的推移而自然衰減因子，廣告作用隨時間的推移而自然衰減的速度，衰減的速度，0；lp 響應系數，表征響應系數，表征A(t) 對對 S(t) 的影響力。的影響力。2021-12-1313l3模型建立模型建立l選擇如下廣告策略，選擇如下廣告策略，t時刻的廣告費用為：時刻的廣告費用為：l建立微分方程：建立微分方程：l模型分析：是否與假設相符

13、？模型分析：是否與假設相符？ttAtA, 00,)()()(1)(tSMtStpAdtdS2021-12-13142.3 微分方程求解方法微分方程求解方法 l2.3.1 解析解法解析解法l解析解法只能解決一些特殊微分方程，這些方解析解法只能解決一些特殊微分方程，這些方法主要針對：法主要針對：l一階特殊的微分方程：如使用分離變量法、方一階特殊的微分方程：如使用分離變量法、方程變換法、線性方程的常數變易法或公式法求程變換法、線性方程的常數變易法或公式法求解。解。l二階或高階常系數線性微分方程的特征根法。二階或高階常系數線性微分方程的特征根法。l在高等數學的教程中有專門介紹。在高等數學的教程中有專門

14、介紹。l下面著重介紹微分方程的數值解法。下面著重介紹微分方程的數值解法。2021-12-13152.3.2 數值解法數值解法l微分方程的數值解法是解決某些實際問題中經常使用的方微分方程的數值解法是解決某些實際問題中經常使用的方法。設待求解的定解問題為法。設待求解的定解問題為l求該問題數值解法的基本過程如下：求該問題數值解法的基本過程如下：l引入自變量取值點序列引入自變量取值點序列 ，定義，定義 為步長，常為步長，常用定步長用定步長( 與與n無關，為常數無關，為常數)，其精確解記為，其精確解記為 ，一般難以得到。為了尋求一般難以得到。為了尋求 的近似值的近似值 ，設想根據一，設想根據一定的原理，

15、結合當前得到近似解，近似地表示該點或前一定的原理，結合當前得到近似解，近似地表示該點或前一點的導數值，由此推出計算點的導數值，由此推出計算 的迭代公式。因此數值解法的迭代公式。因此數值解法一般只能得到微分方程的近似解一般只能得到微分方程的近似解 。下面介紹兩個微分。下面介紹兩個微分方程中最常用的數值解法。方程中最常用的數值解法。00)(),(yxyyxfdxdy nx1nnnxxhnh)(nxy)(nxynyny ny2021-12-13161.歐拉方法歐拉方法l這是一種最簡單的解微分方程的數值方法：就是在小區(qū)間這是一種最簡單的解微分方程的數值方法：就是在小區(qū)間xn, xn+1上用差商代替微商

16、，可以得到近似的表達式上用差商代替微商，可以得到近似的表達式l若若f(x,y)中的中的x取左端點取左端點 ，結合已經得到的，結合已經得到的y(xn)的近似的近似值（數值解）值（數值解）yn，即，即 ，有，有y(xn+1)的近似值為的近似值為l n = 0,1, l這就是求解微分方程的這就是求解微分方程的顯式歐拉公式顯式歐拉公式。也稱。也稱向前歐拉公式向前歐拉公式。l向前歐拉法計算簡單，易于計算，但精度不高，收斂速度向前歐拉法計算簡單，易于計算，但精度不高，收斂速度慢慢),()()(1yxfhxyxynnnnyxy)(),(1nnnnyxfhyynx2021-12-1317l若若f(x,y)中的

17、中的x取右端點，可得取右端點，可得向后歐拉公式向后歐拉公式如下：如下：lyn+1 = yn + h f (xn +1, yn +1) ，n = 0,1,l稱為隱式公式，因為要得出數值解稱為隱式公式，因為要得出數值解yn+1，就必須求解這個，就必須求解這個非線性方程，計算比較困難。非線性方程，計算比較困難。l如果用將向前和向后歐拉公式加以平均，可得到如果用將向前和向后歐拉公式加以平均，可得到梯形公式梯形公式：l該法的計算精度比向前和向后歐拉法都高，但計算和向后該法的計算精度比向前和向后歐拉法都高，但計算和向后歐拉法一樣困難。歐拉法一樣困難。),(),(21111nnnnnnyxfyxfhyy20

18、21-12-1318l改進的歐拉算法：改進的歐拉算法：l（1）先用向前歐拉法算出）先用向前歐拉法算出yn+1的預測值的預測值 ，ll（2）將預測值代入梯形公式的右端作為校正，得到）將預測值代入梯形公式的右端作為校正，得到yn+1l ，n=1,2,l該式稱為該式稱為改進歐拉公式改進歐拉公式。1ny),(1nnnnyxfhyy),(),(21111nnnnnnyxfyxfhyy2021-12-1319l例例2.5 求解微分方程求解微分方程ly = -y +x +1, y(0) = 1, 取步長取步長h = 0.1和和0.001。l分別用三種數值解法求解，并結合其精確解，對求解誤差進分別用三種數值解

19、法求解，并結合其精確解，對求解誤差進行分析比較。行分析比較。l解解 這是一個一階線性微分方程，可用解析解法得到其精確解這是一個一階線性微分方程，可用解析解法得到其精確解y = x + e-x。2021-12-1320l三種數值解如下：三種數值解如下：l 1) 向前歐拉法：迭代公式為向前歐拉法：迭代公式為 yn+1 = （1-h）yn + hxn + h，n=0,1,.。其中。其中y0= y(0) = 1。l2) 后退歐拉法：由后退歐拉法隱式公式得后退歐拉法：由后退歐拉法隱式公式得lyn+1 = yn + h(-yn+1+xn+1+1)，變形為，變形為lyn+1 = (yn + hxn+1 +

20、h)/(1+h)。l3) 梯形法：將隱式梯形公式轉化為顯示迭代公式梯形法：將隱式梯形公式轉化為顯示迭代公式如下：如下：lyn+1 = (yn + (h/2)*(- yn + xn + xn+1 +2)/(1+h/2)。2021-12-1321lx1(1)=0;y1(1)=1;y2(1)=1;y3(1)=1;h=0.1;lfor k=1:10lx1(k+1)=x1(k)+h;ly1(k+1)=(1-h)*y1(k)+h*x1(k)+h;ly2(k+1)=(y2(k)+h*x1(k+1)+h)/(1+h);ly3(k+1)=(y3(k)+(h/2)*(-y3(k)+x1(k)+x1(k+1)+2)

21、/(1+h/2);lendlx=0:0.1:1;ly=x+exp(-x);lx1=x1(1:11),y=y(1:11),y1=y1(1:11),y2=y2(1:11),y3=y3(1:11),lplot(x,y,x1,y1,k:,x1,y2,r-,x1,y3,g*) l程序中，程序中，x1為自變量，為自變量，y為精確解，為精確解，y1、y2、y3分別為向分別為向前歐拉法、后退歐拉法和梯形法的解。結果如下前歐拉法、后退歐拉法和梯形法的解。結果如下：2021-12-1322表表2-1當當h = 0.1時時nx精確解精確解向前歐拉法向前歐拉法后退歐拉法后退歐拉法梯形法梯形法011110.11.004

22、811.00911.00480.21.01871.01001.02641.01860.31.04081.02901.05131.04060.41.07031.05611.08301.07010.51.10651.09051.12091.10630.61.14881.13141.16451.14850.71.19661.17831.21321.19630.81.24931.23051.26651.24900.91.30661.28741.32411.306311.36791.34871.38551.36762021-12-132300.10.20.30.40.50.60.70.80.9111.0

23、51.11.151.21.251.31.351.4圖中，藍色曲線是精確解，黑色曲線是向前歐拉法曲線，圖中，藍色曲線是精確解，黑色曲線是向前歐拉法曲線，紅色曲線是向后歐拉法曲線，綠色紅色曲線是向后歐拉法曲線，綠色“* ”號為梯形法曲線。號為梯形法曲線。24表表2-2 當當h = 0.001時時 nx精確解精確解向前歐拉法向前歐拉法后退歐拉法后退歐拉法梯形法梯形法011110.11.00481.00481.00491.00480.21.01871.01861.01881.01870.31.04081.04071.04091.04080.41.07031.07021.07051.07030.51.1

24、0651.10641.10671.10650.61.14881.14861.14901.14880.71.19661.19641.19681.19660.81.24931.24911.24951.24930.91.30661.30641.30681.306611.36791.36771.36811.36792021-12-132500.20.40.60.811.21.411.051.11.151.21.251.31.351.4 計算結果表明：當步長計算結果表明：當步長h =0.1時，它們的前兩位有效數字時，它們的前兩位有效數字是精確的；當步長是精確的；當步長h =0.001時，它們的前四位有效

25、數字是精確時，它們的前四位有效數字是精確的。說明在迭代中，步長的。說明在迭代中，步長h越小，計算結果越精確。越小，計算結果越精確。 2021-12-1326l衡量求解公式好壞的一個主要標準是求解公式的精度，因衡量求解公式好壞的一個主要標準是求解公式的精度，因此引入局部此引入局部截斷誤差和階數截斷誤差和階數的概念。的概念。l定義定義2.1 假定假定 沒有誤差，即沒有誤差，即 ，用數值方法計，用數值方法計算算 的誤差的誤差 ，稱為該數值方法計算時的，稱為該數值方法計算時的局部截斷誤差。局部截斷誤差。l為估計歐拉公式的局部截斷誤差，先將精確解為估計歐拉公式的局部截斷誤差，先將精確解 在在 處作泰勒級

26、數展開處作泰勒級數展開l對于向前歐拉公式，在對于向前歐拉公式，在 的假定下可記作的假定下可記作l式中式中h的最低階項稱為局部截斷誤差主項，它對于的最低階項稱為局部截斷誤差主項，它對于h是是2階階的。的。ny)(nnxyy 1ny11)(nnyxy)(1nxynx)()(! 2)()()(321hOxyhxyhxyxynnnn )()()(,()(1nnnnnnxyhxyxyxhfxyy)(nnxyy )()()(2)(232111hOhOxyhyxynnnn 2021-12-1327l同理，對于向后歐拉公式的局部截斷誤差同理，對于向后歐拉公式的局部截斷誤差是是 ，梯形公式和改進歐拉公式的局部，

27、梯形公式和改進歐拉公式的局部截斷誤差是截斷誤差是 。l求解微分方程的數值計算方法的精度是由局部截求解微分方程的數值計算方法的精度是由局部截斷誤差中斷誤差中h的階定義的：如果一個方法的局部截的階定義的：如果一個方法的局部截斷誤差為斷誤差為 ，則稱該方法具有，則稱該方法具有p階精度。階精度。l因此向前和向后歐拉方法的精度為因此向前和向后歐拉方法的精度為1階，梯形公階，梯形公式和改進歐拉公式的精度為式和改進歐拉公式的精度為2階。階。)(21hOn)(31hOn)(1phO2021-12-1328l數值方法具有越高階的精度，得到的解就越精確。數值方法具有越高階的精度，得到的解就越精確。我們發(fā)現，向前和

28、向后歐拉公式各用了區(qū)間一個我們發(fā)現，向前和向后歐拉公式各用了區(qū)間一個端點的導數，具有端點的導數，具有1階精度，而梯形和改進歐拉階精度，而梯形和改進歐拉公式用了兩個端點的導數取平均，得到了公式用了兩個端點的導數取平均，得到了2階精階精度。因此就引導人們利用內的若干點的導數，對度。因此就引導人們利用內的若干點的導數，對它們作線性組合得到平均斜率，由此得到更高階它們作線性組合得到平均斜率，由此得到更高階的精度，這就是的精度，這就是龍格龍格-庫塔方法庫塔方法的基本思路。的基本思路。l在在MATLAB軟件中含有數值求解的系統函數，其軟件中含有數值求解的系統函數，其實現原理就是龍格庫塔方法。實現原理就是龍

29、格庫塔方法。 2021-12-13292.4 Matlab軟件求解軟件求解 l2.4.1 解析解解析解l用MATLAB命令ldsolve(eqn1,eqn2, .) 求常微分方程（組）的解析解。l其中eqni表示第i個微分方程，Dny表示y的n階導數,默認的自變量為t。2021-12-1330l1. 微分方程微分方程l例例2.6 求解一階微分方程 l(1) 求通解l輸入：ldsolve(Dy=1+y2) l輸出：lans =ltan(t+C1) 21ydxdy（2）求特解輸入：dsolve(Dy=1+y2,y(0)=1,x) 指定初值為1，自變量為x輸出：ans =tan(x+1/4*pi)

30、2021-12-1331l例例2.7 求解二階微分方程求解二階微分方程 l輸入輸入:ldsolve(D2y+(1/x)*Dy+(1-(1/2)2/x2)*y=0,y(pi/2)=2,Dy(pi/2)=-2/pi,x) lans =2(1/2)*pi(1/2)/x(1/2)*sin(x) l化簡輸出結果，輸入：化簡輸出結果，輸入：lpretty(ans) 1/2 1/2l 2 pi sin(x)l -l 1/2l x2/1/2)2/(2)2/(0)(222 nyyynxyxyxxxysin2即：即： 2021-12-13322微分方程組微分方程組l例2.8 求解 df/dx=3f+4g; dg/

31、dx=-4f+3g。l（1）通解: lf,g=dsolve(Df=3*f+4*g,Dg=-4*f+3*g) lf =lexp(3*t)*(C1*sin(4*t)+C2*cos(4*t)lg =lexp(3*t)*(C1*cos(4*t)-C2*sin(4*t) l特解：lf,g=dsolve(Df=3*f+4*g,Dg=-4*f+3*g,f(0)=0,g(0)=1) lf =lexp(3*t)*sin(4*t)lg =lexp(3*t)*cos(4*t) 2021-12-13332.4.2 數值解數值解l在微分方程在微分方程(組組)難以獲得解析解的情況下，可以用難以獲得解析解的情況下，可以用M

32、atlab方便地求出數值解。格式為方便地求出數值解。格式為:lt,y = ode23(F,ts,y0)， t,y = ode45(F,ts,y0)l注意：注意：l微分方程的形式：微分方程的形式：y = F(t, y)，t為自變量，為自變量，y為因變?yōu)橐蜃兞浚梢允嵌鄠€，如微分方程組）；量（可以是多個，如微分方程組）；lt, y為輸出矩陣，分別表示自變量和因變量的取值；為輸出矩陣，分別表示自變量和因變量的取值；lF代表微分方程組的函數名（代表微分方程組的函數名（m文件，必須返回一個文件，必須返回一個列向量）；列向量）；lts的取法有幾種，（的取法有幾種，（1）ts=t0, tf 表示自變量的取值

33、表示自變量的取值范圍，（范圍，（2）ts=t0,t1,t2,tf,則輸出在指定時刻則輸出在指定時刻t0,t1,t2,tf處給出，（處給出，（3）ts=t0:k:tf,則輸出在區(qū)間則輸出在區(qū)間t0,tf的等分點給出；的等分點給出；ly0為初值條件。為初值條件。2021-12-1334l比如例例2.7的數值解：l解：令y1=y,y2=y1,將二階微分方程轉化為一階微分方程組ly1=y2ly2=-y2/x+(n/x)2-1)y1l首先建立M-文件函數:l function f=jie3(x,y)l f=y(2);-y(2)/x+(1/2)2/x2-1)*y(1);l計算:lx,y=ode23(jie

34、3,pi/2,pi,2,-2/pi) 2/1/2)2/(2)2/(0)(222 nyyynxyxyx2021-12-1335lx =l 1.5708l 1.6074l 1.7645l 1.9215l 2.0786l 2.2357l 2.3928l 2.5499l 2.7069l 2.8640l 3.0211l 3.1416ly =l 2.0000 -0.6366l 1.9758 -0.6869l 1.8518 -0.8879l 1.6982 -1.0631l 1.5192 -1.2108l 1.3193 -1.3293l 1.1032 -1.4174l 0.8756 -1.4744l 0.64

35、16 -1.5002l 0.4060 -1.4951l 0.1735 -1.4602l 0.0002 -1.4140 2021-12-1336y1=y(:,1);y2=y(:,2);plot(x,y1,x,y2,r)2021-12-13371.41.61.822.22.42.62.833.2-2-1.5-1-0.500.511.522.5 范例：地中海鯊魚問題范例：地中海鯊魚問題 l意大利生物學家意大利生物學家Ancona曾致力于魚類種群相互制約關曾致力于魚類種群相互制約關系的研究，他從第一次世界大戰(zhàn)期間（系的研究，他從第一次世界大戰(zhàn)期間（1914年年7月月28日日1918年年11月月11日）

36、日）,地中海各港口捕獲的幾種魚類地中海各港口捕獲的幾種魚類捕獲量百分比的資料中，發(fā)現鯊魚的比例有明顯增加捕獲量百分比的資料中，發(fā)現鯊魚的比例有明顯增加（見下表）。（見下表）。年代年代1914191419151915191619161917191719181918百分比百分比11.911.921.421.422.122.121.221.236.436.4年代年代1919191919201920192119211922192219231923百分比百分比27.327.316.016.015.915.914.814.819.719.72021-12-1338l戰(zhàn)爭為什么使鯊魚數量增加？是什么原因？戰(zhàn)

37、爭為什么使鯊魚數量增加？是什么原因？l因為戰(zhàn)爭使捕魚量下降，食用魚增加，顯然鯊魚也隨之增因為戰(zhàn)爭使捕魚量下降，食用魚增加，顯然鯊魚也隨之增加。加。 l但為何鯊魚的比例大幅增加呢？生物學家但為何鯊魚的比例大幅增加呢？生物學家Ancona無法解無法解釋這個現象，于是求助于著名的意大利數學家釋這個現象，于是求助于著名的意大利數學家V.Volterra，希望建立一個食餌希望建立一個食餌捕食者系統的數學模型，定量地回答捕食者系統的數學模型，定量地回答這個問題。這個問題。l 1、符號說明：、符號說明：lx1(t), x2(t)分別是食餌、捕食者（鯊魚）在分別是食餌、捕食者（鯊魚）在t時刻的數時刻的數量；量

38、； lr1, r2是食餌、捕食者的固有增長率；是食餌、捕食者的固有增長率；l1是捕食者掠取食餌的能力，是捕食者掠取食餌的能力， l 2是食餌對捕食者的供養(yǎng)能力；是食餌對捕食者的供養(yǎng)能力；2021-12-1339l2、基本假設：、基本假設：l 捕食者的存在使食餌的增長率降低，假設捕食者的存在使食餌的增長率降低，假設x1降低的降低的程度與捕食者數量程度與捕食者數量x2成正比，即成正比，即l食餌對捕食者的數量食餌對捕食者的數量x2起到增長的作用，起到增長的作用， 其程度其程度與食餌數量與食餌數量x1成正比，即成正比，即l綜合以上綜合以上和和，得到如下模型：，得到如下模型：l模型一：不考慮人工捕獲的情

39、況模型一：不考慮人工捕獲的情況)(21111xrxdtdx)(12222xrxdtdx)()(1222221111xrxdtdxxrxdtdx2021-12-1340l該模型反映了在沒有人工捕獲的自然環(huán)境中食餌與捕食者該模型反映了在沒有人工捕獲的自然環(huán)境中食餌與捕食者之間的制約關系，沒有考慮食餌和捕食者自身的阻滯作用，之間的制約關系，沒有考慮食餌和捕食者自身的阻滯作用，是是Volterra提出的最簡單的模型。提出的最簡單的模型。l給定一組具體數據，用給定一組具體數據，用matlab軟件求解。軟件求解。 l食餌：食餌： r1= 1, 1= 0.1, x10= 25; l捕食者捕食者(鯊魚鯊魚)：r2=0.5, 2=0.02, x20= 2;l編制程序如下編制程序如下l1、建立、建立m-文件文件shier.m如下：如下：l function dx=shier(t,x) l dx=zeros(2,1); l dx(1)=x(1)*(1-0.1*x(2); l dx(2)=x(2)*(-0.5+0.02*x(1);2021-12-1341l2、在命令窗口執(zhí)行如下程序：、在命令窗口執(zhí)行如下程序： l t,x=ode4