《數(shù)值與算法》第八講-常微分方程初值問(wèn)題

1、數(shù)數(shù) 值值 分分 析析 (8)Numerical AnalysisWenjian Yu2第八章第八章 常微分方程初值問(wèn)題常微分方程初值問(wèn)題n(主要是前主要是前3節(jié)節(jié))常微分方程初值問(wèn)題常微分方程初值問(wèn)題Wenjian Yu3常微分方程基本概念常微分方程基本概念 常微分方程常微分方程nWenjian Yu4常微分方程常微分方程nWenjian Yu5初值問(wèn)題初值問(wèn)題:常微分方程常微分方程 例子例子1n含電容元件的電路問(wèn)題含電容元件的電路問(wèn)題n通過(guò)通過(guò)節(jié)點(diǎn)分析法節(jié)點(diǎn)分析法得到微分方程組得到微分方程組Wenjian Yu6問(wèn)題的解反映了電容充問(wèn)題的解反映了電容充/放電過(guò)程放電過(guò)程C1R2R1R3C2

2、1Ri1Ci2Ri2Ci3Ri1v2v電流電流/ /電壓關(guān)系電壓關(guān)系ddCCviCt1RRivR節(jié)點(diǎn)電流方程節(jié)點(diǎn)電流方程1120CRRiii2320CRRiii112211222232111d0d0111ddvRRRCvtCvvRRRt 1(0)1v2(0)0vt常微分方程常微分方程 例子例子2n雙聯(lián)擺的運(yùn)動(dòng)雙聯(lián)擺的運(yùn)動(dòng)n兩個(gè)擺錘兩個(gè)擺錘(重物重物), 剛性桿的重量可忽略剛性桿的重量可忽略n不考慮摩擦力不考慮摩擦力,運(yùn)動(dòng)不會(huì)停止運(yùn)動(dòng)不會(huì)停止, 且在初始角且在初始角度較大時(shí)擺錘的軌跡呈現(xiàn)混沌現(xiàn)象度較大時(shí)擺錘的軌跡呈現(xiàn)混沌現(xiàn)象Wenjian Yu7求解微分方程初值問(wèn)題求解微分方程初值問(wèn)題, 得到得

3、到擺錘的運(yùn)動(dòng)規(guī)律擺錘的運(yùn)動(dòng)規(guī)律 Matlab演示演示swinger常微分方程常微分方程nWenjian Yu8線性齊次常系數(shù)微分方程線性齊次常系數(shù)微分方程實(shí)際的問(wèn)題基本上都是穩(wěn)定的實(shí)際的問(wèn)題基本上都是穩(wěn)定的!(由于歷史原因由于歷史原因)常微分方程常微分方程nWenjian Yu9局部穩(wěn)定局部穩(wěn)定簡(jiǎn)單方法與有關(guān)概念簡(jiǎn)單方法與有關(guān)概念Wenjian Yu10簡(jiǎn)單的初值問(wèn)題數(shù)值解法簡(jiǎn)單的初值問(wèn)題數(shù)值解法 初值問(wèn)題的數(shù)值解法初值問(wèn)題的數(shù)值解法nWenjian Yu11否則為否則為多步法多步法否則為否則為隱格式方法隱格式方法歐拉法歐拉法nWenjian Yu12“左矩形左矩形”求積公式求積公式h=0.1

4、h=0.050.1 1.000000 1.004837 0.05 1.000000 0.3 1.0350920.2 1.010000 1.018731 0.1 1.002500 0.35 1.0483370.3 1.029000 1.040818 0.15 1.007375 0.4 1.0634200.4 1.056100 1.070320 0.2 1.014506 0.45 1.0802490.5 1.090490 1.106531 0.25 1.023781 0.5 1.098737步長(zhǎng)步長(zhǎng)h=0.1, 和和0.05步長(zhǎng)步長(zhǎng)小的更準(zhǔn)小的更準(zhǔn)數(shù)值解法的穩(wěn)定性數(shù)值解法的穩(wěn)定性nWenjian

5、Yu13-1 0歐拉法解模型問(wèn)題歐拉法解模型問(wèn)題的的穩(wěn)定區(qū)域穩(wěn)定區(qū)域數(shù)值解法的穩(wěn)定性數(shù)值解法的穩(wěn)定性nWenjian Yu14-1 0歐拉法穩(wěn)定歐拉法穩(wěn)定數(shù)值解法的穩(wěn)定性數(shù)值解法的穩(wěn)定性nWenjian Yu1500.0250.050.0750.10.1250.151-1.52.25-3.3755.0625-7.59375 11.39061 0.082085 0.006738 0.000553 4.5410-53.7310-63.0610-7這里設(shè)的這里設(shè)的h太大太大!計(jì)算結(jié)果如下表計(jì)算結(jié)果如下表:數(shù)值解法的局部截?cái)嗾`差數(shù)值解法的局部截?cái)嗾`差nWenjian Yu16整體誤差整體誤差穩(wěn)定的問(wèn)題

6、，整體誤差小穩(wěn)定的問(wèn)題，整體誤差小于局部誤差之和于局部誤差之和不穩(wěn)定的問(wèn)題呢？不穩(wěn)定的問(wèn)題呢？一般僅能控制局部誤差一般僅能控制局部誤差整體誤差？整體誤差？ 局部誤差局部誤差數(shù)值解法的局部截?cái)嗾`差數(shù)值解法的局部截?cái)嗾`差nWenjian Yu17歐拉法是歐拉法是一階方法一階方法我們討論的所有方法都至少有我們討論的所有方法都至少有1階準(zhǔn)確度階準(zhǔn)確度數(shù)值解法的數(shù)值解法的收斂性收斂性：隨著：隨著h0，誤差，誤差0向后歐拉法與梯形法向后歐拉法與梯形法n從數(shù)值積分的角度推導(dǎo)從數(shù)值積分的角度推導(dǎo)n向后歐拉法向后歐拉法:n梯形法梯形法:n兩者均為單步、兩者均為單步、隱格式隱格式方法方法, 每每步步計(jì)算要求解計(jì)算

7、要求解(非線性非線性)方程方程n例例8.6:用向后歐拉法求解用向后歐拉法求解Wenjian Yu18右矩形右矩形梯形梯形00.0250.050.0750.10.1250.1510.006663 0.001904 0.0005441 0.082085 0.006738 0.000553 4.5410-53.7310-63.0610-7向后歐拉法向后歐拉法nWenjian Yu19準(zhǔn)確解準(zhǔn)確解0 1穩(wěn)定區(qū)域穩(wěn)定區(qū)域無(wú)條件穩(wěn)定無(wú)條件穩(wěn)定(unconditionally stable)!向后歐拉法向后歐拉法nWenjian Yu20具有具有1階準(zhǔn)確度階準(zhǔn)確度!向后歐拉法與梯形法向后歐拉法與梯形法nWe

8、njian Yu21穩(wěn)定的條件是穩(wěn)定的條件是: 思考思考無(wú)條件穩(wěn)定無(wú)條件穩(wěn)定!具有具有2階準(zhǔn)確度階準(zhǔn)確度簡(jiǎn)單方法與有關(guān)概念簡(jiǎn)單方法與有關(guān)概念Wenjian Yu22Runge-Kutta方法方法 n在歐拉法基礎(chǔ)上改進(jìn)在歐拉法基礎(chǔ)上改進(jìn)再增加一次函數(shù)求值再增加一次函數(shù)求值: 數(shù)值積分的數(shù)值積分的中矩形中矩形或或梯形公式梯形公式利用歐拉法算半個(gè)步長(zhǎng)的結(jié)果利用歐拉法算半個(gè)步長(zhǎng)的結(jié)果, 估算中點(diǎn)處被積函數(shù)值估算中點(diǎn)處被積函數(shù)值先用歐拉法估計(jì)區(qū)間終點(diǎn)處斜率先用歐拉法估計(jì)區(qū)間終點(diǎn)處斜率, 再用它與再用它與起始點(diǎn)斜率的平均值算一整步起始點(diǎn)斜率的平均值算一整步Runge-Kutta方法方法( , ( )n n

9、n nththnnnnt tyyf s y s dsyyf s y s ds1 1中中矩形矩形梯梯形公式形公式(中中點(diǎn)公式點(diǎn)公式)(Heun方法方法/改進(jìn)的歐拉法改進(jìn)的歐拉法)都比歐拉法準(zhǔn)確都比歐拉法準(zhǔn)確均為均為2級(jí)級(jí)R-K公式公式Wenjian Yu23nRunge-Kutta方法方法( , ( )n nn nththnnnnt tyyf s y s dsyyf s y s ds1 1(, ()r rnninininnininii iyyhc f th y thyyhc f th y th1 11 1()niniy thy th只能估算只能估算() nnnny tyy ty(,) nnnnk

10、f t ykf t y1 1, 令令, , , r rkkkkkk2323r rnni inni ii iyyhc kyyhc k1 11 1Wenjian Yu24nRunge-Kutta方法方法r rnni inni ii iyyhc kyyhc k1 11 1()niniy thy th其他其他, 用所有前面點(diǎn)的信息用所有前面點(diǎn)的信息Wenjian Yu25積分節(jié)點(diǎn)積分節(jié)點(diǎn)n幾種顯式幾種顯式R-K公式公式參數(shù)的值不按具體數(shù)值參數(shù)的值不按具體數(shù)值積分公式設(shè)置積分公式設(shè)置, 而根據(jù)而根據(jù)準(zhǔn)確度階數(shù)準(zhǔn)確度階數(shù)要求設(shè)置要求設(shè)置2級(jí)公式：改進(jìn)歐拉法、級(jí)公式：改進(jìn)歐拉法、中點(diǎn)公式中點(diǎn)公式經(jīng)典經(jīng)典4級(jí)

11、、級(jí)、4階階Runge-Kutta法法 (1905)Runge-Kutta方法方法Wenjian Yu26nRunge-Kutta方法方法此時(shí)此時(shí)局部截局部截?cái)嗾`差斷誤差只要只要對(duì)非模型問(wèn)題也有相同結(jié)論！對(duì)非模型問(wèn)題也有相同結(jié)論！Wenjian Yu27如如n例例8.7: 用用2階改進(jìn)歐拉、階改進(jìn)歐拉、3階階Ralston、4階經(jīng)典階經(jīng)典R-K解問(wèn)題解問(wèn)題, h=0.1, 算到算到y(tǒng)(2)精確解為精確解為Runge-Kutta方法方法r4對(duì)應(yīng)的對(duì)應(yīng)的r級(jí)級(jí)R-K公式公式有有r階準(zhǔn)確度階準(zhǔn)確度高于高于4階的階的公式很少單公式很少單獨(dú)使用獨(dú)使用Wenjian Yu28r4對(duì)應(yīng)的對(duì)應(yīng)的r級(jí)級(jí)R-K公

12、式公式達(dá)不到達(dá)不到r階階準(zhǔn)確度準(zhǔn)確度二階二階Heun 三階三階Ralston 四階四階Runge-Kutta準(zhǔn)確值準(zhǔn)確值10.40.40.40.41.10.4756410.4746260.47463830.47463821.20.5834080.5813640.58138680.58138671.30.7281350.7250340.72506630.72506621.40.9153290.9111370.91117730.91117711.51.1511101.1457851.14583361.14583331.61.4421691.4356641.43572031.43572001.71.

13、7957381.7880041.78806741.78806711.82.2195782.2105612.21063152.21063111.92.7219612.7116062.71168362.71168322.03.3116653.2999163.30000043.3000000nRunge-Kutta方法方法(此時(shí)此時(shí) )一般顯式單步法一般顯式單步法恰好與歐拉法一樣恰好與歐拉法一樣顯格式顯格式, 都都不是不是無(wú)條件無(wú)條件穩(wěn)定的穩(wěn)定的局部局部截?cái)嘟財(cái)嗾`差誤差判斷單步法判斷單步法收斂性收斂性的簡(jiǎn)便方法的簡(jiǎn)便方法Wenjian Yu29簡(jiǎn)單方法與有關(guān)概念簡(jiǎn)單方法與有關(guān)概念Wenjian Yu

14、30多步法多步法 n多步法多步法Wenjian Yu31(線性線性m步法步法)固定步長(zhǎng)固定步長(zhǎng)hnTaylor展開(kāi)法求展開(kāi)法求線性線性m步法步法的系數(shù)的系數(shù) 這些系數(shù)應(yīng)該等于這些系數(shù)應(yīng)該等于0例例8.8: 求兩步法公式求兩步法公式 中參數(shù)值中參數(shù)值多步法多步法Wenjian Yu32滿足它們才可能滿足它們才可能收斂收斂 (相容性相容性), 有有二階準(zhǔn)確度二階準(zhǔn)確度n多步法多步法Wenjian Yu33同同例例8.8的結(jié)果的結(jié)果!多步法公式中包含多步法公式中包含p個(gè)待定參數(shù)，至少可達(dá)到個(gè)待定參數(shù)，至少可達(dá)到p-1階準(zhǔn)確度階準(zhǔn)確度n多步法多步法Wenjian Yu34Vandermonde陣陣T,

15、 非奇異非奇異插值節(jié)點(diǎn)插值節(jié)點(diǎn)函數(shù)值函數(shù)值n常用的多步法公式常用的多步法公式Adams公式公式的推導(dǎo)：的推導(dǎo)：用插值多項(xiàng)式近似被積函數(shù)用插值多項(xiàng)式近似被積函數(shù)例例8.10: 推導(dǎo)推導(dǎo)m=4對(duì)應(yīng)的顯式對(duì)應(yīng)的顯式Adams公式公式多步法多步法Wenjian Yu35可證明滿足最高準(zhǔn)確度階數(shù)可證明滿足最高準(zhǔn)確度階數(shù)類(lèi)似地算其類(lèi)似地算其他系數(shù)，得他系數(shù)，得顯式四階顯式四階Adams-Bashforth公式(單項(xiàng)式函數(shù)代入法單項(xiàng)式函數(shù)代入法)nAdams公式公式幾種顯式公式幾種顯式公式幾種隱式公式幾種隱式公式多步法多步法Wenjian Yu36階數(shù)階數(shù)穩(wěn)定閾值穩(wěn)定閾值 誤差常數(shù)誤差常數(shù)11 -21/22

16、3/2-1/2 -15/12323/12 -16/12 5/12 -6/113/8455/24 -59/24 37/24 -9/24-3/10251/720歐拉法歐拉法階數(shù)階數(shù)穩(wěn)定閾值穩(wěn)定閾值 誤差常數(shù)誤差常數(shù)11 -1/221/21/2 -1/1235/128/12-1/12 -6-1/2449/2419/24 -5/241/24-3-19/720向后歐拉法向后歐拉法梯形法梯形法并非無(wú)條件并非無(wú)條件穩(wěn)定穩(wěn)定!n多步法多步法Wenjian Yu37用用Matlab求解初值問(wèn)題求解初值問(wèn)題Wenjian Yu38用用Matlab解解ODE-IVP Matlab中的中的ODE-IVP求解器求解器n

17、Wenjian Yu39T, Y, TE, YE, IE = solver(odefun, tspan, y0, options) 求解單個(gè)求解單個(gè)ODEn火焰燃燒問(wèn)題火焰燃燒問(wèn)題當(dāng)當(dāng)點(diǎn)燃一根火柴點(diǎn)燃一根火柴時(shí)時(shí), 火焰火焰迅速增大直到一個(gè)迅速增大直到一個(gè)臨界體積臨界體積, 然然后后維持這一體積維持這一體積不變不變, 此時(shí)此時(shí)火焰內(nèi)部燃燒耗費(fèi)的氧氣和火焰內(nèi)部燃燒耗費(fèi)的氧氣和其表面現(xiàn)存的氧氣達(dá)到了一種其表面現(xiàn)存的氧氣達(dá)到了一種平衡平衡.火焰火焰(近似為球近似為球)半徑半徑y(tǒng)滿足滿足ODE f= (t,y) y2-y3; ode23(f, 0, 2.0e4, 1e-4)Wenjian Yu4023

18、-(0)yyyy 設(shè)初始半徑設(shè)初始半徑=0.000100.511.52x 10400.20.40.60.811.21.4例例8.15嘗試用嘗試用ode23s求解它求解它00.10.20.30.40.50.60.70.80.9100.511.522.533.544.5求解求解ODE方程組方程組nWenjian Yu41例例8.1612122 , (0)0,(0)16yyyyyt function ydot = myode2(t, y);ydot= y(2); 6*t; %列向量列向量雙聯(lián)擺問(wèn)題的求解雙聯(lián)擺問(wèn)題的求解Wenjian Yu42隱格式非線性常微分方程初值隱格式非線性常微分方程初值問(wèn)題問(wèn)題, 初值為初值為, , 0, 0Todeset可以設(shè)置質(zhì)量矩陣可以設(shè)置質(zhì)量矩陣質(zhì)量矩陣質(zhì)量矩陣swinger_solve.m很大時(shí)混沌現(xiàn)象很大時(shí)混沌現(xiàn)象: