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

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

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

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

5、x)y0P1P2P3歐 拉 方 法00( , ), ()yf x y y xy基本基本思路思路1,nnx x在小區(qū)間上1( , ),nnf x yx x中的x取內(nèi)的某一點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 x8P3歐 拉 方 法11()()( , ( ),nnnny xy xhf x

6、 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 9向前歐拉公式向前歐拉公式向后歐拉公式向后歐拉公式1(,)nnnnyyhf xy111(,)nnnnyyhfxy二者平均得到二者平均得到梯形公式梯形公式111(,)(,),0,1,2nnnnnnhyyfxyfxyn仍為隱式公

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

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

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

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

11、函數(shù)值ts = t0:k:tf輸出輸出 t0,tf 內(nèi)等分點處的函數(shù)值內(nèi)等分點處的函數(shù)值x0為函數(shù)初值為函數(shù)初值(n維維) 輸出輸出t=ts, x為相應(yīng)函數(shù)值為相應(yīng)函數(shù)值(n維維) optopt為選項，缺省時精度為：相對誤差為選項，缺省時精度為：相對誤差1010-3-3，絕對誤差絕對誤差1010-6-6, , 計算步長按精度要求自動調(diào)整計算步長按精度要求自動調(diào)整14例例 1 求 21 udtdu 的通解.微分方程的符號解微分方程的符號解 求微分方程（組）的解析解(符號解）命令:dsolve(方程方程1, 方程方程2,方程方程n, 初始條件初始條件, 自變量自變量)注意: 在表達(dá)微分方程時，用字

12、母D表示求微分，D2、D3等表示求高階微分.任何D后所跟的字母為因變量，自變量可以指定或由系統(tǒng)規(guī)則選定為確省.例如，微分方程 022dxyd應(yīng)表達(dá)為：D2y=0. 結(jié) 果：u = tan(t-c)t ,u21Dudsolve(的通解. 1dtdu求 1例2輸入命令解：u15例例 2 求微分方程的特解. 15)0( , 0)0(029422yyydxdydxyd 15)0( , 0)0(0294 .求微分方程的特解 2例22yyydxdydxyd 解解 輸入命令: y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x)結(jié) 果 為 : y =3e-2xsin（5x）

13、16zyxdtdzzyxdtdyzyxdtdx244354332 .求微分方程組 3例的解解解 輸入命令 ： x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z, t)； x=simple(x) % 將x化簡 y=simple(y) z=simple(z)結(jié) 果 為：x = (c1-c2+c3+c2e -3t-c3e-3t)e2t y = -c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z = (-c1e-4t+c2e-4t+c1-c2+c3)e2t 17. 12y) 10(1)0(yxexyy

14、xx已知精確解為的數(shù)值解。例：求微分方程解：用解：用MATLABMATLAB編程如下：編程如下：function dx=fun113(t,x)dx=t+x;ts=0:0.1:1;x0=1;t,x=ode45(fun113,ts,x0);y=2*exp(t)-t-1;t,x,y,x-yplot(t,x,r-,t,y,b-.),grid,gtext(x(t),gtext(y(t),title(微分方程數(shù)值解例微分方程數(shù)值解例);xlabel(time t);ylabel(solution x);legend(x,y);微分方程的數(shù)值解微分方程的數(shù)值解18微分方程的數(shù)值解微分方程的數(shù)值解2)2( ;

15、 2)2(022222yyynxdxdyxdxydx例：解解: 令 y1=x，y2=y11、建立m-文件fun114.m如下：function dx=fun114(t,x)dx=zeros(2,1);dx(1)=x(2);dx(2)=-x(2)/t+(1-0.25/t2)*x(1); 2、取初值，輸入命令： t,x=ode45(fun114,pi/2 2,pi/2 -2/pi plot(t,x(:,1),-)19實例實例1海上緝私海上緝私( (續(xù)續(xù)) ) 模型的數(shù)值解模型的數(shù)值解2222()()()()()()dxb cxdtcxatydyb atydtcxaty(0)0,(0)0 xy0yx

16、c(x(t),y(t)ab設(shè)設(shè): 船速船速a=20 (海里海里/小時小時) 艇速艇速b=40 (海里海里/小時小時) 距離距離c=15 (海里海里) 求求: 緝私艇的位置緝私艇的位置x(t),y(t) 緝私艇的航線緝私艇的航線y(x) jisi.m, seajisi.m20 實例實例1海上緝私海上緝私( (續(xù)續(xù)) ) 模型的數(shù)值解模型的數(shù)值解 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

17、.99840.06980.103.98540.29240.155.94450.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.021實例實例1 海上緝私海上緝私( (續(xù)續(xù)) ) 模型的數(shù)值解模型的數(shù)值解 設(shè)設(shè)b，c不變，不變，a變大變大為為30, 35, 接近接近40, 觀察解的變化觀察解的變化 :a=35, b=40, c=15t

18、=? 緝私艇追上走私船緝私艇追上走私船 t x(t) y(t) y1(t)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累積誤差較大累積誤差較大提高精度提高精度! ! 22實例實例1海上緝私海上緝私( (續(xù)續(xù)) ) 模型的數(shù)值解

19、模型的數(shù)值解 a=35, b=40, c=15opt=odeset(RelTol,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.07534

20、417.51.214.99961640.00000542.01.314.99996344.00000545.51.414.99999348.00000549.01.514.99999852.00000552.51.615.00002055.99993156.023實例實例1 海上緝私海上緝私( (續(xù)續(xù)) ) 模型的解析解模型的解析解 2222)()()(,)()()(yatxcyatbdtdyyatxcxcbdtdx)(xcyatdxdyatydxdyxc )(22()d ydtcxadxdxdsbdt22()()dsdxdy222()1 ()d yadycxdxbdxdxdyp 令bakpk

21、dxdpxc,1)(2(0)0pdydsdsdtdxdtdxdtdtdydxdy24實例實例1海上緝私海上緝私( (續(xù)續(xù)) ) 模型的解析解模型的解析解 21)(pkdxdpxc(0)0p21()kcxppc21()kcxppc1()() 2kkcxcxpcc/1ka b(0)0y11211()()2 111kkccxcxkcykckck緝私艇的航線緝私艇的航線y y( (x x) )的解析解的解析解x=c時時 緝私艇追上走私船的緝私艇追上走私船的y y坐標(biāo)坐標(biāo) 緝私艇追上走私船的時間緝私艇追上走私船的時間: 122bctbaa=20, b=40, c=15 t1=0.5 a=35, b=40

22、, c=15 t1=1.6 2221kcabcykba25微分方程數(shù)值解法的誤差分微分方程數(shù)值解法的誤差分析析數(shù)值解法數(shù)值解法: 計算微分方程精確解計算微分方程精確解y(xn)的近似值的近似值yn按照步長按照步長h一步步計算一步步計算, 每步都有誤差每步都有誤差;每一步的誤差會逐步積累每一步的誤差會逐步積累, , 稱稱累積誤差累積誤差. .討論計算一步出現(xiàn)的誤差討論計算一步出現(xiàn)的誤差假定假定yn= y(xn) , 估計估計yn+1的誤差的誤差: y(xn+1)- yn+1局部截斷誤差局部截斷誤差 26誤差分析誤差分析估計歐拉公式的局部截斷誤差估計歐拉公式的局部截斷誤差 y(xn+1)在在xn處

23、作處作Taylor展開展開: )()(2)()()(321hOxyhxyhxyxynnnn 向前歐拉公式向前歐拉公式1(,)nnnnyyhf xyyn= y(xn)()()(,()(1nnnnnnxyhxyxyxhfxyy232111()()()()2nnnnhTy xyyxO hO h局部截斷誤差主項為局部截斷誤差主項為27誤差分析誤差分析估計歐拉公式的局部截斷誤差估計歐拉公式的局部截斷誤差 )()(2)()()(321hOxyhxyhxyxynnnn 向后歐拉公式向后歐拉公式),(111nnnnyxhfyy232111()()()()2nnnnhTy xyyxO hO h 局部截斷誤差主項

24、為局部截斷誤差主項為)(22nxyh 梯形梯形公式公式向前、向后歐拉公式的平均向前、向后歐拉公式的平均3431()()()12nnhTyxO hO h xn xn+1 x y Pn A Q B28誤誤差差分分析析算法算法精度的階精度的階的定義的定義 一個算法的局部截斷誤差為一個算法的局部截斷誤差為O(hp+1) 該算法具有該算法具有p階精度階精度 局部截斷誤差局部截斷誤差精度精度向前歐拉公式向前歐拉公式O(h2)1階階向后歐拉公式向后歐拉公式O(h2)1階階梯形公式梯形公式O(h3)2階階改進歐拉公式改進歐拉公式O(h3)2階階經(jīng)典龍格經(jīng)典龍格- -庫塔公式庫塔公式O(h5)4階階29實實 例

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

26、卻增加，為什么？鯊魚的比例卻增加，為什么？30實實 例例 2弱弱 肉肉 強強 食食 模型模型食餌食餌(甲甲) 的密度的密度x(t), 捕食者捕食者(乙乙)的密度的密度y(t)0,/rrxx 甲獨立生存的增長率甲獨立生存的增長率r r0,/aayrxx 乙使甲的增長率減小乙使甲的增長率減小, ,減小量與減小量與 y y 成正比成正比0,/ddyy 乙獨立生存的死亡率乙獨立生存的死亡率d0),(/bbxdyy 甲使乙的死亡率減小甲使乙的死亡率減小, ,減小量與減小量與 x成正比成正比00()()(0),(0)xray xrxaxyydbx ydybxyxxyyVolterra模型模型 x(t),

27、y(t)無解析解無解析解31實例實例2弱肉強食弱肉強食 模型的數(shù)值解模型的數(shù)值解 00)0(,)0(yyxxbxydyyaxyrxx001,0.50.1,0.0225,2rdabxy051015020406080100 x(t)y(t)020406080100051015202530 xy猜測猜測 x(t), y(t)是周期函數(shù)是周期函數(shù); y(x)是封閉曲線是封閉曲線 數(shù)值積分計算一個周期的平均值：數(shù)值積分計算一個周期的平均值： 25,10 xyshier.m, shier1.m32ybxdxayrdydx)()(實例實例2弱肉強食弱肉強食 模型的解析解模型的解析解 ybxdyxayrx)(

28、)(ceyexayrbxd)(相軌線dyyayrdxxbxdc由初始條件確定由初始條件確定 相軌線是封閉曲線相軌線是封閉曲線(c在一定范圍內(nèi)在一定范圍內(nèi))求求x(t),y(t)一周期的平均值一周期的平均值:yx ,可以可以證明證明ybxdty)()(bdyytx/ )/()(TdttxTx0)(1bdx x(t), y(t)是周期函數(shù)是周期函數(shù)(周期記作周期記作T) )0(ln)(ln(1bdTbyTyT33實例實例2弱肉強食弱肉強食 模型的解析解模型的解析解 ryax(t),y(t)一周期的平均值一周期的平均值:bdx 02. 0, 1 . 0, 5 . 0, 1badr10,25yxr 食

29、餌增長率食餌增長率a 捕食者對食餌的捕獲能力捕食者對食餌的捕獲能力 d 捕食者死亡率捕食者死亡率b 食餌對捕食者的喂養(yǎng)能力食餌對捕食者的喂養(yǎng)能力 結(jié)果解釋結(jié)果解釋ybxdyxayrx)()(與計算結(jié)果同與計算結(jié)果同yar ,xbd ,既相互制約既相互制約又相互依存又相互依存3402040608010012005101520253000yx，00yx，00yx00yx0PT2T3T4T1P024681012020406080100120 x(t)y(t)T1 T2 T3 T4x(t) 的的“相位相位”領(lǐng)先領(lǐng)先 y(t)xayrtx)()(ybxdty)()()()(:1tytxT)()(:2ty

30、txT)()(:3tytxT)()(:4tytxT進一步分析進一步分析)/,/(arbdP),(000yxP初值初值相軌線的方向相軌線的方向35一次大戰(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經(jīng)典經(jīng)典龍格龍格- -庫塔公式庫塔公式2.785/h算法穩(wěn)定算法穩(wěn)定)(,()(,(),(*yyyxfxxyxfyxfyyx0,yy(特征根特征根 - )nnh)1 (1nkn11h2/h111nnhh任意任意xycenyh )1 (40 剛性現(xiàn)象與剛性方剛性現(xiàn)象與剛性方程程bxaxrktfrxxkx)