常微分方程數(shù)值解法歐拉法

會(huì)計(jì)學(xué)1常微分方程數(shù)值解法歐拉法解析解法：（常微分方程理論）只能求解極少一類常微分方程；實(shí)際中給定的問(wèn)題不一定是解析表達(dá)式，而是函數(shù)表，無(wú)法用解析解法。如何求解計(jì)算解函數(shù)y(x)在一系列節(jié)點(diǎn)a=x0<x1<…<xn=b

處的近似值節(jié)點(diǎn)間距為步長(zhǎng)，通常采用等距節(jié)點(diǎn)，即取hi=h

(常數(shù))。數(shù)值解法:求解所有的常微分方程步進(jìn)式：根據(jù)已知的或已求出的節(jié)點(diǎn)上的函數(shù)值計(jì)算當(dāng)前節(jié)點(diǎn)上的函數(shù)值，一步一步向前推進(jìn)。因此只需建立由已知的或已求出的節(jié)點(diǎn)上的函數(shù)值求當(dāng)前節(jié)點(diǎn)函數(shù)值的遞推公式即可。第1頁(yè)/共32頁(yè)第2頁(yè)/共32頁(yè)第3頁(yè)/共32頁(yè)--------Euler’sMethod§1歐拉方法

/*Euler’sMethod*/§1Euler’sMethodTaylor展開(kāi)法第4頁(yè)/共32頁(yè)幾何意義亦稱為歐拉折線法

/*Euler’spolygonalarcmethod*/

幾何直觀是幫助我們尋找解決一個(gè)問(wèn)題的思路的好辦法哦第5頁(yè)/共32頁(yè)定義在假設(shè)yn=y(xn)，即第

n

步計(jì)算是精確的前提下，考慮公式或方法本身帶來(lái)的誤差:Rn=y(xn+1)

yn+1,稱為局部截?cái)嗾`差/*localtruncationerror*/。說(shuō)明

顯然，這種近似有一定誤差，而且步長(zhǎng)越大，誤差越大，如何估計(jì)這種誤差y(xn+1)

yn+1

？§1Euler’sMethod第6頁(yè)/共32頁(yè)截?cái)嗾`差:實(shí)際上，y(xn)

yn，yn也有誤差，它對(duì)yn+1的誤差也有影響，見(jiàn)下圖。但這里不考慮此誤差的影響，僅考慮方法或公式本身帶來(lái)的誤差，因此稱為方法誤差或截?cái)嗾`差。局部截?cái)嗾`差的分析：由于假設(shè)yn=y(xn)

，即yn準(zhǔn)確，因此分析局部截?cái)嗾`差時(shí)將y(xn+1)和

yn+1都用點(diǎn)xn上的信息來(lái)表示，工具：Taylor展開(kāi)。

歐拉法的局部截?cái)嗾`差：Rn+1

的主項(xiàng)/*leadingterm*/§1Euler’sMethod第7頁(yè)/共32頁(yè)定義若某算法的局部截?cái)嗾`差為O(hp+1)，則稱該算法有p

階精度。歐拉法具有1階精度。在xn點(diǎn)用一階向前差商近似一階導(dǎo)數(shù)

在第2章討論牛頓插值公式時(shí)介紹了差商的概念和性質(zhì)，各階差商可以近似各階導(dǎo)數(shù)，具有不同的精度，且可以用函數(shù)值來(lái)表示。上一章中數(shù)值微分的方法之一就是用差商近似導(dǎo)數(shù)Euler’smethod§1Euler’sMethod第8頁(yè)/共32頁(yè)§1Euler’sMethod

歐拉公式的改進(jìn)：隱式歐拉法或后退歐拉法

/*implicitEulermethodorbackwardEulermethod*/xn+1點(diǎn)向后差商近似導(dǎo)數(shù)隱式或后退歐拉公式第9頁(yè)/共32頁(yè)由于未知數(shù)yn+1

同時(shí)出現(xiàn)在等式的兩邊，故稱為隱式/*implicit*/

歐拉公式，而前者稱為顯式/*explicit*/歐拉公式。隱式公式不能直接求解，一般需要用Euler顯式公式得到初值，然后用Euler隱式公式迭代求解。因此隱式公式較顯式公式計(jì)算復(fù)雜，但穩(wěn)定性好收斂性§1Euler’sMethod第10頁(yè)/共32頁(yè)

見(jiàn)上圖，顯然，這種近似也有一定誤差，如何估計(jì)這種誤差y(xn+1)

yn+1

？方法同上，基于Taylor展開(kāi)估計(jì)局部截?cái)嗾`差。但是注意，隱式公式中右邊含有f(xn+1

,yn+1),由于yn+1不準(zhǔn)確，所以不能直接用y'(xn+1)代替f(xn+1

,yn+1)設(shè)已知曲線上一點(diǎn)Pn(xn,yn),過(guò)該點(diǎn)作弦線，斜率為(xn+1

,yn+1)點(diǎn)的方向場(chǎng)f(x,y)方向,若步長(zhǎng)h充分小，可用弦線和垂線x=xn+1的交點(diǎn)近似曲線與垂線的交點(diǎn)。幾何意義xnxn+1PnPn+1xyy(x)§1Euler’sMethod第11頁(yè)/共32頁(yè)

隱式歐拉法的局部截?cái)嗾`差：§1Euler’sMethod第12頁(yè)/共32頁(yè)§1Euler’sMethod第13頁(yè)/共32頁(yè)

隱式歐拉法的局部截?cái)嗾`差：即隱式歐拉公式具有1階精度?！?Euler’sMethod第14頁(yè)/共32頁(yè)比較尤拉顯式公式和隱式公式及其局部截?cái)嗾`差顯式公式隱式公式§1Euler’sMethod第15頁(yè)/共32頁(yè)

若將這兩種方法進(jìn)行算術(shù)平均，即可消除誤差的主要部分/*leadingterm*/而獲得更高的精度,稱為梯形法§1Euler’sMethod梯形公式/*trapezoidformula*/—顯、隱式兩種算法的平均注：的確有局部截?cái)嗾`差，即梯形公式具有2

階精度，比歐拉方法有了進(jìn)步。但注意到該公式是隱式公式，計(jì)算時(shí)不得不用到迭代法，其迭代收斂性與歐拉公式相似。第16頁(yè)/共32頁(yè)梯形法的迭代計(jì)算和收斂性收斂性§1Euler’sMethod第17頁(yè)/共32頁(yè)梯形法的簡(jiǎn)化計(jì)算

迭代計(jì)算量大，且難以預(yù)測(cè)迭代次數(shù)。為了控制計(jì)算量，通常只迭代一次就轉(zhuǎn)入下一點(diǎn)的計(jì)算。用顯式公式作預(yù)測(cè)，梯形公式作校正，得到如下預(yù)測(cè)校正系統(tǒng)，也稱為改進(jìn)尤拉法：改進(jìn)歐拉法/*modifiedEuler’smethod*/Step1:

先用顯式歐拉公式作預(yù)測(cè)，算出),(1nnnnyxfhyy+=+Step2:再將代入隱式梯形公式的右邊作校正，得到1+ny)],(),([2111+++++=nnnnnnyxfyxfhyy§1Euler’sMethod第18頁(yè)/共32頁(yè)注：此法亦稱為預(yù)測(cè)-校正法/*predictor-correctormethod*/。可以證明該算法具有2階精度，同時(shí)可以看到它是個(gè)單步遞推格式，比隱式公式的迭代求解過(guò)程簡(jiǎn)單。后面將看到，它的穩(wěn)定性高于顯式歐拉法?！?Euler’sMethod其它形式第19頁(yè)/共32頁(yè)幾何解釋xnxn+1ABPn+1=(A+B)/2尤拉法后退尤拉法梯形法§1Euler’sMethod第20頁(yè)/共32頁(yè)令x=x1,得Anotherpointofview對(duì)右端積分采用左矩形、右矩形、梯形積分公式，即可得尤拉顯式、隱式、梯形公式§1Euler’sMethod第21頁(yè)/共32頁(yè)中點(diǎn)歐拉公式/*midpointformula*/中心差商近似導(dǎo)數(shù)x0x2x1假設(shè)，則可以導(dǎo)出即中點(diǎn)公式也具有2階精度，且是顯式的。需要2個(gè)初值y0和y1來(lái)啟動(dòng)遞推過(guò)程，這樣的算法稱為雙步法/*double-stepmethod*/，而前面的三種算法都是單步法/*single-stepmethod*/?！?Euler’sMethod第22頁(yè)/共32頁(yè)幾何解釋xnxn+1尤拉法后退尤拉法中點(diǎn)法xn-1令x=x2,得Anotherpointofview對(duì)右端積分采用中矩形公式即得中點(diǎn)公式§1Euler’sMethod第23頁(yè)/共32頁(yè)公式局部截?cái)嗾`差精度顯隱穩(wěn)定性步數(shù)尤拉顯式公式1階顯差單步尤拉隱式公式1階隱好單步梯形公式2階隱差單步中點(diǎn)法2階顯好二步summary第24頁(yè)/共32頁(yè)第25頁(yè)/共32頁(yè)第26頁(yè)/共32頁(yè)

算例

分別用顯式Euler方法，梯形方法和預(yù)估－校正Euler方法解初值問(wèn)題解：取

h=0.1，梯形方法為：續(xù)第27頁(yè)/共32頁(yè)

算例

分別用顯式Euler方法，梯形方法和預(yù)估－校正Euler方法解初值問(wèn)題解：取

h=0.1，梯形方法為：預(yù)估－校正Euler方法：續(xù)第28頁(yè)/共32頁(yè)

Euler方法

梯形方法

預(yù)估－校正方法0.01.0000000.01.0000000.01.0000000.00.11.0000004.8×10－31.0047627.5×10－51.0050001.6×10－40.21.0100008.7×10－31.0185941.4×10－41.0190252.9×10－40.31.0290001.2×10－21.0406331.9×10－41.0412184.0×10－40.41.0561001.4×10－21.0700962.2×10－41.0708004.8×10－40.51.0904901.6×10－21.1062782.5×10－41.1070765.5×10－40.61.1314411.7×10－21.1485372.7×10－41.1494045.9×10－40.71.1782971.8×10－21.1962952.9×10－41.1972106.2×10－40.81.2304671.9×10－21.2490193.0×10－41.2499756.5×10－40.91.2874201.9×10－21.3062643.1×10－41.3072286.6×10－41.01.3486781.9×10－21.3675733.1×