數值分析

NumericalAnalysis第八章常微分方程數值解法鄭州大學研究生課程（2014-2015學年第一學期）

2/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis第八章常微分方程數值解法

§8.1引言§8.2歐拉(Euler)法§8.3改進歐拉(Euler)方法§8.4單步法的穩定性3/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.1引言問題提出

倒葫蘆形狀容器壁上的刻度問題.對于圓柱形狀容器壁上的容積刻度,可以利用圓柱體體積公式其中直徑D為常數.由于體積V與相對于容器底部的任意高度H的函數關系明確,因此在容器上可以方便地標出容器刻度。對于幾何形狀不是規則的容器,比如倒葫蘆形狀容器壁上如何標出刻度呢?4/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.1引言下表是經過測量得到部分容器高度與直徑的關系.H00.20.40.60.81.0D00.110.260.561.041.17根據上表的數據,可以擬合出倒葫蘆形狀容器的圖,建立如圖所示的坐標軸后,問題即為:如何根據任意高度x標出容器體積V的刻度,由微元思想分析可知5/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.1引言其中x表示高度,直徑D是高度x的函數,記為D(x),因此得到如下微分方程初值問題只要求解上述方程,就可求出體積V與高度x之間的函數關系,從而可標出容器壁上容積的刻度,但問題是函數D(x)無解析表達式,我們無法求出其解析解.6/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.1引言

包含自變量x、未知函數y及未知函數的導數或微分的方程稱為微分方程。在微分方程中,自變量的個數只有一個,稱為常微分方程。自變量的個數為兩個或兩個以上的微分方程叫偏微分方程。微分方程中出現的未知函數最高階導數的階數稱為微分方程的階數。如果未知函數y及其各階導數都是一次的,則稱它是線性的,否則稱為非線性的。7/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis常微分方程(ODEs

未知函數是一元函數)

偏微分方程(

PDEs

未知函數是多元函數)

8/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis同一個微分方程,具有不同的初始條件微分方程的定解條件：9/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis當x=0時,y=1,可得c=1時特解當x=0時,y=1,可得c=-1時特解兩邊積分通解本章重點討論一階常微分方程的初值問題，10/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.1引言

在高等數學中，對于常微分方程的求解，給出了一些典型方程求解析解的基本方法，如可分離變量法、常系數齊次線性方程的解法、常系數非齊次線性方程的解法等。但能求解的常微分方程仍然是有限的，大多數的常微分方程是不可能給出解析解。11/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.1引言

待求解的問題：一階常微分方程的初值問題/*Initial-ValueProblem*/:

解的存在唯一性（“常微分方程”理論）：只要f(x,y)在[a,b]R1上連續，且關于y

滿足Lipschitz

條件，即存在與x,y無關的常數L

使則上述IVP存在唯一解。12/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis解析解法：（常微分方程理論）只能求解極少一類常微分方程；實際中給定的問題不一定是解析表達式，而是函數表，無法用解析解法。數值解法：遞推法如何求解13/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis14/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis記號：15/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法推導Euler格式：★Taylor展開法★數值微分★數值積分法對微分方程的離散，可以有多種思路，但最基本的想法是“以直代曲”16/66

鄭州大學研究生2014-2015學年課程數值分析

NumericalAnalysis16/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法(1)Taylor展開法方程初值問題Euler公式設給定等距剖分當步長h充分小時，略去h2項,得17/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法(2)用差商近似導數差分方程初值問題向前Euler方法18/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法若用向后差商近似導數，即向后Euler方法19/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法（3）用數值積分方法20/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法若對積分用梯形公式，則得梯形歐拉公式21/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法

歐拉（Euler）方法是解初值問題的最簡單的數值方法。初值問題的解y=y(x)代表通過點的一條稱之為微分方程的積分曲線。積分曲線上每一點的切線的斜率等于函數在這點的值。

22/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysisEuler法的求解過程是:從初始點P0(即點(x0,y0))出發,作積分曲線y=y(x)在P0點上切線(其斜率為

),與x=x1直線相交于P1點(即點(x1,y1),得到y1作為y(x1)的近似值,如上圖所示。過點(x0,y0),以f(x0,y0)為斜率的切線方程為

當時,得這樣就獲得了P1點的坐標。

23/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis同樣,過點P1(x1,y1),作積分曲線y=y(x)的切線交直線x=x2于P2點,切線的斜率直線方程為當時,得24/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis當時,得由此獲得了P2的坐標。重復以上過程,就可獲得一系列的點:P1,P1,…,Pn。對已求得點以為斜率作直線

取25/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis

從圖形上看,就獲得了一條近似于曲線y=y(x)

的折線。這樣,從x0逐個算出對應的數值解26/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法x0x1x2x3y0hhh歐拉折線法27/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法28/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis例8.2.1

用歐拉法解初值問題

取步長h=0.2,計算過程保留4位小數

解:h=0.2,歐拉迭代格式

當n=0時，已知x0=0，y0=1，有

y(x1)=y(0.2)

y1=0.2×1(4－0×1)＝0.8當n=1時，已知x1=0.2，y1=0.8，有

y(0.4)y2

=0.2×0.8×(4－0.2×0.8)＝0.6144當n=2,時，已知x2=0.4，y2=0.6144，有

y(0.6)y3=0.2×0.6144×(4-0.4×0.6144)=0.461329/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis解：Euler公式為當h=0.5時30/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis當h=0.25時31/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis00.50.751.010.25h=0.5h=0.2532/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法歐拉方法的收斂性假設第n步是準確的33/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法局部截斷誤差稱為局部截斷誤差34/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法歐拉方法的收斂性定義若給定方法的局部截斷誤差滿足則稱該方法是P階的，或稱為具有P階精度。35/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法整體截斷誤差36/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法歐拉方法的收斂性37/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis由此知，當

38/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法

注

整體截斷誤差與局部截斷誤差的關系：

39/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法向后歐拉公式隱式歐拉法或向后歐拉法

/*implicitEulermethodorbackwardEulermethod*/xn+1點向后差商近似導數隱式或后退歐拉公式40/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法向后歐拉公式由于未知數yn+1

同時出現在等式的兩邊，故稱為隱式/*implicit*/

歐拉公式，而前者稱為顯式/*explicit*/歐拉公式。隱式公式不能直接求解，一般需要用Euler顯式公式得到初值，然后用Euler隱式公式迭代求解。因此隱式公式較顯式公式計算復雜，但穩定性好（后面分析）。

隱式歐拉公式中的未知數yn+1

可通過以下迭代法求解：41/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis

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

yn+1

？方法同上，基于Taylor展開估計局部截斷誤差。但是注意，隱式公式中右邊含有f(xn+1

,yn

+1),由于yn

+1不準確，所以不能直接用y'(xn+1)代替f(xn+1

,yn

+1)設已知曲線上一點Pn

(xn

,yn),過該點作弦線，斜率為(xn+1

,yn

+1)點的方向場f(x,y)。若步長h充分小，可用弦線和垂線x=xn+1的交點近似曲線與垂線的交點。幾何意義xnxn+1PnPn+1xyy(x)42/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis

隱式歐拉法的局部截斷誤差：即隱式歐拉公式具有1階精度。43/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.2歐拉(Euler)法向后歐拉公式比較歐拉顯式公式和隱式公式及其局部截斷誤差顯式公式隱式公式44/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis

若將這兩種方法進行算術平均，即可消除誤差的主要部分/*leadingterm*/而獲得更高的精度,稱為梯形法梯形公式/*trapezoidformula*/—顯、隱式兩種算法的平均注：的確有局部截斷誤差，即梯形公式具有2

階精度，比歐拉方法有了進步。但注意到該公式是隱式公式，計算時不得不用到迭代法，其迭代收斂性與歐拉公式相似。45/66鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis例8.2.3

對初值問題

證明用梯形公式求得的近似解為

并證明當步長h0時,yn收斂于精確解證明:解初值問題的梯形公式為∵

∴

整理成顯式

反復迭代,得到∵

∴

46/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis公式局部截斷誤差精度顯隱穩定性步數歐拉顯式公式1階顯差單步歐拉隱式公式1階隱好單步梯形公式2階隱好單步歐拉法小結47/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.3改進歐拉(Euler)方法48/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.3改進歐拉(Euler)方法49/66鄭州大學研究生課程數值分析

NumericalAnalysis§8.3改進歐拉(Euler)方法

顯式歐拉公式計算工作量小,但精度低。梯形公式雖提高了精度,但為隱式公式,需用迭代法求解,計算工作量大。綜合歐拉公式和梯形公式便可得到改進的歐拉公式。

結合已有格式的優點，以得到計算方便、計算量減少且精度保持的數值格式50/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.3改進歐拉(Euler)方法

先用歐拉公式(8.2)求出一個初步的近似值,稱為預測值,它的精度不高,再用梯形公式對它校正一次,即迭代一次,求得yn+1,稱為校正值,這種預測-校正方法稱為改進的歐拉公式：稱為Euler公式與梯形公式的預測—校正系統。51/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.3改進歐拉(Euler)方法實際計算時，常改寫成以下形式幾何解釋xnxn+1ABPn+1=(A+B)/2歐拉法梯形法改進歐拉法52/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysispredictorcorrector53/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.3改進歐拉(Euler)方法

可以證明,改進的歐拉公式的精度為二階。這是一種一步顯式格式,它可以表示為嵌套形式。54/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis例8.3.155/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis56/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.3改進歐拉(Euler)方法57/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis改進歐拉法的算法58/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis§8.4單步法的穩定性穩定性:誤差在以后各步的計算中不會無限制擴大.穩定性在微分方程的數值解法中是一個非常重要的問題。因為微分方程初值問題的數值方法是用差分格式進行計算的，而在差分方程的求解過程中，存在著各種計算誤差，這些計算誤差如舍入誤差等引起的擾動，在傳播過程中，可能會大量積累，對計算結果的準確性將產生影響。這就涉及到算法穩定性問題。

59/66

鄭州大學研究生2011-2012學年課程數值分析

NumericalAnalysis例：考察初值問題在區間[0,0.5]上的解。分別用歐拉顯、隱式格式和改進的歐拉格式計算數值解。0.00.10.20.30.40.5精確解改進歐拉法

歐拉隱式歐拉顯式

節點xi

1.00002.00004.00008.00001.6000101

3.2000101

1.00002.5000101

6.250