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

數(shù)值分析

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

2/66

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

NumericalAnalysis§8.1引言問題提出

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

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

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

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis§8.1引言

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

NumericalAnalysis常微分方程(ODEs

未知函數(shù)是一元函數(shù))

偏微分方程(

PDEs

未知函數(shù)是多元函數(shù))

8/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis§8.1引言

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

NumericalAnalysis§8.1引言

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

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

滿足Lipschitz

條件，即存在與x,y無關(guān)的常數(shù)L

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

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

NumericalAnalysis14/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

鄭州大學(xué)研究生2014-2015學(xué)年課程數(shù)值分析

NumericalAnalysis16/66

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

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

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis§8.2歐拉(Euler)法

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

22/66

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

當(dāng)時,得這樣就獲得了P1點的坐標(biāo)。

23/66

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

取25/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis

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

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

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

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

NumericalAnalysis例8.2.1

用歐拉法解初值問題

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

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

當(dāng)n=0時，已知x0=0，y0=1，有

y(x1)=y(0.2)

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

y(0.4)y2

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

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

NumericalAnalysis解：Euler公式為當(dāng)h=0.5時30/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis當(dāng)h=0.25時31/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

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

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

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

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

NumericalAnalysis由此知，當(dāng)

38/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis§8.2歐拉(Euler)法

注

整體截斷誤差與局部截斷誤差的關(guān)系：

39/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

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

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

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

隱式歐拉公式中的未知數(shù)yn+1

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

NumericalAnalysis

見上圖，顯然，這種近似也有一定誤差，如何估計這種誤差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)42/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis

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

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

NumericalAnalysis

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

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

NumericalAnalysis例8.2.3

對初值問題

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

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

∴

整理成顯式

反復(fù)迭代,得到∵

∴

46/66

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis公式局部截斷誤差精度顯隱穩(wěn)定性步數(shù)歐拉顯式公式1階顯差單步歐拉隱式公式1階隱好單步梯形公式2階隱好單步歐拉法小結(jié)47/66

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

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

結(jié)合已有格式的優(yōu)點，以得到計算方便、計算量減少且精度保持的數(shù)值格式50/66

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysispredictorcorrector53/66

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis例8.3.155/66

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis56/66

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis改進(jìn)歐拉法的算法58/66

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

59/66

鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

歐拉隱式歐拉顯式

節(jié)點xi

1.00002.00004.00008.00001.6000101

3.2000101

1.00002.5000101

6.250