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

數(shù)值分析

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

數(shù)值分析

NumericalAnalysis第八章鄭州大學(xué)2/66

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

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

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

NumericalAnalysis§8.1引言問題提出

倒葫蘆形狀容器壁上的刻度問題.對于圓柱形狀容器壁上的容積刻度,可以利用圓柱體體積公式其中直徑D為常數(shù).由于體積V與相對于容器底部的任意高度H的函數(shù)關(guān)系明確,因此在容器上可以方便地標(biāo)出容器刻度。對于幾何形狀不是規(guī)則的容器,比如倒葫蘆形狀容器壁上如何標(biāo)出刻度呢?3/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分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的刻度,由微元思想分析可知4/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分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á)式,我們無法求出其解析解.5/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分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ù)都是一次的,則稱它是線性的,否則稱為非線性的。6/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值7/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis常微分方程(ODEs未知函數(shù)是一元函數(shù))

偏微分方程(PDEs未知函數(shù)是多元函數(shù))

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

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

NumericalAnalysis當(dāng)x=0時,y=1,可得c=1時特解當(dāng)x=0時,y=1,可得c=-1時特解兩邊積分通解本章重點討論一階常微分方程的初值問題，9/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分10/66

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

NumericalAnalysis§8.1引言

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

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

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

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

NumericalAnalysis§8.2歐拉(Euler)法推導(dǎo)Euler格式：★Taylor展開法★數(shù)值微分★數(shù)值積分法對微分方程的離散，可以有多種思路，但最基本的想法是“以直代曲”15/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(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項,得16/66鄭州大學(xué)研究生2014-2015學(xué)年課17/66

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

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

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

NumericalAnalysis§8.2歐拉(Euler)法（3）用數(shù)值積分方法19/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值20/66

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

NumericalAnalysis§8.2歐拉(Euler)法若對積分用梯形公式，則得梯形歐拉公式20/66鄭州大學(xué)研究生2011-2012學(xué)年課程21/66

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

NumericalAnalysis§8.2歐拉(Euler)法

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

21/66鄭州大學(xué)研究生2011-2012學(xué)年課程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)。

22/66鄭州大學(xué)研究生2011-2012學(xué)年課程23/66

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

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

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

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

NumericalAnalysis

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

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

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

NumericalAnalysis§8.2歐拉(Euler)法27/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值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.461328/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值29/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

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

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

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

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

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

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

NumericalAnalysis由此知，當(dāng)

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

NumericalAnalysis§8.2歐拉(Euler)法

注

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

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

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

/*implicitEulermethodorbackwardEulermethod*/xn+1點向后差商近似導(dǎo)數(shù)隱式或后退歐拉公式39/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(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

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

NumericalAnalysis

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

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

NumericalAnalysis

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

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

NumericalAnalysis例8.2.3

對初值問題

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

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

∴

整理成顯式

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

∴

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

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

NumericalAnalysis公式局部截斷誤差精度顯隱穩(wěn)定性步數(shù)歐拉顯式公式1階顯差單步歐拉隱式公式1階隱好單步梯形公式2階隱好單步歐拉法小結(jié)46/66鄭州大學(xué)研究生2011-2012學(xué)年課47/66

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

NumericalAnalysis§8.3改進(jìn)歐拉(Euler)方法47/66鄭州大學(xué)研究生2011-2012學(xué)年課48/66

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

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

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

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

結(jié)合已有格式的優(yōu)點，以得到計算方便、計算量減少且精度保持的數(shù)值格式49/66鄭州大學(xué)研究生課程數(shù)值分析Numeri50/66

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

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

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

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

NumericalAnalysis§8.3改進(jìn)歐拉(Euler)方法實際計算時，常改寫成以下形式幾何解釋xnxn+1ABPn+1=(A+B)/2歐拉法梯形法改進(jìn)歐拉法51/66鄭州大學(xué)研究生2011-2012學(xué)年課52/66

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

NumericalAnalysispredictorcorrector52/66鄭州大學(xué)研究生2011-2012學(xué)年課53/66

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

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

可以證明,改進(jìn)的歐拉公式的精度為二階。這是一種一步顯式格式,它可以表示為嵌套形式。53/66鄭州大學(xué)研究生2011-2012學(xué)年課54/66

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

NumericalAnalysis例8.3.154/66鄭州大學(xué)研究生2011-2012學(xué)年課55/66

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

NumericalAnalysis55/66鄭州大學(xué)研究生2011-2012學(xué)年課56/66

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

NumericalAnalysis§8.3改進(jìn)歐拉(Euler)方法56/66鄭州大學(xué)研究生2011-2012學(xué)年課57/66

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

NumericalAnalysis改進(jìn)歐拉法的算法57/66鄭州大學(xué)研究生2011-2012學(xué)年課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)定性問題。

58/66鄭州大學(xué)研究生2011-2012學(xué)年課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.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107§8.4單步法的穩(wěn)定性59/66鄭州大學(xué)研究生2011-2012學(xué)年課60/66

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

NumericalAnalysis§8.4單步法的穩(wěn)定性60/66鄭州大學(xué)研究生2011-2012學(xué)年課61/66

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

NumericalAnalysis定義若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計算中都逐步衰減，則稱該算法是絕對穩(wěn)定的/*absolutelystable*/。一般分析某算法的穩(wěn)定性時，為簡單起見，只考慮模型方程或試驗方程/*testequation*/§8.4單步法的穩(wěn)定性61/66鄭州大學(xué)研究生2011-2012學(xué)年課62/66

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

NumericalAnalysis記數(shù)值誤差為：引進(jìn)試驗方程：

62/66鄭州大學(xué)研究生2011-2012學(xué)年課63/66

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

NumericalAnalysis向前歐拉公式的穩(wěn)定性試驗方程的歐拉公式若每步計算有舍入誤差，則

63/66鄭州大學(xué)研究生2011-2012學(xué)年課64/66

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

NumericalAnalysis當(dāng)λ為復(fù)數(shù)時當(dāng)λ為實數(shù)時64/66鄭州大學(xué)研究生2011-2012學(xué)年課65/66

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

NumericalAnalysis本章小結(jié)

本章介紹了常微分方程初值問題的基本數(shù)值解法。包括單步法和多步法。單步法主要有歐拉法、改進(jìn)歐拉法和龍格—庫塔方法。多步法是亞當(dāng)姆斯法。它們都是基于把一個連續(xù)的定解問題離散化為一個差分方程來求解，是一種步進(jìn)式的方法。用多步法求常微分方程的數(shù)值解可獲得較高的精度。實際應(yīng)用時，選擇合適的算法有一定的難度，既要考慮算法的簡易性和計算量，又要考慮截斷誤差和收斂性、穩(wěn)定性。65/66鄭州大學(xué)研究生2011-2012學(xué)年課數(shù)值分析

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

數(shù)值分析

NumericalAnalysis第八章鄭州大學(xué)67/66

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

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

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

NumericalAnalysis§8.1引言問題提出

倒葫蘆形狀容器壁上的刻度問題.對于圓柱形狀容器壁上的容積刻度,可以利用圓柱體體積公式其中直徑D為常數(shù).由于體積V與相對于容器底部的任意高度H的函數(shù)關(guān)系明確,因此在容器上可以方便地標(biāo)出容器刻度。對于幾何形狀不是規(guī)則的容器,比如倒葫蘆形狀容器壁上如何標(biāo)出刻度呢?3/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分69/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的刻度,由微元思想分析可知4/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分70/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á)式,我們無法求出其解析解.5/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分71/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ù)都是一次的,則稱它是線性的,否則稱為非線性的。6/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值72/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis常微分方程(ODEs未知函數(shù)是一元函數(shù))

偏微分方程(PDEs未知函數(shù)是多元函數(shù))

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

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

NumericalAnalysis當(dāng)x=0時,y=1,可得c=1時特解當(dāng)x=0時,y=1,可得c=-1時特解兩邊積分通解本章重點討論一階常微分方程的初值問題，9/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分75/66

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

NumericalAnalysis§8.1引言

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

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

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

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

NumericalAnalysis§8.2歐拉(Euler)法推導(dǎo)Euler格式：★Taylor展開法★數(shù)值微分★數(shù)值積分法對微分方程的離散，可以有多種思路，但最基本的想法是“以直代曲”15/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值81/66

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

NumericalAnalysis81/66

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

NumericalAnalysis§8.2歐拉(Euler)法(1)Taylor展開法方程初值問題Euler公式設(shè)給定等距剖分當(dāng)步長h充分小時，略去h2項,得16/66鄭州大學(xué)研究生2014-2015學(xué)年課82/66

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

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

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

NumericalAnalysis§8.2歐拉(Euler)法（3）用數(shù)值積分方法19/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值85/66

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

NumericalAnalysis§8.2歐拉(Euler)法若對積分用梯形公式，則得梯形歐拉公式20/66鄭州大學(xué)研究生2011-2012學(xué)年課程86/66

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

NumericalAnalysis§8.2歐拉(Euler)法

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

21/66鄭州大學(xué)研究生2011-2012學(xué)年課程87/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)。

22/66鄭州大學(xué)研究生2011-2012學(xué)年課程88/66

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

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

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

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

NumericalAnalysis

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

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

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

NumericalAnalysis§8.2歐拉(Euler)法27/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值93/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.461328/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值94/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

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

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

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

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

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

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

NumericalAnalysis由此知，當(dāng)

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

NumericalAnalysis§8.2歐拉(Euler)法

注

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

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

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

/*implicitEulermethodorbackwardEulermethod*/xn+1點向后差商近似導(dǎo)數(shù)隱式或后退歐拉公式39/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值105/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

可通過以下迭代法求解：40/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值106/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)41/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值107/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis

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

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

NumericalAnalysis

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

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

NumericalAnalysis例8.2.3

對初值問題

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

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

∴

整理成顯式

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

∴

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

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

NumericalAnalysis公式局部截斷誤差精度顯隱穩(wěn)定性步數(shù)歐拉顯式公式1階顯差單步歐拉隱式公式1階隱好單步梯形公式2階隱好單步歐拉法小結(jié)46/66鄭州大學(xué)研究生2011-2012學(xué)年課112/66

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

NumericalAnalysis§8.3改進(jìn)歐拉(Euler)方法47/66鄭州大學(xué)研究生2011-2012學(xué)年課113/66

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

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

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

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

結(jié)合已有格式的優(yōu)點，以得到計算方便、計算量減少且精度保持的數(shù)值格式49/66鄭州大學(xué)研究生課程數(shù)值分析Numeri115/66

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

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

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

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

NumericalAnalysis§8.3改進(jìn)歐拉(Euler)方法實際計算時，常改寫成以下形式幾何解釋xnxn+1ABPn+1=(A+B)/2歐拉法梯形法改進(jìn)歐拉法51/66鄭州大學(xué)研究生2011-2012學(xué)年課117/66

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

NumericalAnalysispredictorcorrector52/66