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

1、數(shù)值分析Numerical Analysis第八章常微分方程數(shù)值解法鄭州大學(xué)研究生課程 （2010-2011學(xué)年第一學(xué)期） ISCM 2007,Beijing China1第八章 常微分方程數(shù)值解法 8.1 引言8.2 歐拉(Euler)法8.3 改進(jìn)歐拉(Euler)方法8.4 單步法的穩(wěn)定性ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.1 引言問題提出 倒葫蘆形狀容器壁上的刻度問題.對于圓柱形狀容器壁上的容積刻度,可以利用圓柱體體積公式其中直徑D為常數(shù).由于體積V與相對于容器底部的任意高度H

2、的函數(shù)關(guān)系明確,因此在容器上可以方便地標(biāo)出容器刻度,而對于幾何形狀不是規(guī)則的容器,比如倒葫蘆形狀容器壁上如何標(biāo)出刻度呢?ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.1 引言下表是經(jīng)過測量得到部分容器高度與直徑的關(guān)系.H 0 0.2 0.4 0.6 0.8 1.0D 0 0.11 0.26 0.56 1.04 1.17根據(jù)上表的數(shù)據(jù),可以擬合出倒葫蘆形狀容器的圖,建立如圖所示的坐標(biāo)軸后,問題即為如何根據(jù)任意高度x標(biāo)出容器體積V的刻度,由微元思想分析可知ISCM 2007,Beijing Chin

3、a/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.1 引言其中x表示高度,直徑D是高度x的函數(shù),記為D(x),因此得到如下微分方程初值問題只要求解上述方程,就可求出體積V與高度x之間的函數(shù)關(guān)系,從而可標(biāo)出容器壁上容積的刻度,但問題是函數(shù)D(x)無解析表達(dá)式,我們無法求出其解析解.ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.1 引言 包含自變量、未知函數(shù)及未知函數(shù)的導(dǎo)數(shù)或微分的方程稱為微分方程。在微分方程中, 自變量的個數(shù)只有一個, 稱為常微分

4、方程。自變量的個數(shù)為兩個或兩個以上的微分方程叫偏微分方程。微分方程中出現(xiàn)的未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)稱為微分方程的階數(shù)。如果未知函數(shù)y及其各階導(dǎo)數(shù)都是一次的,則稱它是線性的,否則稱為非線性的。 ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis常微分方程( ODEs 未知函數(shù)是一元函數(shù)) 偏微分方程( PDEs 未知函數(shù)是多元函數(shù)) ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis同一個微分方程,具有不同的初始條

5、件ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis當(dāng)x=0時,y=1,可得c=1特解當(dāng)x=0時,y=1,可得c=-1特解兩邊積分通解ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.1 引言 在高等數(shù)學(xué)中，對于常微分方程的求解，給出了一些典型方程求解析解的基本方法，如可分離變量法、常系數(shù)齊次線性方程的解法、常系數(shù)非齊次線性方程的解法等。但能求解的常微分方程仍然是有限的，大多數(shù)的常微分方程是不可能給出解析解。

6、 ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.1 引言 待求解的問題：一階常微分方程的初值問題 /* Initial-Value Problem */: 解的存在唯一性（“常微分方程”理論）：只要 f (x, y) 在a, b R1 上連續(xù)，且關(guān)于 y 滿足 Lipschitz 條件，即存在與 x, y 無關(guān)的常數(shù) L 使對任意定義在 a, b 上的 y1(x) 和 y2(x) 都成立，則上述IVP存在唯一解。ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-201

7、1學(xué)年課程 數(shù)值分析 Numerical Analysis解析解法：（常微分方程理論）只能求解極少一類常微分方程；實際中給定的問題不一定是解析表達(dá)式，而是函數(shù)表，無法用解析解法。如何求解ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical AnalysisISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical AnalysisISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical An

8、alysis8.2 歐拉(Euler)法 歐拉（Euler）方法是解初值問題的最簡單的數(shù)值方法。初值問題的解y=y(x)代表通過點 的一條稱之為微分方程的積分曲線。積分曲線上每一點 的切線的斜率 等于函數(shù) 在這點的值。 ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical AnalysisEuler法的求解過程是:從初始點P0(即點(x0,y0)出發(fā),作積分曲線y=y(x)在P0點上切線 (其斜率為 ),與x=x1直線相交于P1點(即點(x1,y1),得到y(tǒng)1作為y(x1)的近似值,如上圖所示。過點(x0,y0),以f(x

9、0,y0)為斜率的切線方程為 當(dāng) 時,得 這樣就獲得了P1點的坐標(biāo)。 ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis同樣, 過點P1(x1,y1),作積分曲線y=y(x)的切線交直線x=x2于P2點,切線 的斜率 直線方程為當(dāng) 時,得 ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis當(dāng) 時,得由此獲得了P2的坐標(biāo)。重復(fù)以上過程,就可獲得一系列的點:P1,P1,Pn。對已求得點以 為斜率作直線 取ISCM 2

10、007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis 從圖形上看,就獲得了一條近似于曲線y=y(x) 的折線 。這樣,從x0逐個算出對應(yīng)的數(shù)值解 ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法通常取 (常數(shù)),則Euler法的計算格式 i=0,1,n ( 8.2 ) ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical An

11、alysis8.2 歐拉(Euler)法還可用以下方法推導(dǎo)Euler格式： 數(shù)值微分 數(shù)值積分法對微分方程的離散，可以有多種思路，但最基本的想法是“以直代曲”ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法(1) 用差商近似導(dǎo)數(shù)差分方程初值問題向前Euler方法ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法若用向后差商近似導(dǎo)數(shù)，即向后Euler方法IS

12、CM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法（2）用數(shù)值積分方法ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法若對積分用梯形公式，則得梯形歐拉公式ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis例8.2.1 用歐拉法解初值問題 取步長h=0.2 ,計算過程保留

13、4位小數(shù) 解: h=0.2, 歐拉迭代格式 當(dāng) k=0, x1=0.2時，已知x0=0，y0=1，有 y(0.2)y1=0.21(401)0.8當(dāng) k=1, x2=0.4時，已知x1 =0.2， y1 =0.8，有 y(0.4) y2 =0.20.8(40.20.8)0.6144當(dāng) k=2, x3 =0.6時，已知x2 =0.4， y2 =0.6144，有 y(0.6) y3=0.20.6144(4-0.40.6144)=0.4613 ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Eu

14、ler)法的解作為微分方程初值問題的數(shù)值解，即以差分方程初值問題ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法x0 x1x2x3y0h h h 歐拉折線法ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis解：Eul

15、er公式為當(dāng)h=0.5時ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis當(dāng)h=0.25時ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis00.50.751.010.25h = 0.5h = 0.25ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法歐拉方法的收斂性ISCM 2007,Beij

16、ing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法局部截斷誤差稱為局部截斷誤差I(lǐng)SCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法歐拉方法的收斂性定義 若給定方法的局部截斷誤差滿足則稱該方法是 P 階的，或稱為具有 P 階精度。ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Eul

17、er)法整體截斷誤差I(lǐng)SCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法歐拉方法的收斂性ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis由此知，當(dāng) ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法 注 整體截斷誤差與局部截斷誤差的關(guān)系： ISCM 2007,Beij

18、ing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法 向后歐拉公式隱式歐拉法或向后歐拉法 /* implicit Euler method or backward Euler method*/xn+1點向后差商近似導(dǎo)數(shù)隱式或后退歐拉公式ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法 向后歐拉公式由于未知數(shù) yn+1 同時出現(xiàn)在等式的兩邊，故稱為隱式 /* implicit */ 歐拉

19、公式，而前者稱為顯式 /* explicit */ 歐拉公式。隱式公式不能直接求解，一般需要用Euler顯式公式得到初值，然后用Euler隱式公式迭代求解。因此隱式公式較顯式公式計算復(fù)雜，但穩(wěn)定性好（后面分析）。 隱式歐拉公式中的未知數(shù) yn+1 可通過以下迭代法求解：ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.2 歐拉(Euler)法 向后歐拉公式迭代法求隱式歐拉格式中yn+1的收斂性ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 N

20、umerical Analysis 見上圖， 顯然，這種近似也有一定誤差，如何估計這種誤差y(xn+1) yn+1 ？方法同上，基于Taylor展開估計局部截斷誤差。但是注意，隱式公式中右邊含有f(xn+1 , yn +1 ) ,由于yn +1不準(zhǔn)確，所以不能直接用y (xn+1)代替f(xn+1 , yn +1 ) 設(shè)已知曲線上一點 Pn (xn , yn ),過該點作弦線，斜率為(xn+1 , yn +1 ) 點的方向場f(x,y)。若步長h充分小，可用弦線和垂線x=xn+1的交點近似曲線與垂線的交點。幾何意義xnxn+1PnPn+1xyy(x)ISCM 2007,Beijing Chin

21、a/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis 隱式歐拉法的局部截斷誤差：ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical AnalysisISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis 隱式歐拉法的局部截斷誤差：即隱式歐拉公式具有 1 階精度。 隱式歐拉法的局部截斷誤差：ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分

22、析 Numerical Analysis8.2 歐拉(Euler)法 向后歐拉公式比較歐拉顯式公式和隱式公式及其局部截斷誤差顯式公式隱式公式ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis 若將這兩種方法進(jìn)行算術(shù)平均，即可消除誤差的主要部分/*leading term*/而獲得更高的精度,稱為梯形法 梯形公式 /* trapezoid formula */ 顯、隱式兩種算法的平均注：的確有局部截斷誤差 ， 即梯形公式具有2 階精度，比歐拉方法有了進(jìn)步。但注意到該公式是隱式公式，計算時不得不用到迭代法，

23、其迭代收斂性與歐拉公式相似。ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis例8.2.3 對初值問題 證明用梯形公式求得的近似解為 并證明當(dāng)步長h0時,yn收斂于精確解證明: 解初值問題的梯形公式為 整理成顯式 反復(fù)迭代,得到 ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis公式局部截斷誤差精度顯隱穩(wěn)定性步數(shù)歐拉顯式公式1階顯差單步歐拉隱式公式1階隱好單步梯形公式2階隱好單步歐拉法小結(jié)ISCM 2007,Be

24、ijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.3 改進(jìn)歐拉(Euler)方法ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.3 改進(jìn)歐拉(Euler)方法ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.3 改進(jìn)歐拉(Euler)方法 顯式歐拉公式計算工作量小,但精度低。梯形公式雖提高了精度,但為隱式公式,需用迭代法求解,計算工

25、作量大。綜合歐拉公式和梯形公式便可得到改進(jìn)的歐拉公式。 結(jié)合已有格式的優(yōu)點，以得到計算方便、計算量減少且精度保持的數(shù)值格式ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.3 改進(jìn)歐拉(Euler)方法 先用歐拉公式(8.2)求出一個初步的近似值,稱為預(yù)測值, 它的精度不高, 再用梯形公式對它校正一次,即迭代一次,求得yn+1,稱為校正值, 這種預(yù)測-校正方法稱為改進(jìn)的歐拉公式：稱為Euler公式與梯形公式的預(yù)測校正系統(tǒng)。ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-

26、2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.3 改進(jìn)歐拉(Euler)方法實際計算時，常改寫成以下形式幾何解釋xnxn+1ABPn+1=(A+B)/2歐拉法改進(jìn)歐拉法梯形法ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical AnalysispredictorcorrectorISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.3 改進(jìn)歐拉(Euler)方法 可以證明,改進(jìn)的歐拉公式的精度為二階。這是一種一

27、步顯式格式,它可以表示為嵌套形式。ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis例8.3.1ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical AnalysisISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.3 改進(jìn)歐拉(Euler)方法ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程

28、 數(shù)值分析 Numerical Analysis改進(jìn)歐拉法的算法ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis8.4 單步法的穩(wěn)定性 穩(wěn)定性在微分方程的數(shù)值解法中是一個非常重要的問題。因為微分方程初值問題的數(shù)值方法是用差分格式進(jìn)行計算的，而在差分方程的求解過程中，存在著各種計算誤差，這些計算誤差如舍入誤差等引起的擾動，在傳播過程中，可能會大量積累，對計算結(jié)果的準(zhǔn)確性將產(chǎn)生影響。這就涉及到算法穩(wěn)定性問題。 ISCM 2007,Beijing China/69 鄭州大學(xué)研究生2010-2011學(xué)年課程 數(shù)值分析 Numerical Analysis例：考察初值問題 在區(qū)間0, 0.5上的解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計算數(shù)值解。0.00.10.20.30.40.5精確解改進(jìn)歐拉法 歐拉隱式歐拉顯式 節(jié)點 xi 1.00002.0000 4.00008.0000 1.6000101 3.2000101 1.00002.5000101 6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.9