計算方法 81 - 83 Euler格式和Runge-Kutta格式

1、則初值問題則初值問題( (* *) )的連續(xù)可微解的連續(xù)可微解y(x)存在存在, ,且唯一且唯一. . 第八章第八章 常微分方程初值問題的數(shù)值解法常微分方程初值問題的數(shù)值解法討論常微分方程（討論常微分方程（ordinary differential equation）的定解問題，這類問題）的定解問題，這類問題的最簡單型是如下一階方程的初值問題：的最簡單型是如下一階方程的初值問題：00( , ), (), .yf x yaxby xyy R ,: ybxaD12,y y, ),(),(2121yyLyxfyxf 滿足滿足Lipschitz (Lipschitz (李普希茲李普希茲) )條件，即存

2、在正數(shù)條件，即存在正數(shù)L L，使得對所有的，使得對所有的和任何和任何 均有均有定理定理 設函數(shù)設函數(shù)f(x,y)在區(qū)域在區(qū)域 上連續(xù)上連續(xù),且在區(qū)域且在區(qū)域D內(nèi)關于內(nèi)關于y由于數(shù)值解是找精確解由于數(shù)值解是找精確解 的近似值，因此，的近似值，因此，總假設總假設方程的解方程的解 存在且唯一，并具有充分的光滑性！存在且唯一，并具有充分的光滑性！( )y x( )y x , xa b(*) 求常微分方程數(shù)值解的必要性求常微分方程數(shù)值解的必要性1 1、方程本身很復雜，不能給出解析解，或難以求出解析解；、方程本身很復雜，不能給出解析解，或難以求出解析解；2 2、即使可以獲得解析解，計算量太大或者計算過程太

3、復雜而不實用；、即使可以獲得解析解，計算量太大或者計算過程太復雜而不實用；例例如：高階常系數(shù)線性常微分方程如：高階常系數(shù)線性常微分方程. .3 3、在實際應用中，只需求得解在某些特殊點上的近似值。、在實際應用中，只需求得解在某些特殊點上的近似值。數(shù)值解數(shù)值解是指：在是指：在解的存在區(qū)間解的存在區(qū)間上取一系列離散節(jié)點上取一系列離散節(jié)點012nxxxx0 1 2(, , ,)iy i 逐個給出逐個給出精確解精確解 的近似值的近似值()iy x 定義定義：相鄰兩個節(jié)點的間距相鄰兩個節(jié)點的間距 稱為稱為步長步長.1,nnnhxx 如何求出？如何求出？考慮考慮等距等距節(jié)點：節(jié)點：0 0 1, ,nxxn

4、h n 從初始條件從初始條件 出發(fā)，按照節(jié)點的排序，依次逐個計算出發(fā)，按照節(jié)點的排序，依次逐個計算 的的 值，稱為值，稱為步進法步進法. 一般有兩種類型：一般有兩種類型：單步法單步法、多步法多步法.00( )y xy iy注：注：bahn 為了考察數(shù)值解是否具有使用價值，必須解決的基本問題：為了考察數(shù)值解是否具有使用價值，必須解決的基本問題：解析解不能用初等函數(shù)及其積分表示！解析解不能用初等函數(shù)及其積分表示！當步長當步長h h取得充分小時，取得充分小時，所得的近似解所得的近似解y yn n能否以足夠的能否以足夠的精度精度逼近逼近初值初值問題的精確解問題的精確解y(y(x xn n) ). .這

5、就是這就是收斂問題收斂問題。即當。即當h 0h 0時時, , y yn n y y( (x xn n) ?) ?在數(shù)值求解的過程中，會產(chǎn)生若干類型的誤差，具體分類如下：在數(shù)值求解的過程中，會產(chǎn)生若干類型的誤差，具體分類如下： （1）局部截斷誤差局部截斷誤差 （2 2）局部舍入誤差）局部舍入誤差 （3 3）整體截斷誤差）整體截斷誤差 （4 4）整體舍入誤差）整體舍入誤差 （5 5）總誤差）總誤差（= =整體截斷誤差整體截斷誤差+ +整體舍入誤差）整體舍入誤差）因此，必須因此，必須估計估計精確解與近似解之間的精確解與近似解之間的誤差誤差。這就是這就是誤差估計問題誤差估計問題。由于初始值由于初始值y

6、 y0 0和右端項和右端項f(f(x,yx,y) )常常是通過測量得到的，所以必須考慮常常是通過測量得到的，所以必須考慮 它們的微小擾動，引起數(shù)值解的變化問題。即最初產(chǎn)生的誤差在以它們的微小擾動，引起數(shù)值解的變化問題。即最初產(chǎn)生的誤差在以 后各步的計算中是否會無限制擴大的問題，這就是后各步的計算中是否會無限制擴大的問題，這就是穩(wěn)定性問題穩(wěn)定性問題。在計算過程中無舍入誤差，只有當問題的數(shù)值解對初始值具有某種連續(xù)在計算過程中無舍入誤差，只有當問題的數(shù)值解對初始值具有某種連續(xù)依賴性時，方法才實用！依賴性時，方法才實用！可證明：可證明：若不考慮初始值誤差，整體截若不考慮初始值誤差，整體截斷誤差的階由局

7、部截斷誤差的階決定！斷誤差的階由局部截斷誤差的階決定！本章著重討論一階本章著重討論一階ODE初值問題的數(shù)值解初值問題的數(shù)值解.對于高階方程（組）對于高階方程（組）的數(shù)值解，其基本思想是完全一樣的的數(shù)值解，其基本思想是完全一樣的.計算一步所產(chǎn)生的誤差。是算法中所固有的，與舍入誤差無關計算一步所產(chǎn)生的誤差。是算法中所固有的，與舍入誤差無關初值問題的初值問題的解析解解析解表示過點表示過點 的一條（的一條（光滑光滑）曲線曲線. .00(,)x y 解析解與其數(shù)值解的解析解與其數(shù)值解的幾何意義幾何意義：oxynx 0 x0 00 0( (x x, ,y y ) ) 1x 初值問題的初值問題的數(shù)值解數(shù)值解

8、表示一組表示一組離散點列離散點列 （或一組數(shù)據(jù)點）（或一組數(shù)據(jù)點）(,)iix y可用可用擬合擬合方法求該組數(shù)據(jù)方法求該組數(shù)據(jù) 的的近似曲線近似曲線( , )iix y2x11(,)x y積分曲線積分曲線近似曲線近似曲線( , )yf x y (,)nnxy22(,)xy本章給出的幾種方法本章給出的幾種方法一、歐拉一、歐拉(Euler)(Euler)方法及其改進形式方法及其改進形式二、龍格二、龍格- -庫塔庫塔( (Runge-KuttaRunge-Kutta) )方法方法三、線性多步法三、線性多步法-Adams-Adams（埃德姆斯）方法（埃德姆斯）方法 得到得到數(shù)值解數(shù)值解有兩個基本途徑有

9、兩個基本途徑： 把近似解表示成有限個獨立函數(shù)之和把近似解表示成有限個獨立函數(shù)之和，例如：截斷的冪（，例如：截斷的冪（TaylorTaylor）級數(shù)或正交函數(shù)展開式中的前幾項級數(shù)或正交函數(shù)展開式中的前幾項. . 涉及到計算高階導，盡管可用涉及到計算高階導，盡管可用“數(shù)數(shù)值微分值微分”技術，但得到的公式太長、太復雜！通常比較適用于手算技術，但得到的公式太長、太復雜！通常比較適用于手算. . 離散化方法離散化方法（也稱為（也稱為差分方法差分方法），它），它提供了提供了用當前節(jié)點上的或前幾用當前節(jié)點上的或前幾個節(jié)點上的近似值來個節(jié)點上的近似值來計算計算下一個節(jié)點上的近似值的下一個節(jié)點上的近似值的規(guī)則規(guī)

10、則. .數(shù)值解所滿足的數(shù)值解所滿足的離散方程統(tǒng)稱為離散方程統(tǒng)稱為差分格式差分格式. . 它是本章中要研究的一種方法它是本章中要研究的一種方法. .微微分分方方程程區(qū)域剖分區(qū)域剖分遞推計算或解線性遞推計算或解線性代數(shù)方程組代數(shù)方程組微分方程離散微分方程離散初始和邊界條件處理初始和邊界條件處理解的存在性、唯一性解的存在性、唯一性解的收斂性和收斂速度解的收斂性和收斂速度解的穩(wěn)定性解的穩(wěn)定性得到數(shù)值解得到數(shù)值解離散系統(tǒng)的離散系統(tǒng)的性態(tài)研究性態(tài)研究現(xiàn)實問題現(xiàn)實問題數(shù)學模型數(shù)學模型離散格式離散格式模型模型誤差誤差舍入舍入誤差誤差觀測觀測誤差誤差截斷截斷誤差誤差數(shù)值解數(shù)值解數(shù)學模型數(shù)學模型離散化離散化計算計

11、算第第2 2節(jié)節(jié) 歐拉（歐拉（EulerEuler）方法）方法 最簡單的一種方法，精度差，不推薦使用！最簡單的一種方法，精度差，不推薦使用！ 歐拉格式的構造歐拉格式的構造解決問題的關鍵：解決問題的關鍵：如何處理方程中的導數(shù)項？如何處理方程中的導數(shù)項？00()(, () ,()nnny xf x y xy xy 在各節(jié)點在各節(jié)點 處處, ,有有求求 的近似值的近似值()ny xny方法：方法：將上述初值問題將上述初值問題化成化成節(jié)點離散方程，在節(jié)點上采用節(jié)點離散方程，在節(jié)點上采用離散化方法離散化方法（也也叫做差分方法叫做差分方法, ,通常用數(shù)值積分、微分、泰勒公式等）通常用數(shù)值積分、微分、泰勒公

12、式等），可逼近，可逼近節(jié)點離散方節(jié)點離散方程程，由此產(chǎn)生，由此產(chǎn)生可計算格式可計算格式，并用計算解，并用計算解 作為解析解作為解析解 的近似值的近似值. .iy( )iy x在在節(jié)點離散方程節(jié)點離散方程中直接用中直接用向前差商向前差商代替微商，得到代替微商，得到nx節(jié)點離散方程節(jié)點離散方程11()()(, ()nnnnnny xy xf x y xxx 1()()(, ()nnnny xy xhf x y x 101(,), , ,nnnnyyhf x yn 稱為稱為歐拉格式歐拉格式.當初值當初值 已知時，已知時，0y可遞推求出可遞推求出12,.,.nyyy( )iiyy x 令令bahn 該

13、切線與直線該切線與直線 的的交點交點的縱坐標的縱坐標現(xiàn)在從現(xiàn)在從 出發(fā)出發(fā)00(,)xyoxy),(00yx0 x),(00yx nx (,)nnxy 2x),(22yx 1x),(11yx 作解曲線的切線，作解曲線的切線，切線方程為：切線方程為：0000(,)()yyf xyxx 切線切線斜率斜率為為1y100010(,)()yyf xyxx 11( ,)x y1111(,)()yyf x yxx 歐拉方法的幾何解釋歐拉方法的幾何解釋00(,),f x y1xx 11( )y xy 再從再從 出發(fā)，以出發(fā)，以 為斜率作直線為斜率作直線11( ,)f x y該直線與該直線與直線直線 的交點的的

14、交點的縱坐標縱坐標2xx 2y211121(,)()yyf xyxx22()y xy 依次類推，依次類推， 的近似值的近似值 的計算公式：的計算公式：1()ny x 1ny 1(,)nnnnyyhf x y 即是即是歐拉公式歐拉公式.又稱又稱用一條折線近似代替積分曲線折線法折線法.歐拉方法歐拉方法的缺陷：誤差比較大！的缺陷：誤差比較大！( )y y x 令令令令( , )yf x y 方法可看做是方法可看做是TaylorTaylor級數(shù)前兩項的近似！級數(shù)前兩項的近似！y(x1)y(x2)y2y(x3)y3y1x1x2x3xn.x0y(x0)EulerEuler格式的整體截斷誤差與局部截斷誤差，

15、示意圖格式的整體截斷誤差與局部截斷誤差，示意圖整體截斷誤差局部截斷誤差在在節(jié)點離散方程節(jié)點離散方程中直接用中直接用向后差商向后差商代替微商，得到代替微商，得到11()()(, ()nnnnnny xy xf x y xxx 1()()(, ()nnnny xy xhf x y x 101(,), , ,nnnnyyhf x yn 稱為稱為隱式歐拉格式隱式歐拉格式. .在在節(jié)點離散方程節(jié)點離散方程中用中用中心差商中心差商代替微商，得到代替微商，得到1111()()(, ()nnnnnny xy xf x y xxx 112()()( , ( )nnnny xy xhf x y x 11201(,

16、), , ,nnnnyyhf x yn 稱為稱為歐拉兩步格式歐拉兩步格式.于是于是11101(,), , ,nnnnyyhf xyn 關于關于 的非線性方程，可用的非線性方程，可用迭代法求解！迭代法求解！1ny 歐拉格式的歐拉格式的局部截斷誤差局部截斷誤差假設假設 ，即第，即第n步的結果步的結果 是準確的，是準確的， 稱為稱為局部截斷誤差局部截斷誤差.()nnyy x ny11()nnyy x 定義定義2.1. 2.1. 如果一種方法的局部截斷誤差為如果一種方法的局部截斷誤差為 ，則稱該數(shù)值方法的精度，則稱該數(shù)值方法的精度 是是 階階的，或簡稱該方法是的，或簡稱該方法是 階階的的. .1()p

17、O h p11(,) ()(, ()() ()nnnnnnnnnnfyxy xyhf xyyhy xy xy xh 顯式顯式EulerEuler格式格式132112( )()(Tayl)()+or()nnnnny xy xhyxO hhxxy 根根據(jù)據(jù)公公式式，有有，13221112()( ) ( ) ()nnnnyy xhO hO hxxy 因因此此有有，計算一步所產(chǎn)生的誤差計算一步所產(chǎn)生的誤差一階方法一階方法( , )yf x y 局部截斷誤差的主項局部截斷誤差的主項p1111111(,) (,) ()()()nnnnnnnnnnf xyyhf xyyhyyxyhxxy 隱式隱式Euler

18、Euler格式格式Taylor根根據(jù)據(jù)公公式式，有有 1111232312111223()()()() ()+()()! =()nnnnnnyy xy xhyxy xh yxh yhOyh 因因此此有有1232211132()()+()!()nnnnnnyy xhyxhxxxxh yy ，2111112()+()()()nnnnnyxh yyxhyxxx ，111(),()(,) nnnnnyy xyxf xy 令令，則則一階方法一階方法 Euler兩步格式兩步格式11111222(,) (, () )()()nnnnnnnnnnf xy xyhf xyyhy xhxyyyx Taylor根根

19、據(jù)據(jù)公公式式，有有3343111133()()()() ! (=)nnnnh yxh yxO hOyy xh 因因此此有有13241213()()+()()!()()nnnnnhy xhyxh yxO hxy xy 13421312()+()()(!)(nnnnnh yy xhyxxy xh yxO h 11(),( )nnnnyy xyy x 令令，則則二階方法二階方法在區(qū)間在區(qū)間 上對微分方程上對微分方程 積分得積分得 歐拉格式的歐拉格式的積分學解釋積分學解釋( )( , ( )y xf x y x 1,kkxx 11( ) d( , ( ) dkkkkxxxxy xxf x y xx 即

20、即11()()( , ( ) dkkxkkxy xy xf x y xx 為了給出迭代格式，只要對積分項提供一種算法！而選擇不同的算法，就為了給出迭代格式，只要對積分項提供一種算法！而選擇不同的算法，就會得到不同的迭代格式會得到不同的迭代格式. . 利用左矩形公式，有利用左矩形公式，有21()()(, ()()kkkky xy xhf xy xO h 略去高階項，便得略去高階項，便得（顯式）歐拉公式（顯式）歐拉公式. 利用右矩形公式，有利用右矩形公式，有2111()()(, ()()kkkky xy xhf xy xO h 略去高階項，便得略去高階項，便得（隱式）歐拉公式（隱式）歐拉公式.略去

21、高階項，然后用略去高階項，然后用 代替代替 ，得，得利用梯形公式，有利用梯形公式，有11()()( , )kkxkkxy xy xf x y dx 1132(, ()(, ()()kkkkhfO hxy xf xy x()iy xiy 1112(,)(,)kkkkkkhyyf xyf xy 單步隱單步隱格式格式稱為稱為梯形格式梯形格式. .二階方法二階方法關于關于 的非線性方程，可用的非線性方程，可用簡單迭代法求解！簡單迭代法求解！1ky 例如：例如： 可構造可構造不動點迭代格式不動點迭代格式1111 20 1()( )(,),(,) ,nnnnnnkkhyyf xyf xky 其中迭代函數(shù)為

22、其中迭代函數(shù)為 1112(,)()(,)nnnnnnhyf xfyxyy 為了保證迭代法收斂，我們分析為了保證迭代法收斂，我們分析111111111112( )()( )()( )|(,)(,)| |()()|nnnnnnnnkkkkkkyyhf xf xyyyy 1112( )()|nnkkL yyh 于是，對任意正整數(shù)于是，對任意正整數(shù)p，有，有111111111121()( )()()()()()| |() ()(| )nnnnnnnnk pk pkk pkk pk pkyyyyyyyy 1111111121()()()()()( )|nnnnnnk pk pkk ppkkyyyyyy

23、1111111112112111221222( )()( )()()( )()( )|nnnnnnnnkkkkkppkkppkhhLyLyhhhLLyyyyyLy 111112()( )|nnkkyhyL 11101122( )( )|nnkLyyhhL 012hL 只只要要這表明：這表明：序列序列 是是CauchyCauchy列！列！ 10( )nkkky 將將顯式歐拉公式顯式歐拉公式和和梯形公式梯形公式聯(lián)合使用聯(lián)合使用11112 (,)(,)(,)nnnnnnnnnnyhf xyhyyyyf xyf x 即先用歐拉格式即先用歐拉格式預估預估一個粗糙的近似值，然后用梯形公式對其進行一個粗糙的

24、近似值，然后用梯形公式對其進行校正校正.這種這種預估預估-校正格式校正格式稱為稱為改進的改進的Euler格式格式.它可表為如下等價形式：它可表為如下等價形式：12112()nnkykyh 1(,)nnf xyk 0 1, ,n 12(,)nnf xh ykkh 單步顯單步顯格式格式是是 的斜率值的斜率值.nx是是 的斜率值，它是利用的斜率值，它是利用已知信息已知信息 通過通過Euler格式格式得到的得到的.1nx ny EulerEuler格式的格式的微分學解釋微分學解釋依據(jù)微分中值定理，依據(jù)微分中值定理，有有11 ()( ),( ,().( )nnnny xy xhxxhyfy 1*( ,

25、( ),( ),nnKfyy xxx 稱稱之之為為在在區(qū)區(qū)記記平平間間上上的的均均斜斜率率. .是是 兩點的斜率值兩點的斜率值 的算術平均作為的算術平均作為平均平均斜率斜率 的近似值的近似值.1,nnxx 12,kk*K局部截斷誤差局部截斷誤差為為O( (h3) )二階方法二階方法一般地，我們考慮：一般地，我們考慮：如果設法在如果設法在 內(nèi)多預報幾個點的斜率值，然內(nèi)多預報幾個點的斜率值，然后將它們加權平均作為后將它們加權平均作為平均斜率平均斜率 的近似值，則有可能構造出更高精度的近似值，則有可能構造出更高精度的計算格式的計算格式 這就是這就是Runge-KuttaRunge-Kutta方法的基

26、本思想方法的基本思想. .1,nnxx *K第第3 3節(jié)節(jié) Runge-KuttaRunge-Kutta方法方法 - - 一種提高截斷誤差階的方法一種提高截斷誤差階的方法 二階二階Runge-KuttaRunge-Kutta格式格式121122111221 01 *,()n pnnn pnnnnxxphpxxKKxKKKyyhxKK 在在區(qū)區(qū)間間內(nèi)內(nèi)選選取取一一點點并并用用兩兩個個點點的的斜斜率率值值的的加加權權平平均均平平均均斜斜率率的的 近近似似值值，使使所所構構造造出出的的計計算算格格式式問問題題的的精精確確解解問問，題題且且 與與與與作作為為逼逼近近3 3具具有有二二階階精精度度. .

27、12(),(,)：n pn pnn pn py xyyhKKf xyp 改改進進的的E Eu ul le er r格格式式仿仿照照，用用E Eu ul le er r格格式式預預測測的的值值并并用用它它來來估估計計斜斜率率，于于是是得得到到如如下下計計算算格格式式：11211221 ()(,),(,).nnnnn pnyyhKKKf xyKf xyphK 1)()n pnn pny xy xKxx 問題：問題：適當選取待定適當選取待定參數(shù)，使得格式具有參數(shù)，使得格式具有二階精度二階精度. .工具：工具：用用TaylorTaylor公式，公式，分析格式的局部截斷分析格式的局部截斷誤差誤差nx1n

28、x n px ( , )yf x y 為此，考慮格式的局部截斷誤差為此，考慮格式的局部截斷誤差. . 1211212 =(,)(,)(,)(,)(,)(,)()(, ()()(,)()nnn pnnnnnnnnxnnynnnnnnnnKf xyKf xyhKf xph yphf xyph fxyfxyOyy xKKxy xpKf xyyxhf xy x 注注意意到到：設設，分分別別將將與與在在處處作作T Ta ay yl lo or r展展開開，有有 112122312122122 = 12=(,)()(,)()()()()()()()()()()()()nnxnnnynnnnnnnnnnnn

29、nfxyy xfxyKKy xy xphyxO hy xy xyxy xy xyxyyhy xhhhphO hy xhh 于于是是和和3)()O h 比較系數(shù)可知，若要格式具有二階精度，只需比較系數(shù)可知，若要格式具有二階精度，只需122112,.p 易易見見，滿滿足足條條件件的的參參數(shù)數(shù)不不止止一一組組，而而是是一一族族. .二二階階把把滿滿R Ru u足足這這一一條條n ng ge e- -件件的的一一族族格格式式統(tǒng)統(tǒng)稱稱K Ku ut tt t為為a a格格式式. .( , )yf x y 特別，若取特別，若取 有有 格式就是格式就是改進的改進的EulerEuler格式格式. .12211

30、2,.p 1,p 1212. 若取若取 有有 格式稱為格式稱為變形的變形的EulerEuler格式（也稱為中點格式）格式（也稱為中點格式）. .其形式為其形式為12,p 1201,. 1212112 2(,),(,).nnnnnnyyhKKf xyhKf xyK 三階三階Runge-KuttaRunge-Kutta格式格式12311223311122313 1 *,(,)，nn pn qnnn pn qnnnnxxxxqh pqxxxKKKKKKKyyxxhKKK 在在區(qū)區(qū)間間內(nèi)內(nèi)除除了了和和外外再再選選取取一一點點并并用用三三個個點點的的斜斜率率值值的的加加權權平平均均平平均均斜斜率率 ,

31、,問問題題4 4的的近近似似值值，使使所所構構造造出出的的計計算算格格式式問問題題的的精精確確解解，且且具具有有三三, ,作作為為逼逼近近階階精精度度. .112112212, ,(,)(),(,)., ,()(),：n pnnn pnnn pn pnn pn pn qnnn qn qn qnxxhx xKf x yy xyyhKKf xyxxhx xK Ky xyyhpqKpqK ，點點此此時時，在在區(qū)區(qū)間間內(nèi)內(nèi)有有一一個個斜斜率率值值用用E Eu ul le er r格格式式預預測測的的值值并并用用它它來來估估計計斜斜率率，點點此此時時，在在區(qū)區(qū)間間內(nèi)內(nèi)有有兩兩個個斜斜率率值值，同同樣樣，

32、用用E Eu ul le er r格格式式預預測測注注意意的的值值, ,到到注注意意到到, ,并并用用它它來來3(,).n qn qKf xy 估估計計斜斜率率nx1nx n px n qx 123123123,K K KKKK 利利用用這這三三個個斜斜率率值值作作加加權權平平均均得得出出上上的的平平均均斜斜率率，得得到到如如下下計計算算格格式式：123121123121312 ()(,),(,),(,().nnnnn pnn qnyyhKKKKf xyKf xyhKKf xyhKKpq 問題：問題：適當選取適當選取7 7個個待定參數(shù)，使得格式待定參數(shù)，使得格式具有三階精度具有三階精度. .類似地，利用類似地，