數(shù)值分析第五章常微分方程數(shù)值解法

1、第五章第五章 常微分方程數(shù)值解常微分方程數(shù)值解/* Numerical Methods for Ordinary Differential Equations */ 待求解的問題待求解的問題：一階一階常微分方程的常微分方程的初值問題初值問題 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy 解的存在唯一性解的存在唯一性（“常微分方程常微分方程”理論）：只要理論）：只要 f (x, y) 在在a, b R1 上連續(xù)，且關(guān)于上連續(xù)，且關(guān)于 y 滿足滿足 Lipschitz 條件條件，即存在與，即存在與 x, y 無關(guān)的常數(shù)無關(guān)的常數(shù) L 使使對任意

2、定義在對任意定義在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，則上述都成立，則上述IVP存存在唯一解在唯一解。| ),(),(|2121yyLyxfyxf 解析解法解析解法：（常微分方程理論）：（常微分方程理論）只能求解極少一類常微分方程；實際中給定的問題不一只能求解極少一類常微分方程；實際中給定的問題不一定是解析表達式，而是函數(shù)表，無法用解析解法。定是解析表達式，而是函數(shù)表，無法用解析解法。如何求解如何求解計算解函數(shù)計算解函數(shù) y(x) 在一系列節(jié)點在一系列節(jié)點 a = x0 x1 xn= b 處的近似值處的近似值),., 1()(nixyyii 節(jié)點間距節(jié)點間距 為步長，通

3、常采用為步長，通常采用等距節(jié)點等距節(jié)點，即取即取 hi = h (常數(shù)常數(shù))。) 1,., 0(1 nixxhiii數(shù)值解法數(shù)值解法: 求解所有的常微分方程求解所有的常微分方程步進式步進式：根據(jù)已知的或已求出的節(jié)點上的函數(shù)值計算：根據(jù)已知的或已求出的節(jié)點上的函數(shù)值計算當前節(jié)點上的函數(shù)值，一步一步向前推進。因此只需當前節(jié)點上的函數(shù)值，一步一步向前推進。因此只需建立由已知的或已求出的節(jié)點上的函數(shù)值求當前節(jié)點建立由已知的或已求出的節(jié)點上的函數(shù)值求當前節(jié)點函數(shù)值的遞推公式即可。函數(shù)值的遞推公式即可。111()()() ()()(,)nnnnnnnnnny xy xhy xy xyy xyyh f xy

4、1(,) 0, 1,.nnnnyyh f xyn-Eulers Method1 歐拉方法歐拉方法 /* Eulers Method */1 Eulers MethodTaylor展開法展開法幾何意義幾何意義亦稱為亦稱為歐拉折線法歐拉折線法 /* Eulers polygonal arc method*/ 幾何直觀是幫助我們尋幾何直觀是幫助我們尋找解決一個問題的思路找解決一個問題的思路的好辦法哦的好辦法哦定義定義在假設(shè)在假設(shè) yn = y(xn)，即第，即第 n 步計算是精確的前提下，考步計算是精確的前提下，考慮公式或方法本身帶來的誤差慮公式或方法本身帶來的誤差: Rn = y(xn+1) yn

5、+1 , 稱為稱為局部局部截斷誤差截斷誤差 /* local truncation error */。說明 顯然，這種近似有一定誤差，顯然，這種近似有一定誤差，而且步長越大，誤差越大，而且步長越大，誤差越大，如何估計這種誤差如何估計這種誤差y(xn+1) yn+1 ？1 Eulers Method截斷誤差截斷誤差: 實際上，實際上，y(xn) yn， yn 也有誤差，它對也有誤差，它對yn+1的誤差的誤差也有影響，見下圖。但這里不考慮此誤差的影響，僅考慮方也有影響，見下圖。但這里不考慮此誤差的影響，僅考慮方法或公式本身帶來的誤差，因此稱為方法誤差或截斷誤差。法或公式本身帶來的誤差，因此稱為方法

6、誤差或截斷誤差。局部截斷誤差的分析局部截斷誤差的分析：由于假設(shè)：由于假設(shè)yn = y(xn) ，即，即yn準確，因此準確，因此分析局部截斷誤差時將分析局部截斷誤差時將y(xn+1) 和和 yn+1都用點都用點xn上的信息來表上的信息來表示，工具：示，工具：Taylor展開。展開。 歐拉法的局部截斷誤差：歐拉法的局部截斷誤差：223111232() ( )( )( )( ) ( ,) ( )( )hnnnnnnnnnhnRy xyy xhy xy xO hyhf x yy xO hRn+1 的的主項主項/* leading term */1 Eulers Method定義定義若某算法的局部截斷誤

7、差為若某算法的局部截斷誤差為O(hp+1)，則稱該算法有，則稱該算法有p 階精度。階精度。 歐拉法具有歐拉法具有 1 階精度。階精度。在在xn點用一階向前差點用一階向前差商近似一階導(dǎo)數(shù)商近似一階導(dǎo)數(shù)1()()()nnny xy xyxh 在第二章討論牛頓插值公式時在第二章討論牛頓插值公式時介紹了介紹了差商差商的概念和性質(zhì)，的概念和性質(zhì)，各階差商可以近似各階導(dǎo)數(shù)，具有不同的精度，各階差商可以近似各階導(dǎo)數(shù)，具有不同的精度，且可以用函數(shù)值來表示。且可以用函數(shù)值來表示。上一章中數(shù)值微分的方法之一上一章中數(shù)值微分的方法之一就是用差商近似導(dǎo)數(shù)就是用差商近似導(dǎo)數(shù)111()()() ()()(,)nnnnnn

8、nnnny xy xhy xy xyy xyyh f xyEulers method1 Eulers Method1 Eulers Method 歐拉公式的改進歐拉公式的改進：隱式歐拉法或后退歐拉法隱式歐拉法或后退歐拉法 /* implicit Euler method or backward Euler method*/11()()()nnny xy xy xhxn+1點向后差商近似導(dǎo)數(shù)點向后差商近似導(dǎo)數(shù)111111()()() ()()(,)nnnnnnnnnny xy xhy xy xyy xyyh f xy隱式或后退歐拉公式隱式或后退歐拉公式由于未知數(shù)由于未知數(shù) yn+1 同時出現(xiàn)在等

9、式的兩邊，故稱為同時出現(xiàn)在等式的兩邊，故稱為隱式隱式 /* implicit */ 歐拉公式，而前者稱為歐拉公式，而前者稱為顯式顯式 /* explicit */ 歐拉公歐拉公式。隱式公式不能直接求解，一般需要用式。隱式公式不能直接求解，一般需要用Euler顯式公式顯式公式得到初值，然后用得到初值，然后用Euler隱式公式迭代求解。因此隱式公隱式公式迭代求解。因此隱式公式較顯式公式計算復(fù)雜，但穩(wěn)定性好（后面分析）。式較顯式公式計算復(fù)雜，但穩(wěn)定性好（后面分析）。01(1)( )111(,)(,)nnnnkknnnnyyh f xyyyh f xy(1)()1111111()(0 )1111(1)

10、11111()1(,)(,) 1, ()(,)kknnnnnnkknnnnknnnnnnknyyh fxyfxyhL yyhLyyhLyykyyh fxyy 在 迭 代 公 式 中 取 極 限 ， 有因 此的 極 限 就 是 隱 式 方 程 的 解收斂性收斂性1 Eulers Method 見上圖，見上圖， 顯然，這種近似也有一定誤差，顯然，這種近似也有一定誤差，如何估計這種誤差如何估計這種誤差y(xn+1) yn+1 ？方法同上，基于方法同上，基于Taylor展開估計局部截斷誤差。展開估計局部截斷誤差。但是注意，隱式公式中右邊含有但是注意，隱式公式中右邊含有f(xn+1 , yn +1 )

11、,由于由于yn +1不準確，所以不能直接用不準確，所以不能直接用y (xn+1)代替代替f(xn+1 , yn +1 ) 設(shè)已知曲線上一點設(shè)已知曲線上一點 Pn (xn , yn ),過該過該點作弦線，斜率為點作弦線，斜率為(xn+1 , yn +1 ) 點的點的方向場方向場f(x,y)方向方向,若步長若步長h充分小，充分小，可用弦線和垂線可用弦線和垂線x=xn+1的交點近似的交點近似曲線與垂線的交點。曲線與垂線的交點。幾何意義xnxn+1PnPn+1xyy(x)1 Eulers Method 隱式隱式歐拉法的局部截斷誤差：歐拉法的局部截斷誤差：11111111121111111321, ,2

12、 , 2nnnnynnnnnnnnnnnnynnnnnnnnfxyfxy xfxyy xyy xhfxy xyxyxhyxyxyhfxyy xy xhhyxh yxyxy xy x由微分中值定理，得在，之間；又而 2326nnnnhhhyxyxyx1 Eulers Method111111232311(), 231,23nnnynnnnnynnnnRy xyhfxy xyhhyxyxhhhfxRyxyx 從而即2211121,1,(1)1ynynynhfxhfxhfxxxx 111 Eulers Method 2322111231221 1,23 3,226 ()2nynynnnnynnnnn

13、hhRhfxhfxy xyxhhy xfxy xy xhRy xo h 隱式隱式歐拉法的局部截斷誤差：歐拉法的局部截斷誤差：111()nnnRy xy232()()hny xO h即隱式歐拉公式具有即隱式歐拉公式具有 1 階精度。階精度。1 Eulers Method1(,) 0, 1,.nnnnyyh f xyn比較尤拉顯式公式和隱式公式及其局部截斷誤差231112()()()hnnnnRy xyy xO h顯式公式111(,)nnnnyyh f xy隱式公式231112()()()hnnnnRy xyy xO h1 Eulers Method 若將這兩種方法進行算術(shù)平均，即可消除誤差若將這

14、兩種方法進行算術(shù)平均，即可消除誤差的主要部分的主要部分/*leading term*/而獲得更高的精度而獲得更高的精度,稱為梯形法稱為梯形法1 Eulers Method 梯形公式梯形公式 /* trapezoid formula */ 顯、隱式兩種算法的顯、隱式兩種算法的平均平均111 (,)(,)2nnnnnnhyyf xyf xy注：注：的確有局部截斷誤差的確有局部截斷誤差 ， 即梯形公式具有即梯形公式具有2 階精度，比歐拉方法有了進步。但階精度，比歐拉方法有了進步。但注意到該公式是注意到該公式是隱式隱式公式，計算時不得不用到迭代公式，計算時不得不用到迭代法，其迭代收斂性與歐拉公式相似。

15、法，其迭代收斂性與歐拉公式相似。3111()()nnnRy xyO h梯形法的迭代計算和收斂性梯形法的迭代計算和收斂性01(1)( )111(,)(,)(,)2nnnnkknnnnnnyyh f xyhyyf xyf xy(1)()1111111()(0 )1111(1)11111()1(,)(,)2222 1, ()2(,)(,)2kknnnnnnkknnnnknnnnnnnnknhyyfxyfxyhhLL yyyyhLhyykLhyyfxyfxyy 當 時 ，在 迭 代 公 式 中 取 極 限 ， 有因 此的 極 限 就 是 隱 式 方 程 的 解收斂性收斂性1 Eulers Method

16、梯形法的簡化計算梯形法的簡化計算 迭代計算量大，且難以預(yù)測迭代次數(shù)。為了控制計算量，通常只迭代迭代計算量大，且難以預(yù)測迭代次數(shù)。為了控制計算量，通常只迭代一次就轉(zhuǎn)入下一點的計算。用顯式公式作預(yù)測，梯形公式作校正，得到如下一次就轉(zhuǎn)入下一點的計算。用顯式公式作預(yù)測，梯形公式作校正，得到如下預(yù)測校正系統(tǒng)，也稱為改進尤拉法預(yù)測校正系統(tǒng)，也稱為改進尤拉法： 改進歐拉法改進歐拉法 /* modified Eulers method */Step 1: 先用先用顯式顯式歐拉公式作歐拉公式作預(yù)測預(yù)測，算出，算出),(1nnnnyxfhyy Step 2: 再將再將 代入代入隱式隱式梯形公式的右邊作梯形公式的右

17、邊作校正校正，得到，得到1 ny),(),(2111 nnnnnnyxfyxfhyy11( ,),( ,)2nnnnnnnnhyyf x yf xyh f x y1 Eulers Method注：注：此法亦稱為此法亦稱為預(yù)測預(yù)測-校正法校正法 /* predictor-corrector method */?？梢宰C明該算法具有可以證明該算法具有 2 階精度，同時可以看到它是個階精度，同時可以看到它是個單單步步遞推格式，比隱式公式的迭代求解過程遞推格式，比隱式公式的迭代求解過程簡單簡單。后面將。后面將看到，它的看到，它的穩(wěn)定性高穩(wěn)定性高于顯式歐拉法。于顯式歐拉法。1 Eulers Method+

18、1+11+1 (,)(,(,)21 (,)(,)2nnnnnnnnnnnnnnnhyyf xyf xh yh f xyyyhf xyyhf xy或其它形式其它形式1+1(,) (,)12pnnncnnpnpcyyh f xyyyh f xyyyy或幾何解釋幾何解釋xnxn+1ABPn+1=(A+B)/2尤拉法尤拉法后退尤拉法后退尤拉法梯形法梯形法1 Eulers Method 0)(,),(yaybaxyxfdxdy 00( ),xxy xyft y tdt令令x=x1,得得 1010(),xxy xyft y tdtAnother point of view對右端積分采用左矩形、右矩形、梯形

19、積分公式，即對右端積分采用左矩形、右矩形、梯形積分公式，即可得尤拉顯式、隱式、梯形公式可得尤拉顯式、隱式、梯形公式1 Eulers Method公式公式局部截斷誤差局部截斷誤差精精度度顯顯隱隱穩(wěn)穩(wěn)定定性性步數(shù)步數(shù)尤拉顯尤拉顯式公式式公式1 1階階顯顯差差單步單步尤拉隱尤拉隱式公式式公式1 1階階隱隱好好單步單步梯形公梯形公式式2 2階階隱隱差差單步單步中點法中點法2 2階階顯顯好好二步二步3(3)3nhyx3(3)12nhyx2(2)2nhyx2(2)2nhyxsummary例例 HW: p.201 #1-5證明中點法和梯形公式的精度為證明中點法和梯形公式的精度為2階階2 龍格龍格 - 庫塔法

20、庫塔法 /* Runge-Kutta Method */建立高精度的單步遞推格式建立高精度的單步遞推格式: :在改進尤拉法和尤拉在改進尤拉法和尤拉兩步法預(yù)測兩步法預(yù)測- -校正系統(tǒng)中校正系統(tǒng)中, ,預(yù)測公式都是單步法預(yù)測公式都是單步法, ,如果預(yù)如果預(yù)測誤差很小測誤差很小, ,則通過校正后得到的近似值誤差會更小則通過校正后得到的近似值誤差會更小, ,因因此需要研究高精度的單步法此需要研究高精度的單步法. . 1. Taylor級數(shù)法級數(shù)法 0)(,),(yaybaxyxfdxdyIVP:設(shè)其解為設(shè)其解為y=y(x)23126nnnnnhhy xy xhyxyxyx 由由Taylor展開，有展開

21、，有(1)2 Runge-Kutta Method(0)(0)(0)(1)(1)(1)(2)(2)(2)( )(1)( ,)jjjjyf x yfff dyffyffxy dxxyffyffxyffyffxy (2)2 Runge-Kutta Method要使公式具有要使公式具有p階精度，則在階精度，則在(1)式中截取前式中截取前p+1項，項，用用(2)式計算各階導(dǎo)數(shù)，即得下面式計算各階導(dǎo)數(shù)，即得下面Taylor公式：公式：2()12!ppnnnnnhhyyhyyyp 局部截斷誤差局部截斷誤差(3) 1(1)111, 1 !ppnnnnhy xyyxxp 2 Runge-Kutta Metho

22、d2.Taylor公式公式(3)表面上看形式簡單，但具體構(gòu)造時表面上看形式簡單，但具體構(gòu)造時往往很困難，因為按往往很困難，因為按(2)式求導(dǎo)，這一過程可能很復(fù)式求導(dǎo)，這一過程可能很復(fù)雜。因此通常不直接用雜。因此通常不直接用Taylor公式，而借鑒其思想公式，而借鑒其思想提出其它公式。提出其它公式。1. 由此看出，一種方法具有由此看出，一種方法具有p階精度階精度公式對不公式對不超過超過p次的多項式準確成立（局部截斷誤差為次的多項式準確成立（局部截斷誤差為0）。）。這一等價條件也可以用來判斷一種方法的精度。這一等價條件也可以用來判斷一種方法的精度。2 Runge-Kutta Method單步遞推法

23、的單步遞推法的基本思想基本思想是從是從 ( xn , yn ) 點出發(fā)，以點出發(fā)，以某一某一斜率斜率沿直線達到沿直線達到 ( xn+1 , yn+1 ) 點。歐拉法及其各種變形點。歐拉法及其各種變形所能達到的最高精度為所能達到的最高精度為2階階。 2. RungeKutta Method由微分中值定理，有由微分中值定理，有*1*1()( )(), 01 ()( )nnnnnnny xy xy xhf xh y xhkhy xy xhk k*稱為區(qū)間稱為區(qū)間xn, xn+1上的上的平均斜率平均斜率，只要知道平均斜，只要知道平均斜率，就可計算率，就可計算y(xn+1).因此只要對平均斜率提供一種因

24、此只要對平均斜率提供一種近似算法，則由近似算法，則由(4)式可導(dǎo)出一種相應(yīng)的求解公式。式可導(dǎo)出一種相應(yīng)的求解公式。(4) 2 Runge-Kutta Method *1*112*121. 0, ,Eulers formula2. 1, ,Eulers implicit formula3. Trapezoid formula2nnnnkf x y xkkf xy xkkkk例例121*12(,)4. (,) 2 modified Eulers method nnnnkf xykf xh yhkkkk2 Runge-Kutta Method 由此看出，改進的尤拉公式用由此看出，改進的尤拉公式用xn

25、與與xn+1兩個節(jié)點兩個節(jié)點的斜率的算術(shù)平均作為平均斜率，的斜率的算術(shù)平均作為平均斜率， xn+1點的斜率通點的斜率通過已知信息過已知信息yn來預(yù)測。來預(yù)測。 考察改進的歐拉法，可以將其改寫為：考察改進的歐拉法，可以將其改寫為：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii 斜率斜率一定取一定取K1， K2 的的平均值平均值嗎？嗎？步長一定是步長一定是h 嗎？即嗎？即第二個節(jié)點一定是第二個節(jié)點一定是xn+1嗎？嗎？2 Runge-Kutta Method2 Runge-Kutta Method首先希望能確定系數(shù)首先希望能確定系數(shù) 1、 2、p，使得到的算法格式有，

26、使得到的算法格式有2階階精度，即在精度，即在 的前提假設(shè)下，使得的前提假設(shè)下，使得 )(iixyy )()(311hOyxyRiii Step 1: 將將 K2 在在 ( xi , yi ) 點作點作 Taylor 展開展開)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 將改進歐拉法推廣為：將改進歐拉法推廣為：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step

27、2: 將將 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii 2階階RungeKutta Method2 Runge-Kutta MethodStep 3: 將將 yi+1 與與 y( xi+1 ) 在在 xi 點的點的泰勒泰勒展開作比較展開作比較)()()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhxyhxyxyiiii 要求要求 ，則必須有：，則必須有：)()(311hOyxyRiii21,1221 p 這里有這里有 個未知個未知數(shù)，數(shù)， 個方程

28、。個方程。32存在存在無窮多個解無窮多個解。所有滿足上式的格式統(tǒng)稱為。所有滿足上式的格式統(tǒng)稱為2階龍格階龍格 - 庫庫塔格式塔格式。21, 121 p注意到，注意到， 就是改進的歐拉法；就是改進的歐拉法；p=1/2, 1=0, 2=1, 變形尤拉公式變形尤拉公式。Q: 為獲得更高的精度，應(yīng)該如何進一步推廣？為獲得更高的精度，應(yīng)該如何進一步推廣？改進的改進的Euler 公式推廣為二階公式推廣為二階Runge-Kutta公式公式帶來這樣的啟示：帶來這樣的啟示：若在若在xn, xn+1上多預(yù)測幾個點的斜率值，然后將上多預(yù)測幾個點的斜率值，然后將它們的算術(shù)平均作為平均斜率，則有可能構(gòu)造出它們的算術(shù)平均

29、作為平均斜率，則有可能構(gòu)造出具有更高精度的計算公式。具有更高精度的計算公式。-Runge-Kutta方法的基本思想。方法的基本思想。注：二階注：二階Runge-Kutta公式用多算一次函數(shù)值公式用多算一次函數(shù)值f 的的辦法避開了二階辦法避開了二階Taylor級數(shù)法所要計算的級數(shù)法所要計算的f 的導(dǎo)數(shù)。的導(dǎo)數(shù)。在這種意義上，可以說在這種意義上，可以說Runge-Kutta方法實質(zhì)上是方法實質(zhì)上是Taylor級數(shù)法的變形。級數(shù)法的變形。2 Runge-Kutta Method其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j

30、= 1, , i 1 ) 均為待定系數(shù)，確均為待定系數(shù)，確定這些系數(shù)的步驟與前面相似。定這些系數(shù)的步驟與前面相似。 2 Runge-Kutta Method).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 高階高階RungeKutta Method Gill公式：公式：4階經(jīng)典龍格階經(jīng)典龍格-庫塔公式的一種改進庫塔公式的一種改進112346121222 1231222222423222222( ,)(,)(,1)(,1)hiiiihhiihiiiiy

31、yKKKKKf xyKf xyKKf xyhKhKKf xh yhKhK2 Runge-Kutta Method 最常用為四級最常用為四級4階階經(jīng)典龍格經(jīng)典龍格-庫塔法庫塔法 /* Classical Runge-Kutta Method */ ：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 2 Runge-Kutta Method注：注： 龍格龍格-庫塔法庫塔法的主要運算在于計算的主要運算在于計算 Ki 的值，即計算的值，即計算 f 的的值。值。Butcher 于于1965年給出了計算量

32、與可達到的最高精年給出了計算量與可達到的最高精度階數(shù)的關(guān)系：度階數(shù)的關(guān)系：753可達到的最高精度可達到的最高精度642每步須算每步須算Ki 的個數(shù)的個數(shù))(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n 由于龍格由于龍格-庫塔法的導(dǎo)出基于泰勒展開，故精度主要受庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對于光滑性不太好的解，最好解函數(shù)的光滑性影響。對于光滑性不太好的解，最好采用采用低階算法低階算法而將步長而將步長h 取小取小。2 Runge-Kutta Method 變步長的變步長的RungeKutta Method1(1)1ppnnRchyxQ: 由局

33、部截斷誤差可以看出，步長由局部截斷誤差可以看出，步長 h 越小，局部截斷越小，局部截斷誤差越??；但步長減小，在一定求解范圍（區(qū)間）內(nèi)誤差越??；但步長減小，在一定求解范圍（區(qū)間）內(nèi)要完成的步數(shù)就增加了，步數(shù)增加會引起計算量增大，要完成的步數(shù)就增加了，步數(shù)增加會引起計算量增大，導(dǎo)致舍入誤差積累。因此要選取適當?shù)牟介L。導(dǎo)致舍入誤差積累。因此要選取適當?shù)牟介L。選擇步長時要考慮兩個問題：選擇步長時要考慮兩個問題： 1.如何衡量和檢驗計算結(jié)果的精度？如何衡量和檢驗計算結(jié)果的精度？ 2.如何根據(jù)所獲得的精度處理步長？如何根據(jù)所獲得的精度處理步長？3 單步法的單步法的收斂性與穩(wěn)定性收斂性與穩(wěn)定性 /* Con

34、vergency and Stability */ 前面介紹了兩大類微分方程數(shù)值解法：一類是前面介紹了兩大類微分方程數(shù)值解法：一類是用差商近似導(dǎo)數(shù)得到的尤拉系列公式，另一類是用差商近似導(dǎo)數(shù)得到的尤拉系列公式，另一類是基于平均斜率概念的基于平均斜率概念的RungeKutta公式。公式?；舅蓟舅枷攵际峭ㄟ^某種離散化手續(xù)，將微分方程轉(zhuǎn)化為想都是通過某種離散化手續(xù)，將微分方程轉(zhuǎn)化為差分方程（代數(shù)方程）來求解差分方程（代數(shù)方程）來求解。 Q1. 這種轉(zhuǎn)化是否合理？要看差分問題的解這種轉(zhuǎn)化是否合理？要看差分問題的解yn當當h0時是否收斂到微分方程的解時是否收斂到微分方程的解y(xn),即是否成立即是否

35、成立 yn y(xn)， h0. -收斂性問題收斂性問題 Q2. 實際計算時，由于舍入誤差的影響，差分方實際計算時，由于舍入誤差的影響，差分方程的解本身也有誤差，這種誤差在計算過程中會程的解本身也有誤差，這種誤差在計算過程中會不會擴大不會擴大 -穩(wěn)定性問題穩(wěn)定性問題3 Convergency and Stability 收斂性收斂性 /* Convergency */定義定義 若某算法對于任意固定的若某算法對于任意固定的 x = xi = x0 + i h，當，當 h0 ( 同時同時 i ) 時有時有 yi y( xi )，則稱該算法是，則稱該算法是收斂收斂的。的。 例：例：就初值問題就初值問

36、題 考察歐拉顯式格式的收斂性?？疾鞖W拉顯式格式的收斂性。 0)0(yyyy 解：解：該問題的精確解為該問題的精確解為 xeyxy 0)( 歐拉公式為歐拉公式為iiiiyhyhyy)1 (1 0)1 (yhyii 對任意固定的對任意固定的 x = xi = i h ，有，有iixhhxihyhyy )1()1(/10/0 ehhh /10)1(lim)(0ixxyeyi 3 Convergency and Stability顯式單步法的收斂性顯式單步法的收斂性(1)3 Convergency and Stability而整體截斷誤差為而整體截斷誤差為1111111nnnnnnney xyy xy

37、yy證明證明: 設(shè)設(shè) 表示當表示當yn =y(xn)時，時， 由公式由公式(1)求得的求得的結(jié)果，即結(jié)果，即1,nnnnyy xhxy xh1ny則局部截斷誤差為則局部截斷誤差為111pnny xych(2)(2)-(1),得得 11,11nnnnnnnnnnnnnnnyyy xyhx y xhx y hy xyhL y xyhLy xyhLe3 Convergency and Stability從而從而 11111211112010100111111111111111pppnnnppnnnnpnnnpnnpehLechhLhLechchhLehLchchehLehLhLhLchhLhLech

38、hLhLhLechL3 Convergency and Stability000 ()1 ,1,1011 0 1 ()nnnhLTLnxnxTLTLpnTLpnpnnxxnhTbahLeexRxexxeeee echeLceheLy xyo h 當時，有當時， 而 ，3 Convergency and Stability Euler法的收斂性法的收斂性 ： (x,y)=f(x,y),故當故當f(x,y)滿足滿足Lipschitz條件時，尤拉條件時，尤拉法收斂；法收斂；3 Convergency and Stability3 Convergency and Stability 穩(wěn)定性穩(wěn)定性 /*

39、 Stability */例：例：考察初值問題考察初值問題 在區(qū)間在區(qū)間0, 0.5上的解。上的解。分別用歐拉顯、隱式格式和改進的歐拉格式計算數(shù)值解。分別用歐拉顯、隱式格式和改進的歐拉格式計算數(shù)值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精確解精確解改進歐拉法改進歐拉法 歐拉隱式歐拉隱式歐拉顯式歐拉顯式 節(jié)點節(jié)點 xixey30 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6.2500 10 21.5625 10 23.9063 10 39.7656 10 41.00002

40、.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 7What is wrong ?! An Engineer complains: Math theorems are so unstable that a small perturbation on the conditions will cause a crash on the conclusions!3 Convergency and Stability3 Convergency and S

41、tability定義定義若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計算中都算中都逐步衰減逐步衰減，則稱該算法是，則稱該算法是絕對穩(wěn)定的絕對穩(wěn)定的 /*absolutely stable */。一般分析某算法的穩(wěn)定性時，為簡單起見，只考慮一般分析某算法的穩(wěn)定性時，為簡單起見，只考慮模型方程模型方程或或試驗方程試驗方程 /* test equation */yy 常數(shù)，可以常數(shù)，可以是復(fù)數(shù)是復(fù)數(shù)當步長取為當步長取為 h 時，將某算法應(yīng)用于上式，并假設(shè)只在初值時，將某算法應(yīng)用于上式，并假設(shè)只在初值產(chǎn)生誤差產(chǎn)生誤差 ，則若此誤差以后逐步衰減，就稱該，則若

42、此誤差以后逐步衰減，就稱該算法相對于算法相對于 絕對穩(wěn)定絕對穩(wěn)定， 的全體構(gòu)成的全體構(gòu)成絕對穩(wěn)定區(qū)域絕對穩(wěn)定區(qū)域。我們稱我們稱算法算法A 比算法比算法B 穩(wěn)定穩(wěn)定，就是指，就是指 A 的絕對穩(wěn)定區(qū)域比的絕對穩(wěn)定區(qū)域比 B 的的大大。000yy h h h3 Convergency and Stability3 Convergency and Stability3 Convergency and Stability3 Convergency and Stability例：例：隱式龍格隱式龍格-庫塔法庫塔法 ),., 1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmii

43、 而而顯式顯式 1 4 階方法的絕對穩(wěn)定階方法的絕對穩(wěn)定區(qū)域為區(qū)域為 )2,2(1111KhyhxfKhKyyiiii其中其中2階方法階方法 的絕對穩(wěn)定區(qū)域為的絕對穩(wěn)定區(qū)域為0ReImgk=1k=2k=3k=4-1-2-3-123ReImg無條件穩(wěn)定無條件穩(wěn)定注：注：一般來說，隱式歐拉法的絕對穩(wěn)定性比同階的顯式法的好。一般來說，隱式歐拉法的絕對穩(wěn)定性比同階的顯式法的好。小結(jié)小結(jié) 尤拉兩步公式與尤拉單步公式相比，使用兩尤拉兩步公式與尤拉單步公式相比，使用兩個節(jié)點上的已知信息將精度提高一階?？梢栽O(shè)想，個節(jié)點上的已知信息將精度提高一階?？梢栽O(shè)想，計算計算y(xn+1)時，充分利用前面已經(jīng)求出的節(jié)點上

44、時，充分利用前面已經(jīng)求出的節(jié)點上的的 y 及及 y 值的值的線性組合線性組合來近似來近似y(xn+1)，精度會大，精度會大大提高。大提高。-線性多步法的基本思想線性多步法的基本思想。 構(gòu)造線性多步法有多種途徑，這里介紹兩種：構(gòu)造線性多步法有多種途徑，這里介紹兩種： 基于數(shù)值積分的構(gòu)造方法基于數(shù)值積分的構(gòu)造方法； 基于基于TaylorTaylor展開的構(gòu)造方法展開的構(gòu)造方法。4 線性多步法線性多步法 /* Multistep Method */).(.110111101rnrnnnrnrnnnffffhyyyy 線性多步法的通式可寫為：線性多步法的通式可寫為：),(jjjyxff 當當 1 0

45、時，為時，為隱式公隱式公式式; 1=0 則為則為顯式公式顯式公式。 基于數(shù)值積分的構(gòu)造法基于數(shù)值積分的構(gòu)造法將將 在在 上積分，得到上積分，得到),(yxfy 1,nnx x11()()( , ( )nnxnnxy xy xf x y x dx只要只要近似地算出右邊的積分近似地算出右邊的積分 ，則可通，則可通過過 近似近似y(xn+1) 。而。而選用不同近似式選用不同近似式 I，可得到不，可得到不同的計算公式同的計算公式。例如利用左矩形積分公式得到尤拉公式；。例如利用左矩形積分公式得到尤拉公式；梯形積分公式得到梯形公式。梯形積分公式得到梯形公式。1( , ( )nnxxIf x y x dx1

46、nnyyI一般地，利用插值原理所建立的一系列數(shù)值積分一般地，利用插值原理所建立的一系列數(shù)值積分方法也可以導(dǎo)出解微分方程的一系列計算公式。方法也可以導(dǎo)出解微分方程的一系列計算公式。運用插值方法的關(guān)鍵在于選取合適的插值節(jié)點。運用插值方法的關(guān)鍵在于選取合適的插值節(jié)點。假設(shè)已構(gòu)造出假設(shè)已構(gòu)造出f(x,y(x)的插值多項式的插值多項式Pr(x),則則1111( , ( )( )( )nnnnnnxxrxxxnnrxf x y x dxP x dxyyP x dx4 Multistep Method 亞當姆斯顯式公式亞當姆斯顯式公式 /* Adams explicit formulae */利用利用r+1

47、 個節(jié)點上的被積函數(shù)值個節(jié)點上的被積函數(shù)值構(gòu)造構(gòu)造 r 階牛頓階牛頓后插后插多項式多項式, 有有 11,nnnnn rn rxfxfxf11100( , ( )()()nnnx xthxrnrnxf x y x dxN xth hdtR xth hdtNewton插值余項插值余項1110001000()11rjjnnrnnn jjrrjjjnn jnjn jjjtyyhN xth dtyhfdtjtyhdtfyhfj0()11, rjjrnnjjjntNxthfjxxtjh ，其中0表示 階向前差分，/*Adams 顯式公式顯式公式 */外推技術(shù)外推技術(shù) /* extrapolation */

48、10 is independent of n and r, 1jjjtdtjwhere table 5-6, p.181j j0111/225/1233/8實際計算時，將差分展開實際計算時，將差分展開 000000111jijnjn iijrrrrriijjnjjn ijn irin ijjiij iijffijjffffii ijrj=i10rnnrin iiyyhf局部截斷誤差為：局部截斷誤差為： 1110(2)12(2)102(2)()()()1 !()rnnrnrnrrrrnrrrnRy xyhR xth dtyht dtB hyrB hyx注：注：Br 與與yn+1 計算公式中計算公

49、式中 fn , , fn k 各項的各項的系數(shù)系數(shù) 均可查表均可查表得到得到 。 見下表。見下表。ri10123r2123122324552112162459125243724912583720251fnfn 1fn 2fn 3Br例：例：1110,1,(3)2nnnnnnnryyhfhryyff35()12nRh yx常用的是常用的是 r = 3 的的4階亞當姆斯顯式公式階亞當姆斯顯式公式)9375955(243211 iiiiiiffffhyy5(5)251()720nRh yx21()2nRh y xAdams顯式公式用顯式公式用 作為插作為插值節(jié)點，在求積區(qū)間值節(jié)點，在求積區(qū)間xn,

50、xn+1上用插值函數(shù)上用插值函數(shù)Nr(x)近近似似f(x,y(x)，而，而xn+1不在插值節(jié)點內(nèi)，因此是一個不在插值節(jié)點內(nèi)，因此是一個外推外推的過程。雖然公式是顯式的，便于計算，但的過程。雖然公式是顯式的，便于計算，但效果并不理想，比如穩(wěn)定性較差等。效果并不理想，比如穩(wěn)定性較差等。因此改用通過因此改用通過r+1個節(jié)點個節(jié)點的插值多項式的插值多項式Nr(x)近似近似f(x,y(x)，由于，由于xn+1是其中是其中一個插值點，因此是一個插值點，因此是內(nèi)插多項式內(nèi)插多項式，但導(dǎo)出的公式，但導(dǎo)出的公式是是隱式的隱式的。 11,nnnnn rn rxfxfxf 1111,nnnnn rn rxfxfxf

51、 4 Multistep Method 亞當姆斯隱式公式亞當姆斯隱式公式 /* Adams implicit formulae */利用利用r+1 個節(jié)點上的被積函數(shù)值個節(jié)點上的被積函數(shù)值 fn+1 , fn , , fn r+1 構(gòu)造構(gòu)造r 階階牛頓牛頓前插前插多項式。與顯式多項式完全類似地可得到一系列多項式。與顯式多項式完全類似地可得到一系列隱式公式。隱式公式。0*1110*110, 1, 1rjjnnjnjjjrrinnrinirijij ityyhfdtjjyyhfi 差分展開，得j0 11 -1/22 -1/123 -1/24*j局部截斷誤差為：局部截斷誤差為：2(2)112(2)(

52、)()()rrrnnrnrrrnRy xyB hyB hyx10123k21 21125249211282419121 245 241121 241 72019 fi+1fifi 1fi 2Br小于小于Br注：注： 與與yn+1 計算公式中計算公式中 fn+1 , , fn+1 k 各項的各項的系數(shù)系數(shù) 均均可查表得到可查表得到 。 見下表。見下表。*rirB常用的是常用的是 k = 3 的的4階亞當姆斯隱式公式階亞當姆斯隱式公式)5199(242111 iiiiiiffffhyy較同階顯式較同階顯式穩(wěn)定穩(wěn)定例：例：11110,1,()2nnnnnnnryyhfhryyff31()12nRh

53、yx 21()2nRh y x 5(5)19()720nRh yx 4 Multistep Method 亞當姆斯預(yù)測亞當姆斯預(yù)測-校正系統(tǒng)校正系統(tǒng) /* Adams predictor-corrector system */Step 1: 用用Runge-Kutta 法法計算前計算前 k 個初值；個初值；Step 2: 用用Adams 顯式顯式計算計算預(yù)測預(yù)測值；值；Step 3: 用同階用同階Adams 隱式隱式計算計算校正校正值。值。注意：注意：三步所用公三步所用公式的精度必須相同。式的精度必須相同。通常用通常用經(jīng)典經(jīng)典Runge-Kutta 法法配合配合4階階Adams 公式公式。4階

54、階Adams隱式公式的截斷誤差為隱式公式的截斷誤差為5(5)1119()()720nnny xyh yx 4階階Adams顯式公式的截斷誤差為顯式公式的截斷誤差為5(5)11251()()720nnny xyh yx1111()19()251nnnny xyy xy Predicted value pn+1Corrected value cn+11111251()()270nnnny xyyy 111119()()270nnnny xyyyModified value mn+1預(yù)測值與校正預(yù)測值與校正值的事后誤差值的事后誤差估計估計1111251()()270nnnny xyyy111119(

55、)()270nnnny xyyy11231111121111(5559379)24251()270(9195)2419()270nnnnnnnnnnnnnnnnnnnnhpyffffmppchcyffffycpc預(yù)測：改進：校正：改進：Adams顯隱預(yù)測顯隱預(yù)測-校正校正-改進系統(tǒng)改進系統(tǒng)4 Multistep Method Adams 4th-Order predictor-corrector AlgorithmTo approximate the the solution of the initial-value problemAt (N+1) equally spaced numbers

56、 in the interval a, b.Input: endpoints a, b; integer N; initial value y0 .Output: approximation y at the (N+1) values of x.Step 1 Set h = (b a) / N ; x0 = a; y0 = y0; Output ( x0, y0 );Step 2 For i = 1, 2, 3 Compute yi using classical Runge-Kutta method; Output ( xi , yi );Step 3 For i = 4, , N do s

57、teps 4-10Step 5 ; /* predict */ Step 6 ; /* modify */Step 7 ; /* correct */Step 8 ; /* modify the final value */Step 9 Output ( xi+1 , yi+1 ); Step 10 For j = 0, 1, 2, 3 Set xi = xi+1 ; yi = yi+1 ; /* Prepare for next iteration */Step 11 STOP. 0)(,),(yaybxayxfy 24/ )9375955(3211 iiiiiiffffhyp270/ )(

58、25111iiiipcpm 24/ )519),(9(21111 iiiiiiifffmxfhyc270/ )(191111 iiiipccy應(yīng)為應(yīng)為( ci+1 pi+1 ), 但因但因ci+1 尚未尚未算出，只好用算出，只好用( ci pi )取代之。取代之。4 Multistep Method 基于泰勒展開的構(gòu)造法基于泰勒展開的構(gòu)造法).(.110111101rnrnnnrnrnnnffffhyyyy 將通式中的右端各項將通式中的右端各項 yn 1, , yn k ; fn+1, fn 1, , fn k (yn+1, , yn k)分別在分別在 xn 點作點作泰勒泰勒展開展開，與精確解

59、，與精確解 y(xn+1) 在在 xn 點的泰勒展開作點的泰勒展開作比較比較。通過令通過令同類項系數(shù)相等同類項系數(shù)相等，得到足以確定待定系數(shù)，得到足以確定待定系數(shù) 0, , r ; 1, 0, , r 的等式，則可構(gòu)造出線性的等式，則可構(gòu)造出線性多步法的公式。多步法的公式。),(jjjyxff 當當 1 0 時，為時，為隱式公隱式公式式; 1=0 則為則為顯式公式顯式公式。101011( )(1)01( )(1)110 ,Taylor!1 !1 !rrrrnkn kkn kkn kkn kkkkkn kn kn kn kjppjpn knnjjppjpn knnjrnknkyyhfyhyyy

60、xyyxkhkhyyyjpkhkhyyyjpyy設(shè)由展開，有1( )11111(1)111( )(1)10!11 !1 !jprrjjjkknjkkprrpppkknkkjppjpnnnjhkjkyjhkpkyphhy xyyjp111(),1pnnpo hyy xp要使公式具有 階精度，即局部截斷誤差為只要令與的前項相同即可，即011111 1 2,rkkrrjjkkkkkjkjp， ，11(1)1111111 !prrpppnnkknkkhy xykpkyp局部截斷誤差為局部截斷誤差為例：例：設(shè)設(shè))(3322110221101 nnnnnnnnyyyyhyyyy 確定式中待定系數(shù)確定式中待