數(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引言問題提出

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

NumericalAnalysis§8.1引言下表是經(jīng)過測(cè)量得到部分容器高度與直徑的關(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ù)或微分的方程稱為微分方程。在微分方程中,自變量的個(gè)數(shù)只有一個(gè),稱為常微分方程。自變量的個(gè)數(shù)為兩個(gè)或兩個(gè)以上的微分方程叫偏微分方程。微分方程中出現(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同一個(gè)微分方程,具有不同的初始條件微分方程的定解條件：9/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

NumericalAnalysis§8.1引言

在高等數(shù)學(xué)中，對(duì)于常微分方程的求解，給出了一些典型方程求解析解的基本方法，如可分離變量法、常系數(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解析解法：（常微分方程理論）只能求解極少一類常微分方程；實(shí)際中給定的問題不一定是解析表達(dá)式，而是函數(shù)表，無法用解析解法。數(shù)值解法：遞推法如何求解13/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

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

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

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

NumericalAnalysis16/66

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

NumericalAnalysis§8.2歐拉(Euler)法(1)Taylor展開法方程初值問題Euler公式設(shè)給定等距剖分當(dāng)步長(zhǎng)h充分小時(shí)，略去h2項(xiàng),得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)法若對(duì)積分用梯形公式，則得梯形歐拉公式21/66

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

NumericalAnalysis§8.2歐拉(Euler)法

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

22/66

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

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

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

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

23/66

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

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

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

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

NumericalAnalysis

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

的折線。這樣,從x0逐個(gè)算出對(duì)應(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

用歐拉法解初值問題

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

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

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

y(x1)=y(0.2)

y1=0.2×1(4－0×1)＝0.8當(dāng)n=1時(shí)，已知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,時(shí)，已知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時(shí)30/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis當(dāng)h=0.25時(shí)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)法局部截?cái)嗾`差稱為局部截?cái)嗾`差34/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

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

NumericalAnalysis§8.2歐拉(Euler)法整體截?cái)嗾`差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)法

注

整體截?cái)嗾`差與局部截?cái)嗾`差的關(guān)系：

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

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

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

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

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

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

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

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

NumericalAnalysis

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

yn+1

？方法同上，基于Taylor展開估計(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),過該點(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)42/66鄭州大學(xué)研究生2011-2012學(xué)年課程數(shù)值分析

NumericalAnalysis

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

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

NumericalAnalysis

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

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

NumericalAnalysis例8.2.3

對(duì)初值問題

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

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

∴

整理成顯式

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

∴

46/66

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

NumericalAnalysis公式局部截?cái)嗾`差精度顯隱穩(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ì)算工作量小,但精度低。梯形公式雖提高了精度,但為隱式公式,需用迭代法求解,計(jì)算工作量大。綜合歐拉公式和梯形公式便可得到改進(jìn)的歐拉公式。

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

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

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

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

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

NumericalAnalysis§8.3改進(jìn)歐拉(Euler)方法實(shí)際計(jì)算時(shí)，常改寫成以下形式幾何解釋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)定性:誤差在以后各步的計(jì)算中不會(huì)無限制擴(kuò)大.穩(wěn)定性在微分方程的數(shù)值解法中是一個(gè)非常重要的問題。因?yàn)槲⒎址匠坛踔祮栴}的數(shù)值方法是用差分格式進(jìn)行計(jì)算的，而在差分方程的求解過程中，存在著各種計(jì)算誤差，這些計(jì)算誤差如舍入誤差等引起的擾動(dòng)，在傳播過程中，可能會(huì)大量積累，對(duì)計(jì)算結(jié)果的準(zhǔn)確性將產(chǎn)生影響。這就涉及到算法穩(wěn)定性問題。

59/66

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

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

歐拉隱式歐拉顯式

節(jié)點(diǎn)xi

1.00002.00004.00008.00001.6000101

3.2000101

1.00002.5000101

6.250