常微分方程數值解

1、大學數學實驗大學數學實驗Mathematical Experiments 實驗實驗4 常微分方程數值解常微分方程數值解為什么要學習微分方程數值解為什么要學習微分方程數值解 微分方程是研究函數變化規(guī)律的重要工具，有著廣泛微分方程是研究函數變化規(guī)律的重要工具，有著廣泛的應用。如的應用。如: 物體的運動物體的運動, 電路的電壓電路的電壓, 人口增長的預測人口增長的預測 許多微分方程沒有解析解，數值解法是求解的重要手許多微分方程沒有解析解，數值解法是求解的重要手段，如段，如2,dyyxdx( )( )x taxyy taxyby 實驗實驗4 4的基本內容的基本內容3. 實際問題用微分方程建模，并求解實

2、際問題用微分方程建模，并求解2. 龍格龍格-庫塔方法的庫塔方法的MATLAB實現實現*4. 數值算法的收斂性、穩(wěn)定性與剛性方程數值算法的收斂性、穩(wěn)定性與剛性方程1. 兩個最常用的數值算法兩個最常用的數值算法: 歐拉（歐拉（Euler）方法）方法 龍格龍格-庫塔（庫塔（Runge-Kutta）方法）方法實例實例1 1 海上緝私海上緝私 海防某部緝私艇上的雷達發(fā)現正東方向海防某部緝私艇上的雷達發(fā)現正東方向c海里處有一艘走私船海里處有一艘走私船正以速度正以速度a向正北方向行駛，緝私艇立即以最大速度向正北方向行駛，緝私艇立即以最大速度b(a)前往前往攔截。如果用雷達進行跟蹤時，可保持緝私艇的速度方向始

3、終攔截。如果用雷達進行跟蹤時，可保持緝私艇的速度方向始終指向走私船。指向走私船。 建立任意時刻緝私艇位置及建立任意時刻緝私艇位置及 航線的數學模型航線的數學模型,并求解并求解; 求出緝私艇追上走私船的時間。求出緝私艇追上走私船的時間。 a北北bc艇艇船船實例實例1 1 海上緝私海上緝私 建立坐標系如圖建立坐標系如圖: t=0 艇在艇在(0, 0), 船在船在(c, 0); 船速船速a, 艇速艇速b 時刻時刻 t 艇位于艇位于P(x, y), 船到達船到達 Q(c, at)模型模型:0yxcR(c,y ) Q(c,at)P(x,y)bcos ,sindxdybbdtdt2222()()()()(

4、)()dxb cxdtcxatydyb atydtcxaty由方程無法得到由方程無法得到x(t), y(t)的解析解的解析解需要用數值解法求解需要用數值解法求解 “常微分方程初值問題數值解常微分方程初值問題數值解”的提法的提法00( , ), ()yf x yy xy設的解y=y(x)存在且唯一不求解析解不求解析解 y = y(x) (無解析解或求解困難無解析解或求解困難)12nxxx而在一系列離散點而在一系列離散點()(1, 2 ,)nyxn n求的 近 似 值 y通常取等步長通常取等步長h0nxxnh0 x0y)(xyy yxy1y2yn)(nxy1x2xnxyP0 x0 x1x2x3 x

5、y=y(x)y0P1P2P3歐 拉 方 法00( , ), ()yf x y y xy基本基本思路思路1,nnx x在小區(qū)間上1( , ),nnf x yx x中的x取內的某一點1 ()()/ ,nnyy xy xh11()( )( , ( ), ,nnnny xy xhf x y xxx xx取不同點取不同點各種歐拉公式各種歐拉公式nxx取左端點1()()( , ()nnnny xy xhf x y x1( ,),0,1,nnnnyyhf x yn向前歐拉公式向前歐拉公式顯式公式顯式公式11(),()nnnnyy xyy xP3歐 拉 方 法11()()( , ( ),nnnny xy xh

6、f x y xxx x11(),()nnnnxyy xyy x取右端點，111(,),0,1,nnnnyyhf xyn向后歐拉公式向后歐拉公式 隱式公式隱式公式y(tǒng)P0 x0 x1x2x3 xy0y=y(x)P1P2n+1y右端未知，需迭代求解(0)1(,)nnnnyyhfxy初值(1)()111(,)0,1,2,0,1,2,kknnnnyyhfxykn()11limknnkyy 向前歐拉公式向前歐拉公式向后歐拉公式向后歐拉公式1(,)nnnnyyhf xy111(,)nnnnyyhfxy二者平均得到二者平均得到梯形公式梯形公式111(,)(,),0,1,2nnnnnnhyyfxyfxyn仍為隱

7、式公式，需迭代求解仍為隱式公式，需迭代求解將梯形公式的迭代過程簡化為兩步將梯形公式的迭代過程簡化為兩步1(,)nnnnyyhf x y預測預測111(,)(,),0,1,2nnnnnnhyyfxyfxyn校正校正改進歐拉公式改進歐拉公式龍格龍格- -庫塔方法庫塔方法11()()( , ( ),nnnny xy xhf x y xxx x 向前，向后歐拉公式：向前，向后歐拉公式： 用用xn, xn+1內個點的導內個點的導 數代替數代替 f(x, y(x) 梯形公式，改進歐拉公式：梯形公式，改進歐拉公式：用用xn, xn+1內個點導數的平均值代替內個點導數的平均值代替 f(x, y(x)龍格龍格-

8、 -庫塔方法的基本思想庫塔方法的基本思想在在xn, xn+1內多取幾個點，將它們的導數加權平內多取幾個點，將它們的導數加權平均代替均代替 f(x, y(x)，設法構造出設法構造出精度更高精度更高的計算公式。的計算公式。 1, 10, 1,111ijijiLiiijiiacac滿足常用的常用的(經典經典)龍格龍格庫塔庫塔公式公式不足不足:收斂速度較慢收斂速度較慢111222111(,)(,)(,),3,4,Lnniiinnnniininiijjjyyhkkf xykf xc h yc hkkf xc h yc ha kiLL級?階龍格龍格-庫塔庫塔方法的方法的一般形式一般形式1123412132

9、43(22)6(,)(/2,/2)(/2,/2)(,)nnnnnnnnnnhyykkkkkf xykf xhyhkkf xhyhkkf xh yhk使精度盡量高使精度盡量高4級4階微分方程組和高階方程初值問題的數值解微分方程組和高階方程初值問題的數值解 歐拉方法和龍格歐拉方法和龍格-庫塔方法可直接推廣到微分方程組庫塔方法可直接推廣到微分方程組0000( , , )( , , )(), ()yf x y zzg x y zy xyz xz向前歐拉公式向前歐拉公式 ),(yyxfy 1yy),(21221yyxfyyy高階方程需要先降階為一階微分方程組高階方程需要先降階為一階微分方程組 , 2 ,

10、 1 , 0),(),(11nzyxhgzzzyxhfyynnnnnnnnnn龍格龍格庫塔方法的庫塔方法的 MATLAB 實現實現 TnTnfffxxxxtxxtftx),(,),(,)(),()(1100t,x=ode23(f,ts,x0,opt) 3級級2階龍格階龍格-庫塔公式庫塔公式 t,x=ode45(f,ts,x0,opt) 5級級4階龍格階龍格-庫塔公式庫塔公式 f是待解方程寫成是待解方程寫成的函數的函數m文件：文件： function dx=f(t,x)dx=f1; f2; fn;ts = t0,t1, ,tf輸出指定時刻輸出指定時刻 t0,t1, ,tf 的函數值的函數值ts

11、= t0:k:tf輸出輸出 t0,tf 內等分點處的函數值內等分點處的函數值x0為函數初值為函數初值(n維維) 輸出輸出t=ts, x為相應函數值為相應函數值(n維維) optopt為選項，缺省時精度為：相對誤差為選項，缺省時精度為：相對誤差1010-3-3，絕對誤差絕對誤差1010-6-6, , 計算步長按精度要求自動調整計算步長按精度要求自動調整實例實例1海上緝私海上緝私( (續(xù)續(xù)) ) 模型的數值解模型的數值解2222()()()()()()dxb cxdtcxatydyb atydtcxaty(0)0,(0)0 xy0yxc(x(t),y(t)ab設設: 船速船速a=20 (海里海里/

12、小時小時) 艇速艇速b=40 (海里海里/小時小時) 距離距離c=15 (海里海里) 求求: 緝私艇的位置緝私艇的位置x(t),y(t) 緝私艇的航線緝私艇的航線y(x) jisi.m, seajisi.m 實例實例1海上緝私海上緝私( (續(xù)續(xù)) ) 模型的數值解模型的數值解 a=20, b=40, c=15走私船的位置走私船的位置x1(t)= c=15y1(t)=at=20t t=0.5時緝私艇追上走私船時緝私艇追上走私船 051015200246810 xy緝私艇的航線緝私艇的航線y(x)tx(t)y(t)0000.051.99840.06980.103.98540.29240.155.9

13、4450.69060.207.85151.28990.259.67052.11780.3011.34963.20050.3512.81704.55520.4013.98066.17730.4514.74518.02730.5015.00469.9979y1(t)01.02.03.04.05.06.07.08.09.010.0實例實例1 海上緝私海上緝私( (續(xù)續(xù)) ) 模型的數值解模型的數值解 設設b，c不變，不變，a變大變大為為30, 35, 接近接近40, 觀察解的變化觀察解的變化 :a=35, b=40, c=15t=? 緝私艇追上走私船緝私艇追上走私船 t x(t) y(t) y1(t

14、)0 0 0 00.13.95610.50583.50.27.59282.13087.00.310.52404.828310.50.412.53848.275514.00.513.755112.083017.51.214.998640.016442.01.314.999644.016545.51.415.011748.018349.01.515.002352.014652.51.614.986655.948656.0累積誤差較大累積誤差較大提高精度提高精度! ! 實例實例1海上緝私海上緝私( (續(xù)續(xù)) ) 模型的數值解模型的數值解 a=35, b=40, c=15opt=odeset(RelT

15、ol,1e-6, AbsTol,1e-9);t,x=ode45(jisi,ts,x0,opt);t=1.6時緝私艇追上走私船時緝私艇追上走私船 051015200204060 xy緝私艇的航線緝私艇的航線y(x)判斷判斷“追上追上”的有效方的有效方法法?t x(t) y(t)y1(t)0 0 0 00.13.9561040.5058133.50.27.5928222.1306787.00.310.5219214.82930810.50.412.5394548.26984014.00.513.75397412.07534417.51.214.99961640.00000542.01.314.99

16、996344.00000545.51.414.99999348.00000549.01.514.99999852.00000552.51.615.00002055.99993156.0實例實例1 海上緝私海上緝私( (續(xù)續(xù)) ) 模型的解析解模型的解析解 2222)()()(,)()()(yatxcyatbdtdyyatxcxcbdtdx)(xcyatdxdyatydxdyxc )(22()d ydtcxadxdxdsbdt22()()dsdxdy222()1 ()d yadycxdxbdxdxdyp 令bakpkdxdpxc,1)(2(0)0pdydsdsdtdxdtdxdtdtdydxdy

17、實例實例1海上緝私海上緝私( (續(xù)續(xù)) ) 模型的解析解模型的解析解 21)(pkdxdpxc(0)0p21()kcxppc21()kcxppc1()() 2kkcxcxpcc/1ka b(0)0y11211()()2 111kkccxcxkcykckck緝私艇的航線緝私艇的航線y y( (x x) )的解析解的解析解x=c時時 緝私艇追上走私船的緝私艇追上走私船的y y坐標坐標 緝私艇追上走私船的時間緝私艇追上走私船的時間: 122bctbaa=20, b=40, c=15 t1=0.5 a=35, b=40, c=15 t1=1.6 2221kcabcykba微分方程數值解法的誤差分微分方

18、程數值解法的誤差分析析數值解法數值解法: 計算微分方程精確解計算微分方程精確解y(xn)的近似值的近似值yn按照步長按照步長h一步步計算一步步計算, 每步都有誤差每步都有誤差;每一步的誤差會逐步積累每一步的誤差會逐步積累, , 稱稱累積誤差累積誤差. .討論計算一步出現的誤差討論計算一步出現的誤差假定假定yn= y(xn) , 估計估計yn+1的誤差的誤差: y(xn+1)- yn+1局部截斷誤差局部截斷誤差 誤差分析誤差分析估計歐拉公式的局部截斷誤差估計歐拉公式的局部截斷誤差 y(xn+1)在在xn處作處作Taylor展開展開: )()(2)()()(321hOxyhxyhxyxynnnn

19、向前歐拉公式向前歐拉公式1(,)nnnnyyhf xyyn= y(xn)()()(,()(1nnnnnnxyhxyxyxhfxyy232111()()()()2nnnnhTy xyyxO hO h局部截斷誤差主項為局部截斷誤差主項為誤差分析誤差分析估計歐拉公式的局部截斷誤差估計歐拉公式的局部截斷誤差 )()(2)()()(321hOxyhxyhxyxynnnn 向后歐拉公式向后歐拉公式),(111nnnnyxhfyy232111()()()()2nnnnhTy xyyxO hO h 局部截斷誤差主項為局部截斷誤差主項為)(22nxyh 梯形梯形公式公式向前、向后歐拉公式的平均向前、向后歐拉公式

20、的平均3431()()()12nnhTyxO hO h xn xn+1 x y Pn A Q B誤誤差差分分析析算法算法精度的階精度的階的定義的定義 一個算法的局部截斷誤差為一個算法的局部截斷誤差為O(hp+1) 該算法具有該算法具有p階精度階精度 局部截斷誤差局部截斷誤差精度精度向前歐拉公式向前歐拉公式O(h2)1階階向后歐拉公式向后歐拉公式O(h2)1階階梯形公式梯形公式O(h3)2階階改進歐拉公式改進歐拉公式O(h3)2階階經典龍格經典龍格- -庫塔公式庫塔公式O(h5)4階階實實 例例 2弱弱 肉肉 強強 食食 問問題題自然界中同一環(huán)境下兩個種群之間的生存方式自然界中同一環(huán)境下兩個種群

21、之間的生存方式 相互競爭相互競爭 相互依存相互依存 弱肉強食弱肉強食 弱肉弱肉強食強食種群甲靠豐富的自然資源生存種群甲靠豐富的自然資源生存 食餌食餌(Prey) (Prey) 種群乙靠捕食種群甲為生種群乙靠捕食種群甲為生 捕食者捕食者(Predator) (Predator) 兩個種群的數量如何演變兩個種群的數量如何演變? ? 歷史背景歷史背景一次世界大戰(zhàn)期間地中海漁業(yè)一次世界大戰(zhàn)期間地中海漁業(yè)的捕撈量下降的捕撈量下降(食用魚和鯊魚同時捕撈食用魚和鯊魚同時捕撈)，但，但是其中是其中鯊魚的比例卻增加，為什么？鯊魚的比例卻增加，為什么？實實 例例 2弱弱 肉肉 強強 食食 模型模型食餌食餌(甲甲)

22、 的密度的密度x(t), 捕食者捕食者(乙乙)的密度的密度y(t)0,/rrxx 甲獨立生存的增長率甲獨立生存的增長率r r0,/aayrxx 乙使甲的增長率減小乙使甲的增長率減小, ,減小量與減小量與 y y 成正比成正比0,/ddyy 乙獨立生存的死亡率乙獨立生存的死亡率d0),(/bbxdyy 甲使乙的死亡率減小甲使乙的死亡率減小, ,減小量與減小量與 x成正比成正比00()()(0),(0)xray xrxaxyydbx ydybxyxxyyVolterra模型模型 x(t), y(t)無解析解無解析解實例實例2弱肉強食弱肉強食 模型的數值解模型的數值解 00)0(,)0(yyxxbx

23、ydyyaxyrxx001,0.50.1,0.0225,2rdabxy051015020406080100 x(t)y(t)020406080100051015202530 xy猜測猜測 x(t), y(t)是周期函數是周期函數; y(x)是封閉曲線是封閉曲線 數值積分計算一個周期的平均值：數值積分計算一個周期的平均值： 25,10 xyshier.m, shier1.mybxdxayrdydx)()(實例實例2弱肉強食弱肉強食 模型的解析解模型的解析解 ybxdyxayrx)()(ceyexayrbxd)(相軌線dyyayrdxxbxdc由初始條件確定由初始條件確定 相軌線是封閉曲線相軌線是

24、封閉曲線(c在一定范圍內在一定范圍內)求求x(t),y(t)一周期的平均值一周期的平均值:yx ,可以可以證明證明ybxdty)()(bdyytx/ )/()(TdttxTx0)(1bdx x(t), y(t)是周期函數是周期函數(周期記作周期記作T) )0(ln)(ln(1bdTbyTyT實例實例2弱肉強食弱肉強食 模型的解析解模型的解析解 ryax(t),y(t)一周期的平均值一周期的平均值:bdx 02. 0, 1 . 0, 5 . 0, 1badr10,25yxr 食餌增長率食餌增長率a 捕食者對食餌的捕獲能力捕食者對食餌的捕獲能力 d 捕食者死亡率捕食者死亡率b 食餌對捕食者的喂養(yǎng)能

25、力食餌對捕食者的喂養(yǎng)能力 結果解釋結果解釋ybxdyxayrx)()(與計算結果同與計算結果同yar ,xbd ,既相互制約既相互制約又相互依存又相互依存02040608010012005101520253000yx，00yx，00yx00yx0PT2T3T4T1P024681012020406080100120 x(t)y(t)T1 T2 T3 T4x(t) 的的“相位相位”領先領先 y(t)xayrtx)()(ybxdty)()()()(:1tytxT)()(:2tytxT)()(:3tytxT)()(:4tytxT進一步分析進一步分析)/,/(arbdP),(000yxP初值初值相軌線的

26、方向相軌線的方向一次大戰(zhàn)期間地中海漁業(yè)的捕撈量下降，一次大戰(zhàn)期間地中海漁業(yè)的捕撈量下降，但是其中但是其中鯊魚的比例卻在增加，為什么？鯊魚的比例卻在增加，為什么？rr- 1, dd+ 1捕撈捕撈戰(zhàn)時戰(zhàn)時捕撈捕撈rr- 2, dd+ 2 , 2 0微分方程穩(wěn)定微分方程穩(wěn)定 ),(1nnnnyxhfyy向前歐拉公式向前歐拉公式向后歐拉公式向后歐拉公式11nnnyyh y經典經典龍格龍格- -庫塔公式庫塔公式2.785/h算法穩(wěn)定算法穩(wěn)定)(,()(,(),(*yyyxfxxyxfyxfyyx0,yy(特征根特征根 - )nnh)1 (1nkn11h2/h111nnhh任意任意xycenyh )1 ( 剛性現象與剛性方剛性現象與剛性方程程bxaxrktfrxxkx) 0 (,) 0 (),0,()( 振動系統(tǒng)或電路系統(tǒng)的數學模型振動系統(tǒng)或電路系統(tǒng)的數學模型 現象現象 k=2000.5, r=1000, a=1, b=-1999.5, f(t)=1 瞬態(tài)解與穩(wěn)態(tài)解瞬態(tài)解與穩(wěn)態(tài)解 e-2000t 快瞬態(tài)解快瞬態(tài)解 e