常微分方程數(shù)值解法課件-002

第7章常微分方程數(shù)值解法7.1引言7.2簡(jiǎn)單的數(shù)值方法與基本概念（歐拉法）7.3-庫塔方法7.5方程組和高階方程7.1引言

科學(xué)技術(shù)中常常需要求解常微分方程的定解問題.這類問題最簡(jiǎn)單的形式，是本章將著重考察的一階方程的初值問題

我們知道，只有f(x,y)適當(dāng)光滑—譬如關(guān)于y滿足萊布尼茲(Lipschitz)條件理論上就可以保證初值問題的解y＝f(x)存在并且唯一.

雖然求解常微分方程有各種各樣的解析方法，但解析方法只能用來求解一些特殊類型的方程，實(shí)際問題中歸結(jié)出來的微分方程主要靠數(shù)值解法.

所謂數(shù)值解法,就是尋求解y(x)在一系列離散節(jié)點(diǎn)上的近似值y1,y2,,yn,yn+1,.相鄰兩個(gè)節(jié)點(diǎn)的間距hn=xn+1-xn稱為步長.今后如不特別說明，總是假定hi=h(i=1,2,)為常數(shù),這時(shí)節(jié)點(diǎn)為xn=x0+nh(i=0,1,2,)(等距節(jié)點(diǎn)).截去最后一項(xiàng),可得這就是著名的(顯式)歐拉(Euler)公式.若初值y0已知，則依公式(2.1)可逐次逐步算出各點(diǎn)數(shù)值解.即7.2簡(jiǎn)單的數(shù)值方法與基本概念7.2.1歐拉法與后退歐拉法用向前差商代替導(dǎo)數(shù)

例1

用歐拉公式求解初值問題

解取步長h=0.1，歐拉公式的具體形式為其中xn=nh=0.1n(n=0,1,,10),已知y0=1,由此式可得依次計(jì)算下去，部分計(jì)算結(jié)果見下表.與準(zhǔn)確解相比，可看出歐拉公式的計(jì)算結(jié)果精度很差.

xn

歐拉公式數(shù)值解yn準(zhǔn)確解y(xn)

誤差0.20.40.60.81.01.1918181.3582131.5089661.6497831.7847701.1832161.3416411.4832401.6124521.7320510.0086020.0165720.0257260.0373310.052719

歐拉公式具有明顯的幾何意義,就是用折線近似代替方程的解曲線，因而常稱公式(2.1)為歐拉折線法.

還可以通過幾何直觀來考察歐拉方法的精度.假設(shè)yn=y(xn),即頂點(diǎn)Pn落在積分曲線y=y(x)上，那么，按歐拉方法做出的折線PnPn+1便是y=y(x)過點(diǎn)Pn的切線.從圖形上看,這樣定出的頂點(diǎn)Pn+1顯著地偏離了原來的積分曲線，可見歐拉方法是相當(dāng)粗糙的.

為了分析計(jì)算公式的精度，通?？捎锰├照归_將y(xn+1)在xn處展開，則有在yn=y(xn)的前提下，f(xn,yn)=f(xn,y(xn))=y(xn).于是可得歐拉法(2.1)的公式誤差為稱為此方法的局部截?cái)嗾`差.截去最后一項(xiàng),可得這就是著名的(隱式)后退歐拉(Euler)公式.若初值y0已知，則依公式(2.5)可逐次逐步算出各點(diǎn)數(shù)值解.即用向后差商代替導(dǎo)數(shù)

后退的歐拉公式與歐拉公式有著本質(zhì)的區(qū)別,后者是關(guān)于yn+1的一個(gè)直接計(jì)算公式，這類公式稱作是顯式的；前者公式的右端含有未知的yn+1，它實(shí)際上是關(guān)于yn+1的一個(gè)函數(shù)方程,這類方程稱作是隱式的.

顯式與隱式兩類方法各有特點(diǎn)，考了到數(shù)值穩(wěn)定性等其他因素，人們有時(shí)需要選用隱式方法，但使用顯式算法遠(yuǎn)比隱式方便.

隱式方程通常用迭代法求解，而迭代過程的實(shí)質(zhì)是逐步顯式化.7.2.2梯形方法

為得到比歐拉法精度高的計(jì)算公式，將歐拉公式和后退的歐拉公式進(jìn)行算數(shù)平均，可得稱為梯形公式.

梯形公式是隱式單步法，用迭代法求解，同后退的歐拉方法一樣，仍用歐拉法提供迭代初值。7.2.3單步法的局部截?cái)嗾`差與階

初值問題(1.1),(1.2)的單步法可用一般形式表示為其中多元函數(shù)與f(x,y

)有關(guān)，當(dāng)含有yn+1時(shí)，方法是隱式的，若不含yn+1則為顯式方法，所以顯式單步法可表示為(x,y,h)稱為增量函數(shù)，例如對(duì)歐拉法(2.1)有它的局部截?cái)嗾`差(2.3)已由給出,對(duì)一般顯式單步法則可如下定義.

定義1

設(shè)y(x)是初值問題(1.1),(1.2)的準(zhǔn)確解,稱為顯式單步法(2.10)的局部截?cái)嗾`差.

Tn+1之所以稱為局部的，是假設(shè)在xn前各步?jīng)]有誤差.當(dāng)yn=y(xn)時(shí)，計(jì)算一步，則有

所以，局部截?cái)嗾`差可理解為用方法(2.10)計(jì)算一步的誤差，也即公式(2.10)中用準(zhǔn)確解y(x)代替數(shù)值解產(chǎn)生的公式誤差.根據(jù)定義,顯然歐拉法的局部截?cái)嗾`差為顯然Tn+1=O(h2).一般情形的定義如下

定義2

設(shè)y(x)是初值問題的準(zhǔn)確解，若存在最大整數(shù)p使顯式單步法(2.10)的局部截?cái)嗾`差滿足則稱方法(2.10)具有p階精度.若將(2.10)展開式寫成則稱為局部截?cái)嗾`差主項(xiàng).

以上定義對(duì)隱式單步法(2.9)也是適用的.例如，對(duì)后退歐拉法(2.5)其局部截?cái)嗾`差為這里p=1是1階方法，局部截?cái)嗾`差主項(xiàng)為同樣對(duì)梯形法(2.7)有所以梯形方法(2.7)是二階的.其局部截?cái)嗾`差主項(xiàng)為7.2.4改進(jìn)的歐拉公式我們看到，梯形方法雖然提高了精度，但其算法復(fù)雜，在應(yīng)用迭代公式(2.9)進(jìn)行實(shí)際計(jì)算時(shí)，每迭代一次，都要重新計(jì)算函數(shù)f(x,y

)的值，而迭代又要反復(fù)進(jìn)行若干次，計(jì)算量很大，而且往往難以預(yù)測(cè).為了控制計(jì)算量，通常只迭代一兩次就轉(zhuǎn)入下一步的計(jì)算，這就簡(jiǎn)化了算法.具體地說，我們先用歐拉公式求得一個(gè)初步的近似值，稱之為預(yù)測(cè)值，此預(yù)測(cè)值的精度可能很差，再用梯形公式(2.7)將它校正一次，即按(2.8)式迭代一次，這個(gè)結(jié)果稱之為校正值.這樣建立的預(yù)測(cè)—校正系統(tǒng)通常稱為改進(jìn)的歐拉公式:或表示為下列平均化形式(2.13)預(yù)測(cè)校正

例2

用改進(jìn)的歐拉法解例1中的初值問題(2.2).

解仍取步長h=0.1,改進(jìn)的歐拉公式為部分計(jì)算結(jié)果見下表

xnyn

誤差

xnyn

誤差00.20.41.1.1840961.34336000.0000880.0017190.60.81.01.4859561.6164761.7378690.0027160.0040240.005818同例1中的歐拉法的計(jì)算結(jié)果比較，明顯改善了精度.

xn

歐拉公式

改進(jìn)歐拉公式

yn

誤差

yn

誤差0.00.20.40.60.81.01.1.1918181.3582131.5089661.6497831.784770

00.0086020.016572

0.0257260.0373310.05271911.1840961.3433601.4859561.6164761.737869

00.0000880.0017190.0027160.0040240.0058187.3龍格—庫塔方法

對(duì)許多實(shí)際問題來說，歐拉公式與改進(jìn)歐拉公式精度還不能滿足要求，為此從另一個(gè)角度來分析這兩個(gè)公式的特點(diǎn)，從而探索一條構(gòu)造高精度方法的途徑.

7.3.1顯式龍格—庫塔法的一般形式

上節(jié)給出了顯式單步法的表達(dá)式(2.10),其局部截?cái)嗾`差為O(hp+1)，對(duì)歐拉法Tn+1=O(h2)，即方法為p=1階，若用改進(jìn)歐拉法(2.13)，它可表為此時(shí)增量函數(shù)為一般說來，點(diǎn)數(shù)r越多，精度越高，上式右端相當(dāng)于增量函數(shù)(xn,yn,h)，為得到便于計(jì)算的顯式方法，可類似于改進(jìn)歐拉法(3.1),(3.2)，將公式表示為其中這里ci,ai,bij均為常數(shù).(3.4)和(3.5)稱為r級(jí)顯式龍格-庫塔(Runge-Kutta)法,簡(jiǎn)稱R-K方法.7.3.2二階顯式R-K方法

對(duì)r=2的R-K方法，由(3.4),(3.5)式可得如下計(jì)算公式這里c1,c2,a2,b21均為待定常數(shù)，我們希望適當(dāng)選取這些系數(shù)，使公式階數(shù)p盡量高.

則由此可以看出在改進(jìn)的歐拉公式中相當(dāng)于取(xn,yn),(xn+1,yn+1)兩點(diǎn)處斜率的平均值，近似代替平均斜率，其精度比歐拉公式提高了.c1=c2=1/2,a2=b21=a=1時(shí)，即為改進(jìn)的歐拉公式：稱為變形的歐拉公式,也可以表示為:

如取a=1，則c1=0,c2=1,a2=b21=1/2.得計(jì)算公式

對(duì)r=2的R-K公式能否使局部誤差提高到O(h4)?為此需把k2多展開一項(xiàng)，展開式中的項(xiàng)是不能通過選擇參數(shù)削掉的，實(shí)際上要使h3

的項(xiàng)為零，需增加3個(gè)方程，要確定4個(gè)參數(shù)c1,c2,λ2,μ21，這是不可能的.故r＝2的顯式R-K方法的階只能是p=2，而不能得到三階公式.

四階R-K方法的每一步需要計(jì)算四次函數(shù)值f，可以證明其局部截?cái)嗾`差為O(h5).

7.3.3標(biāo)準(zhǔn)四階顯式R-K方法

然而值得指出的是，龍格-庫塔方法的推導(dǎo)基于泰勒展開方法，因而它要求所求的解具有較好的光滑性質(zhì).反之，如果解的光滑性差，那么，使用龍格-庫塔方法求得的數(shù)值解，其精度可能反而不如改進(jìn)的歐拉方法.實(shí)際計(jì)算時(shí)，我們應(yīng)當(dāng)針對(duì)問題的具體特點(diǎn)選擇合適的算法.7.3.4變步長的龍格-庫塔方法

單從每一步看，步長越小，截?cái)嗾`差就越小，但隨著步長的縮小，在一定求解范圍內(nèi)所要完成的步數(shù)就增加了.步數(shù)的增加不但引起計(jì)算量的增大，而且可能導(dǎo)致舍入誤差的嚴(yán)重積累.因此同積分的數(shù)值計(jì)算一樣，微分方程的數(shù)值解法也有個(gè)選擇步長的問題.

在選擇步長時(shí)，需要考慮兩個(gè)問題：

1.怎樣衡量和檢驗(yàn)計(jì)算結(jié)果的精度？

2.如何依據(jù)所獲得的精度處理步長？我們考察四階R-K公式(3.13)，從節(jié)點(diǎn)xn出發(fā)，先以h為步長求出一個(gè)近似值，記為，由于公式的局部截?cái)嗾`差為O(h5)，故然后將步長折半，即取為步長,從xn跨兩步到xn+1，再求得一個(gè)近似值，每跨一步的局部截?cái)嗾`差是

，因此有比較(3.14)式和(3.15)式我們看到，步長折半后，誤差大約減少到1/16，即有由此易得下列事后估計(jì)式這樣，我們可以通過檢查步長，折半前后兩次計(jì)算結(jié)果的偏差來判定所選的步長是否合適，具體地說，將區(qū)分以下兩種情況處理：這種通過加倍或折半處理步長的方法稱為變步長方法.表面上看，為了選擇步長，每一步的計(jì)算量增加了，但總體考慮往往是合算的.

1.對(duì)于給定的精度ε,如果Δ>ε,我們反復(fù)將步長折半計(jì)算,直至Δ<ε為止,這時(shí)取最終得到的作為結(jié)果；

2.如果Δ<ε，我們將反復(fù)將步長作加倍計(jì)算，直至Δ>ε為止，這時(shí)再將步長折半計(jì)算一次，就得到所要的結(jié)果.7.5方程組和高階方程7.5.1一階方程組

前面我們研究了單個(gè)方程y=f

的數(shù)值解，只要把y和f理解為向量，那么，所提供的各種計(jì)算公式即可應(yīng)用到一階方程組的情形.

考察一階方程組的初值問題，初始條件給為若采用向量的記號(hào)，記(向量)則上述方程組的初值問題可表示為求解這一初值問題的四階龍格-庫塔公式為(向量)式中(向量)或表示為(分量)其中(分量)這里yin是第i個(gè)因變量yi(x)在節(jié)點(diǎn)xn=x0+nh的近似值.為了幫助理解這一公式的計(jì)算過程，我們?cè)倏疾靸蓚€(gè)方程的特殊情形這時(shí)四階龍格-庫塔公式具有形式其中