工程數(shù)理建模第九章

1、第九章第九章 常微分方程數(shù)值解法常微分方程數(shù)值解法理學(xué)院理學(xué)院 2013年年4月月工程數(shù)學(xué)模型及數(shù)值方法工程數(shù)學(xué)模型及數(shù)值方法優(yōu)質(zhì)課程優(yōu)質(zhì)課程 許多氣井都不同程度地含有液體許多氣井都不同程度地含有液體.對(duì)于存在對(duì)于存在底水或邊底水或邊水水的氣藏，在開(kāi)采過(guò)程中液氣比將逐漸增高，會(huì)明顯地的氣藏，在開(kāi)采過(guò)程中液氣比將逐漸增高，會(huì)明顯地影響氣井的產(chǎn)能，甚至將氣井淹死影響氣井的產(chǎn)能，甚至將氣井淹死.因此，正確地預(yù)測(cè)氣因此，正確地預(yù)測(cè)氣井在井在較高含液程度下的舉升能力較高含液程度下的舉升能力，對(duì)于氣井動(dòng)態(tài)分析和，對(duì)于氣井動(dòng)態(tài)分析和排水采氣（如氣舉）設(shè)計(jì)具有重要的實(shí)際意義排水采氣（如氣舉）設(shè)計(jì)具有重要的實(shí)際

2、意義. 盡管流體力學(xué)的基本方程也適用于油氣水多相流動(dòng)，盡管流體力學(xué)的基本方程也適用于油氣水多相流動(dòng)，不過(guò)在解決采油或采氣工程技術(shù)問(wèn)題時(shí)，一般不過(guò)在解決采油或采氣工程技術(shù)問(wèn)題時(shí)，一般把油水兩把油水兩種液體視為液相種液體視為液相，著重考慮，著重考慮氣液兩相間的作用氣液兩相間的作用.描述兩相描述兩相管流的數(shù)學(xué)模型比單相管流管流的數(shù)學(xué)模型比單相管流復(fù)雜復(fù)雜得多得多.引言兩相管流數(shù)值模型兩相管流數(shù)值模型 由于流體的非均質(zhì)性，在氣液兩相管流中，氣液各相由于流體的非均質(zhì)性，在氣液兩相管流中，氣液各相的分布狀況可能是多種多樣的，存在著各種不同的流動(dòng)的分布狀況可能是多種多樣的，存在著各種不同的流動(dòng)形態(tài)，而氣液界

3、面又很復(fù)雜和多變形態(tài)，而氣液界面又很復(fù)雜和多變.因此，尋求實(shí)用的、因此，尋求實(shí)用的、嚴(yán)格的數(shù)學(xué)解是很困難的嚴(yán)格的數(shù)學(xué)解是很困難的.對(duì)于采氣工程中的氣液兩相管對(duì)于采氣工程中的氣液兩相管流，其流，其核心問(wèn)題是探討沿程壓力損失及影響因素核心問(wèn)題是探討沿程壓力損失及影響因素. 60-70年代，一般的處理方法是從物理概念和基本方程出發(fā)，年代，一般的處理方法是從物理概念和基本方程出發(fā)，采用采用實(shí)驗(yàn)和因次分析法實(shí)驗(yàn)和因次分析法得到描述某一特定兩相管流過(guò)程得到描述某一特定兩相管流過(guò)程的一些無(wú)因次參數(shù)，然后根據(jù)實(shí)驗(yàn)數(shù)據(jù)得出的一些無(wú)因次參數(shù)，然后根據(jù)實(shí)驗(yàn)數(shù)據(jù)得出經(jīng)驗(yàn)關(guān)系式經(jīng)驗(yàn)關(guān)系式.引言兩相管流數(shù)值模型兩相管流數(shù)

4、值模型 MukherjeeMukherjee和和BrillBrill（19851985）在前人研究工作的基礎(chǔ)上，）在前人研究工作的基礎(chǔ)上，改進(jìn)實(shí)驗(yàn)條件，提出了更為實(shí)用的傾斜管（包括水平管改進(jìn)實(shí)驗(yàn)條件，提出了更為實(shí)用的傾斜管（包括水平管）兩相流的流型判別準(zhǔn)則和應(yīng)用方便的持液率及摩阻系）兩相流的流型判別準(zhǔn)則和應(yīng)用方便的持液率及摩阻系數(shù)經(jīng)驗(yàn)公式數(shù)經(jīng)驗(yàn)公式.MB.MB模型的壓降梯度方程為模型的壓降梯度方程為引言221sin/mmmmmmsggfvDdpdzv vp 式中，式中， 為油管內(nèi)徑為油管內(nèi)徑.對(duì)于油套環(huán)空流動(dòng)，對(duì)于油套環(huán)空流動(dòng)， 為水力當(dāng)量為水力當(dāng)量直徑（套管內(nèi)徑和油管外徑之差）直徑（套管內(nèi)徑

5、和油管外徑之差）. 為兩相摩阻系數(shù)，為兩相摩阻系數(shù)， 為氣液混合物平均密度，為氣液混合物平均密度， 為持液率為持液率. 氣相表觀速度氣相表觀速度.DDmf1 ()mllglHH lHsgv兩相管流數(shù)值模型兩相管流數(shù)值模型引言數(shù)學(xué)上的一般表達(dá)式為數(shù)學(xué)上的一般表達(dá)式為若已知起始點(diǎn)若已知起始點(diǎn) （井口或井底）處的流壓（井口或井底）處的流壓 ，聯(lián)合，聯(lián)合上述方程，就構(gòu)成了一個(gè)常微分方程的初值問(wèn)題上述方程，就構(gòu)成了一個(gè)常微分方程的初值問(wèn)題. .其一般其一般形式為形式為0z0p00( ,)()dpF z pdzp zp 00( , )()dyf x ydxy xy 兩相管流數(shù)值模型兩相管流數(shù)值模型Nume

6、rical Method for Ordinary Differential Equations常用的一些常用的一些解析解法解析解法：常數(shù)常數(shù)變易變易法、法、LapalaceLapalace變換等變換等分離變量分離變量法、變量法、變量代換代換、一階一階常微分方程初值問(wèn)題：常微分方程初值問(wèn)題：00( , );()dyf x y axbdxy xy ( ) 常微分方程數(shù)值解法 在科學(xué)研究及工程技術(shù)領(lǐng)域中，常常會(huì)遇到大量的如在科學(xué)研究及工程技術(shù)領(lǐng)域中，常常會(huì)遇到大量的如上述的常微分方程的求解問(wèn)題上述的常微分方程的求解問(wèn)題.除了一些簡(jiǎn)單的方程外，除了一些簡(jiǎn)單的方程外，要用傳統(tǒng)的數(shù)學(xué)分析方法找出復(fù)雜的變

7、系數(shù)或非線性問(wèn)要用傳統(tǒng)的數(shù)學(xué)分析方法找出復(fù)雜的變系數(shù)或非線性問(wèn)題的解析表達(dá)式是困難的，有時(shí)甚至是不可能的題的解析表達(dá)式是困難的，有時(shí)甚至是不可能的.同時(shí)許同時(shí)許多實(shí)際問(wèn)題也只需要獲得解在若干個(gè)點(diǎn)上的近似值即可多實(shí)際問(wèn)題也只需要獲得解在若干個(gè)點(diǎn)上的近似值即可.因此，研究和掌握常常微分方程數(shù)值解法，即求出解在因此，研究和掌握常常微分方程數(shù)值解法，即求出解在一系列離散點(diǎn)上的解的近似值的方法，是很有必要的一系列離散點(diǎn)上的解的近似值的方法，是很有必要的.本章主要介紹本章主要介紹一階一階常微分方程初值問(wèn)題的常微分方程初值問(wèn)題的數(shù)值解法數(shù)值解法。常微分方程數(shù)值解法對(duì)于初值問(wèn)題對(duì)于初值問(wèn)題 ，如果，如果 在下

8、列區(qū)域內(nèi)連續(xù)：在下列區(qū)域內(nèi)連續(xù)：( , )f x y( ) （解的（解的存在唯一存在唯一性）性）;|Gaxby 且關(guān)于且關(guān)于 滿足滿足LipschitzLipschitz條件，即存在常數(shù)條件，即存在常數(shù) ，使，使y0L 1212|( ,)( ,)|;,f x yf x yL yyx yG則初值問(wèn)題則初值問(wèn)題 存在唯一解，且解是存在唯一解，且解是連續(xù)可微連續(xù)可微的。的。( ) 所謂所謂數(shù)值解數(shù)值解是指：在解的是指：在解的存在區(qū)間存在區(qū)間上取一系列點(diǎn)上取一系列點(diǎn)012.nxxxx逐個(gè)求出逐個(gè)求出 的近似值的近似值1 2 3(, , ,.)iy i ()iy x0;ixxih 等距等距節(jié)點(diǎn)節(jié)點(diǎn)：:h

9、步長(zhǎng)步長(zhǎng)常微分方程數(shù)值解法定理定理1定義定義1初值問(wèn)題初值問(wèn)題 的解析解及其數(shù)值解的的解析解及其數(shù)值解的幾何幾何意義：意義：( ) oxy初值問(wèn)題初值問(wèn)題 的解表示過(guò)點(diǎn)的解表示過(guò)點(diǎn) 的一條的一條曲線曲線( ) 00(,)xynx (,)nnxy ),(00yx( )yy x 0 x),(00yx 2x),(22yx 1x),(11yx 初值問(wèn)題初值問(wèn)題 的數(shù)值解表示一組的數(shù)值解表示一組離散點(diǎn)列離散點(diǎn)列( ) (,)iixy可用可用擬合擬合方法求該組數(shù)據(jù)方法求該組數(shù)據(jù) 的的近似曲線近似曲線(,)iixy積分積分曲線曲線常微分方程數(shù)值解法EulerEuler方法的方法的導(dǎo)出導(dǎo)出212()()()(

10、)!nnnnh yy xy xhy x 將將 在點(diǎn)在點(diǎn) 處進(jìn)行處進(jìn)行Taylor展開(kāi)展開(kāi)1()ny x nx略去略去 項(xiàng)：項(xiàng)：2h然后用然后用 代替代替 ，即得，即得ny()ny x1()()(, ()nnnny xy xhf xy x 10 1 21(,), , ,nnnnyyhf xynN 稱上述公式為稱上述公式為向前向前Euler 公式。公式。Euler法2112()()()()!nnnnh yy xy xhy x 若將若將 在點(diǎn)在點(diǎn) 處進(jìn)行處進(jìn)行Taylor展開(kāi)展開(kāi)()ny x1nx 略去略去 項(xiàng)：項(xiàng)：2h然后用然后用 代替代替 ，即得，即得ny()ny x111()()(, ()nn

11、nny xy xhf xy x 1110 1 21, ,(,),nnnnyyhfynxN稱上述公式為稱上述公式為向后向后Euler 公式。公式。向后向后Euler 公式為公式為隱式隱式格式，需要利用格式，需要利用迭代法迭代法求解求解Euler法解向前向前EulerEuler公式：公式：1(,)nnnnyyhf xy （取步長(zhǎng)為取步長(zhǎng)為 ）101( )dyxydxy 例例1 分別利用向前和向后分別利用向前和向后EulerEuler方法求解初值問(wèn)題的方法求解初值問(wèn)題的數(shù)值解數(shù)值解0 1 .h 10 10 90 1.nnnyxy 向后向后EulerEuler公式：公式：111(,)nnnnyyhf

12、xy1110 10 11 1( . ).nnnyxyEuler法Euler法解向前向前EulerEuler公式：公式：1(,)nnnnyyhf xy （取步長(zhǎng)為取步長(zhǎng)為 ）201( )dyxydxyy 例例2 分別利用向前和向后分別利用向前和向后EulerEuler方法求解初值問(wèn)題的方法求解初值問(wèn)題的數(shù)值解數(shù)值解0 1 .h 12()nnnnnxyyh yy 向后向后EulerEuler公式：公式：111(,)nnnnyyhf xy11112()nnnnnxyyh yy Euler法取步長(zhǎng)取步長(zhǎng) ，這時(shí)，這時(shí) ，計(jì)算到，計(jì)算到 ，誤差，誤差 采用歐拉法計(jì)算采用歐拉法計(jì)算 0 1 .h 00 1

13、 .nxxnhn 10n ()nnny xy 0.11.10001.09540.00460.61.50901.48320.02580.21.19181.18320.00860.71.58031.54920.03110.31.27741.26490.01250.81.64981.61250.03730.41.35821.34160.01660.91.71781.67330.04450.51.43511.141420.02091.01.78481.73210.0527nxny()ny x|n nxny()ny x|n Euler法隱式格式通常采用迭代法計(jì)算，迭代過(guò)程的實(shí)質(zhì)就是將隱式格式通常采用迭代

14、法計(jì)算，迭代過(guò)程的實(shí)質(zhì)就是將隱式格式逐步顯式化。隱式格式逐步顯式化。兩種格式計(jì)算差異兩種格式計(jì)算差異： EulerEuler法和后退法和后退EulerEuler法有本質(zhì)的區(qū)別。前者是關(guān)于法有本質(zhì)的區(qū)別。前者是關(guān)于 的一個(gè)直接計(jì)算公式，這種格式稱為的一個(gè)直接計(jì)算公式，這種格式稱為顯式顯式的；后者一般的；后者一般是關(guān)于是關(guān)于 的一個(gè)非線性方程，這種格式稱為的一個(gè)非線性方程，這種格式稱為隱式隱式的。的。1ny 1ny Euler法稱為稱為梯形法梯形法。易見(jiàn)梯形法也是隱式法，需要采用迭代計(jì)。易見(jiàn)梯形法也是隱式法，需要采用迭代計(jì)算。算。梯形法將將EulerEuler法和后退法和后退EulerEuler法

15、加權(quán)平均得到：法加權(quán)平均得到：111(,)(,)2nnnnnnhyyf xyf xy 采用顯式的采用顯式的EulerEuler法提供迭代初值，梯形法的迭代公式為法提供迭代初值，梯形法的迭代公式為01111120 1 2( )()( )(,)(,)(,)(, , ,)nnnnkknnnnnnyyhf xyhyyf xyf xyk 預(yù)測(cè)：預(yù)測(cè)：改進(jìn)Euler法將利用梯形法迭代一次的格式稱為將利用梯形法迭代一次的格式稱為改進(jìn)改進(jìn)EulerEuler法，也稱為預(yù)測(cè)法，也稱為預(yù)測(cè)校正法校正法：校正：校正：1(,)nnnnyyhf xy 1112(,)(,)nnnnnnhyyf xyf xy 或者表示成：

16、或者表示成：11(,)(,)1()2pnnncnnpnpcyyhf xyyyhf xyyyy 解向前向前EulerEuler公式：公式： （取步長(zhǎng)為取步長(zhǎng)為 ）101( )dyxydxy 例例3：分別利用分別利用四種方法四種方法求解初值問(wèn)題的求解初值問(wèn)題的數(shù)值解數(shù)值解0 1 .h 10 10 90 1.nnnyxy 向后向后EulerEuler公式：公式：1110 10 11 1( . ).nnnyxy梯形法：梯形法：10 0950 9050 1.nnnyxy 改進(jìn)改進(jìn)EulerEuler法：法：110 10 950 1051 05( .).nnnyxy Euler法取步長(zhǎng)取步長(zhǎng) ，這時(shí)，這時(shí)

17、 ，計(jì)算到，計(jì)算到 ，誤差，誤差 . . 計(jì)算結(jié)果及誤差分別為計(jì)算結(jié)果及誤差分別為: : 0 1 .h 0 1 .nxn 5n ()nnny xy EulerEuler法法后退后退EulerEuler法法梯形法梯形法改進(jìn)的改進(jìn)的EulerEuler法法準(zhǔn)確值準(zhǔn)確值0.10.11 10000000000001 10090910090911 10047620047621 10050000050001 10048370048370.20.21 10100000100001 10264460264461 10185940185941 10190250190251 10187310187310.30.3

18、1 10290000290001 10513150513151 10406330406331 10412180412181 10408180408180.40.41 10561000561001 10830130830131 10700960700961 10708020708021 10703200703200.50.51 10904900904901 11209211209211 11062781062781 11070761070761 1106531106531nxEuler法nx0.10.10.0048370.0048370.0042540.0042540.0000750.00007

19、50.0001630.0001630.20.20.0087310.0087310.0077150.0077150.0001370.0001370.0002940.0002940.30.30.0118180.0118180.0104970.0104970.0001850.0001850.0004000.0004000.40.40.0142200.0142200.0126930.0126930.0002240.0002240.0004820.0004820.50.50.0160410.0160410.0143900.0143900.0002530.0002530.0005450.000545|n

20、|n |n |n Euler法1.1.梯形法和改進(jìn)梯形法和改進(jìn)Euler Euler 法的誤差較小，法的誤差較小，EulerEuler法和和后退法和和后退EulerEuler法誤差較大。法誤差較大。注：注：2.2.梯形法計(jì)算中為什么沒(méi)有采用迭代法計(jì)算？梯形法計(jì)算中為什么沒(méi)有采用迭代法計(jì)算？3.3.從理論上如何刻畫(huà)某種方法的誤差？從理論上如何刻畫(huà)某種方法的誤差？Euler法11(, )nnnnnyyhxyyh 1(, )nnnnyyhxyh 隱式隱式單步法單步法通常稱通常稱 為為增量增量函數(shù)函數(shù)( , , )x y h 顯式顯式單步法單步法定義定義2 若計(jì)算若計(jì)算 時(shí)只用到前一點(diǎn)的值時(shí)只用到前一

21、點(diǎn)的值 ，這類算法稱為單，這類算法稱為單步法步法. .若計(jì)算若計(jì)算 時(shí)需要用到時(shí)需要用到 前面前面 點(diǎn)的值點(diǎn)的值 ，這類算法稱，這類算法稱為為 步法步法. .1ny ny1ny 1ny kk單步法的一般形式：?jiǎn)尾椒ǖ囊话阈问剑簡(jiǎn)尾椒椒ǖ挠嘘P(guān)概念設(shè)設(shè) 是是準(zhǔn)確準(zhǔn)確的，用某種方法計(jì)算的，用某種方法計(jì)算 時(shí)產(chǎn)生的截時(shí)產(chǎn)生的截ny稱稱 為某方法在點(diǎn)為某方法在點(diǎn) 的的整體截?cái)嗾w截?cái)嗾`差誤差()nnney xy nx1ny 斷誤差，稱為該方法的斷誤差，稱為該方法的局部截?cái)嗑植拷財(cái)嗾`差，即誤差，即111()nnnTy xy1()(, )nnnny xyhxyh 1()()(, (), )nnnny xy

22、xhxy xh 注：注： 是準(zhǔn)確的，即是準(zhǔn)確的，即ny()nnyy x 單步方法的有關(guān)概念定義定義3其中其中 為自然數(shù)，則稱該方法是為自然數(shù)，則稱該方法是 階的或具有階的或具有 階精度。階精度。 定義定義4 如果給定方法的如果給定方法的局部截?cái)嗑植拷財(cái)嗾`差為誤差為如果一個(gè)如果一個(gè) 階單步方法的階單步方法的局部截?cái)嗑植拷財(cái)嗾`差為誤差為11()pnTO h pppp121(, ()()ppnnnTg xy xhO h 則稱則稱 為該方法的局部截?cái)嗾`差的為該方法的局部截?cái)嗾`差的主項(xiàng)主項(xiàng)。1(, ()pnng xy xh 如向前如向前EulerEuler方法的方法的局部截?cái)嗑植拷財(cái)嗾`差誤差11()()

23、(, ()nnnnnTy xy xhf xy x22()()()!nnnhhy xyxhy x2323()()!nnhhyxyx2()O h 一階一階方法方法單步方法的有關(guān)概念如向后如向后EulerEuler方法的方法的局部截?cái)嗑植拷財(cái)嗾`差誤差1111()()(, ()nnnnnTy xy xhf xy x232()()!nhyxO h 2()O h 一階一階方法方法 上述局部截?cái)嗾`差的定義對(duì)上述局部截?cái)嗾`差的定義對(duì)隱式單步法隱式單步法也適用。也適用。如梯形如梯形法的法的局部截?cái)嗑植拷財(cái)嗾`差誤差1111()()(, ()(, ()2nnnnnnnhTy xy xf xy xf xy x 341

24、2()()nhyxO h 3()O h 二階二階方法方法單步方法的有關(guān)概念定義定義5 對(duì)于單步法，它在對(duì)于單步法，它在 處的解為處的解為 ，若對(duì)任意，若對(duì)任意固定的固定的 ，均有，均有 ，稱，稱單步法收斂單步法收斂。0nxxnhnynx0lim()nnhyy x 定理定理2 設(shè)單步法具有設(shè)單步法具有 階精度，增量函數(shù)階精度，增量函數(shù) 關(guān)于關(guān)于 滿足滿足Lipschitz條件，即存在常數(shù)條件，即存在常數(shù) ，使得對(duì)任意，使得對(duì)任意 成立成立.又假設(shè)初值是準(zhǔn)確的，即又假設(shè)初值是準(zhǔn)確的，即 .則則單步單步法收斂且法收斂且 1()p ( , , )x y h yL 12,y y1212|( , )( ,

25、 )|x y hx y hLyy 00()yy x ()()pnny xyO h上述定義和定理是針對(duì)顯式單步法，隱式可類似定義，上述定義和定理是針對(duì)顯式單步法，隱式可類似定義，有類似結(jié)果。有類似結(jié)果。單步方法的有關(guān)概念Runge-KuttaRunge-Kutta方法的方法的基本思想基本思想顯式顯式單步法的一般形式：?jiǎn)尾椒ǖ囊话阈问剑?(, )nnnnyyhxyf h R-K方法是利用一些點(diǎn)的線性組合構(gòu)造方法是利用一些點(diǎn)的線性組合構(gòu)造增量函數(shù)增量函數(shù)，使得相應(yīng)方法的使得相應(yīng)方法的局部局部截?cái)嗾`差的截?cái)嗾`差的階數(shù)階數(shù)盡可能盡可能高高。二階二階Runge-KuttaRunge-Kutta方法方法12

26、221( , , , )( , )(,( , )x y f hc f x yc f xa h yb hf x y 確定參數(shù)確定參數(shù) ，使得，使得12221,c c a b與與 在點(diǎn)在點(diǎn) 的的Taylor展開(kāi)式有盡可能多的展開(kāi)式有盡可能多的相同項(xiàng)相同項(xiàng)。 ( )( , , )y xhx y f h ()y xh x龍格-庫(kù)塔（Runge-Kutta）方法12221( )( , , )( )( , )(,( , )y xh c f x yc f xa h yy xhx y fb hfhx y 122221( )( , ) ( , )( , )( , )( , )()xyy xh c f x ycf

27、 x ya hfx yb hf x y fx yO h 232( )( )()(hy xhy xyxxOhyh 22( ) ( , )( , )( , ) ( , )( )xyhy xh f x yf x yf x y f x yOh 比較兩式的比較兩式的相同項(xiàng)相同項(xiàng)得得2 2112c b 121cc 2212c a 方程組有無(wú)窮多解方程組有無(wú)窮多解龍格-庫(kù)塔（Runge-Kutta）方法若取其一組解若取其一組解122211112,ccab 21(,)nnKf xh yhK 1122()nnhyyKK 1()nnKf x y 則得到則得到改進(jìn)改進(jìn)的的EulerEuler公式（公式（二階二階方法

28、）方法）若取其另一組解若取其另一組解1222110,1,2ccab則得到則得到二階二階的中點(diǎn)公式。的中點(diǎn)公式。二階龍格-庫(kù)塔方法定義定義5 設(shè)設(shè) 是一個(gè)正整數(shù)，代表使用函數(shù)值是一個(gè)正整數(shù)，代表使用函數(shù)值 的個(gè)數(shù)，的個(gè)數(shù)，和和fr是一些特定的權(quán)因子（均為是一些特定的權(quán)因子（均為2 31 21,(, , ;, ,iija bir ji 1 2(, , )ic ir 實(shí)數(shù)），則稱下列方法（公式）實(shí)數(shù)），則稱下列方法（公式）111()nnrryyh c Kc K 為初值問(wèn)題為初值問(wèn)題 的的r r級(jí)級(jí)顯式顯式RungeKuttaRungeKutta公式公式. .其中其中11(,)rrnrnriiiKf

29、xa h yhb K 1(),nnKf x y 22211(,),nnKf xa h yhb K 顯示龍格-庫(kù)塔方法類似前面的處理方法，可以得到類似前面的處理方法，可以得到四級(jí)四級(jí)方法：方法：r =45()O h局部截?cái)嗑植拷財(cái)嗾`差誤差最常用的一種最常用的一種四階四階方法：經(jīng)典顯式方法：經(jīng)典顯式Runge-Kutta公式公式11234226()nnhyykkkk 1(,)nnkf xy 2122(,)nnhhkf xyk 3222(,)nnhhkf xyk 43(,)nnkf xh yhk 四階Runge-Kutta方法解101( )dyxydxy 0 1 , x 例例4 用經(jīng)典的四階用經(jīng)典的

30、四階Runge-Kutta方法求解下列初值問(wèn)題方法求解下列初值問(wèn)題 。0 1 .h 經(jīng)典的四階經(jīng)典的四階Runge-KuttaRunge-Kutta公式：公式：11234226()nnhyykkkk 11;nnkxy 21122nnhhkxyk 431.nnkxhyhk32122;nnhhkxyk 四階Runge-Kutta方法四階Runge-Kutta方法四階Runge-Kutta方法解向前向前EulerEuler法步長(zhǎng)取法步長(zhǎng)取100( )dyydxy 例例5：分別利用分別利用三種方法三種方法求解初值問(wèn)題的求解初值問(wèn)題的數(shù)值解數(shù)值解0 025.h 改進(jìn)改進(jìn)EulerEuler法步長(zhǎng)取法步長(zhǎng)

31、取0 05.h 經(jīng)典經(jīng)典4 4階階RKRK法步長(zhǎng)取法步長(zhǎng)取0 1 .h 計(jì)算結(jié)果為：計(jì)算結(jié)果為： 常微分方程數(shù)值解法EulerEuler方法方法改進(jìn)的改進(jìn)的EulerEuler法法經(jīng)典的經(jīng)典的R RK K方法方法準(zhǔn)確值準(zhǔn)確值 000000.10.0963120.0951230.095162500.095162580.20.1833480.1811930.181269100.181269250.30.2620010.2590850.259181580.259181780.40.3330790.3295630.329679710.329679950.50.3973120.3933370.39346

32、9060.39346934nx常微分方程數(shù)值解法注：注：1. 三種方法的步長(zhǎng)取法有何特點(diǎn)？三種方法的步長(zhǎng)取法有何特點(diǎn)？2. 比較三種方法在節(jié)點(diǎn)處的誤差？比較三種方法在節(jié)點(diǎn)處的誤差？3. 如在節(jié)點(diǎn)如在節(jié)點(diǎn)0.5處，三種方法的誤差分別為處，三種方法的誤差分別為3473.8 101.3 102.8 10常微分方程數(shù)值解法k步步線性多步法線性多步法 線性多步法 所謂的線性所謂的線性多步法多步法，指的是某一步解的公式不僅與前一，指的是某一步解的公式不僅與前一步的值有關(guān)，而且與前面若干步解的值有關(guān)的方法。步的值有關(guān)，而且與前面若干步解的值有關(guān)的方法。 利用前面多步的信息，則可以期望得到較高的精度。利用前面

33、多步的信息，則可以期望得到較高的精度。 構(gòu)造多步法的主要途徑是基于構(gòu)造多步法的主要途徑是基于數(shù)值積分的方法和基于數(shù)值積分的方法和基于TaylorTaylor展開(kāi)法展開(kāi)法. .前者是將常微分方程兩端積分后利用插值求前者是將常微分方程兩端積分后利用插值求積公式得到，后者是利用局部截?cái)嗾`差定義和積公式得到，后者是利用局部截?cái)嗾`差定義和TaylorTaylor展開(kāi)展開(kāi)得到得到。Linear Mutistep Methodk步線性多步法的步線性多步法的一般形式一般形式100kkn kin iin iiiyyhf 其中其中 為為 的的近似近似， ， 為常數(shù)，為常數(shù)， 不全為零不全為零. .由于由于上式上式

34、給出了給出了 之間之間的線性關(guān)系，的線性關(guān)系，故稱為故稱為線性線性k k步法步法. .n iy ()n iy x (,)n in in iff xy ,ii 00, ,n in iyf ,稱為稱為顯式顯式k步法步法， 稱為稱為隱式隱式k步法。步法。0k 0k 可根據(jù)可根據(jù)局部截?cái)嗾`差以及階確定。局部截?cái)嗾`差以及階確定。,ii 線性多步法 將將 在在 處進(jìn)行處進(jìn)行Taylor展開(kāi)展開(kāi)n kT 100()()()kkn kn kin iin iiiTy xy xhy x 在在 處處局部截?cái)嗾`差局部截?cái)嗾`差n kx nx()()n iny xy xih 23()()()()()()23!nnnnih

35、ihy xihy xyxyx 線性多步法()()n iny xy xih 2()()()()2nnnihy xihyxyx 代入局部截?cái)嗾`差的表達(dá)式代入局部截?cái)嗾`差的表達(dá)式2012()()()n knnnTc y xc hy xc h yx ()()pppnc h yx 其中其中 線性多步法00111 ()kc 11212(1)kckk 11121(2)(1)!rrkkr 01()k121121()!rrrrkckkr 則有則有 選擇選擇 滿足滿足,ii 0110,0ppcccc 1121()()()pppn kpnTchyxO h 從而該多步法是從而該多步法是p階方法階方法 線性多步法基于基

36、于數(shù)值積分?jǐn)?shù)值積分方法方法積分可得積分可得將原初值問(wèn)題中的方程在區(qū)間將原初值問(wèn)題中的方程在區(qū)間 上積分上積分4,nnxx 44()()( , ( )nnxnnxy xy xf x y x dx 被積函數(shù)用在三個(gè)節(jié)點(diǎn)被積函數(shù)用在三個(gè)節(jié)點(diǎn) 上的上的Lagrange插值代替插值代替123,nnnxxx 21011( )njijninjj ixxl xxx 線性多步法421122334( )2 (, ()(, ()2 (, ()3nnxnnnnnnxL x dxhf xy xf xy xf xy x 稱為米爾尼（稱為米爾尼（Milne）方法：）方法：顯式四步四階方法顯式四步四階方法被積函數(shù)用在三個(gè)節(jié)點(diǎn)

37、被積函數(shù)用在三個(gè)節(jié)點(diǎn) 上的上的Lagrange插值插值代替代替123,nnnxxx 用用 表示表示 的近似值，記的近似值，記 ，從而得到，從而得到j(luò)y()jy x(,)jjjff xy 41234223()nnnnnhyyfff 線性多步法基于基于數(shù)值積分?jǐn)?shù)值積分方法方法稱為稱為Simpson方法：方法：隱式兩步四階方法隱式兩步四階方法將原初值問(wèn)題中的方程在區(qū)間將原初值問(wèn)題中的方程在區(qū)間 上積分上積分2,nnxx 22()()( , ( )nnxnnxy xy xf x y x dx 被積函數(shù)用在三個(gè)節(jié)點(diǎn)被積函數(shù)用在三個(gè)節(jié)點(diǎn) 上的上的Lagrange插值代替，計(jì)算后可得插值代替，計(jì)算后可得12

38、,nnnxxx 212()3nnnnnhyyfff 線性多步法基于基于Taylor展開(kāi)展開(kāi)方法方法考慮形如考慮形如的的K步法，稱為阿當(dāng)姆斯（步法，稱為阿當(dāng)姆斯（Adams）方法）方法10kn kn kin iiyyhf 0k 為顯式方法，為顯式方法， 為隱式方法為隱式方法0k 線性多步法K=3時(shí)，時(shí)，Adams顯式三步三階顯式三步三階方法方法3221(23165)12nnnnnhyyfff K=3時(shí)，時(shí)，Adams隱式三步四階隱式三步四階方法方法32321919524()nnnnnnhyyffff 43321555937924()nnnnnnhyyffff K=4時(shí)，時(shí)，Adams顯式四步四階

39、顯式四步四階方法方法 線性多步法利用前面所述的利用前面所述的Taylor展開(kāi)方法可得：展開(kāi)方法可得：基于基于Taylor展開(kāi)方法構(gòu)造線性多步法比較靈活，可以構(gòu)造任意展開(kāi)方法構(gòu)造線性多步法比較靈活，可以構(gòu)造任意多步法公式。多步法公式。例例5 構(gòu)造初值問(wèn)題構(gòu)造初值問(wèn)題試確定參數(shù)試確定參數(shù) 使方法的使方法的局部截?cái)嗾`差局部截?cái)嗾`差的階的階盡可能盡可能高高，并求局部截?cái)嗾`差。，并求局部截?cái)嗾`差。,ii 的顯式二步公式的顯式二步公式00( , )()yf x yy xy 1011011()nnnnnyyyhff 線性多步法解 由局部截?cái)嗾`差定義并利用由局部截?cái)嗾`差定義并利用TaylorTaylor公式得

40、到公式得到101()()()nnnnTy xhy xy xh 01()()nnhy xy xh0110111() ()()()nny xhy x 23111111111()()()22662nh yxh 4451111124246( )()()()()nnyxh yxO h 線性多步法令令0110 10110 1111022 111110662得到得到04 15 04 12 局部截?cái)嗾`差為局部截?cái)嗾`差為4(4)511()()6nnTh yxO h 線性多步法從而所求二步法為從而所求二步法為111452 (2)nnnnnyyyhff 顯式多步法計(jì)算簡(jiǎn)單，但是其顯式多步法計(jì)算簡(jiǎn)單，但是其精度及計(jì)算

41、的穩(wěn)定精度及計(jì)算的穩(wěn)定性性沒(méi)有隱式方法好沒(méi)有隱式方法好. .隱式多步法一般需采用迭代隱式多步法一般需采用迭代計(jì)算計(jì)算，計(jì)算量大，計(jì)算量大. .在實(shí)際應(yīng)用中，隱式法一般不在實(shí)際應(yīng)用中，隱式法一般不單獨(dú)使用，而是用于改善用顯式方法計(jì)算得到的單獨(dú)使用，而是用于改善用顯式方法計(jì)算得到的近似值近似值. .由顯式方法給出預(yù)測(cè)，再用隱式方法校由顯式方法給出預(yù)測(cè)，再用隱式方法校正該預(yù)測(cè)值，這樣的得到的算法稱為正該預(yù)測(cè)值，這樣的得到的算法稱為預(yù)測(cè)預(yù)測(cè)校正校正方法方法. . Adams預(yù)測(cè)-校正法 從一般情況下，預(yù)測(cè)公式與校正公式都取同階的從一般情況下，預(yù)測(cè)公式與校正公式都取同階的顯式方法和隱式方法相顯式方法和隱式方法相匹配匹配。0433322555924( )(,)(,)nnnnnnhyyf xyf xy 043443391924( )(,)(,)nnnnnnhyyf xyf xy 11379(,)(,)nnnnf xyf xy 22115 (,)(,)nnnnf xyf xy 顯式顯式 隱式隱式初始迭代值由初始迭代值由4 4階階R-KR-K方法計(jì)算方法計(jì)算 Adams預(yù)測(cè)-校正法201( )dyxydxyy 0 1 , x 例例5 用用Adams預(yù)報(bào)預(yù)報(bào)- -校正公式校正公式求解下