計(jì)算方法第五章1-3節(jié)

第五章常微分方程初值問(wèn)題數(shù)值解法15.1

引言本章要著重考察的一階方程的初值問(wèn)題（1.1）（1.2）只要函數(shù)適當(dāng)光滑——譬如關(guān)于滿足利普希茨(Lipschitz）條件（1.3）理論上就可以保證初值問(wèn)題（1.1），（1.2）的解存在并且惟一.2所謂數(shù)值解法，就是尋求解在一系列離散節(jié)點(diǎn)上的近似值.相鄰兩個(gè)節(jié)點(diǎn)的間距稱為步長(zhǎng).如不特別說(shuō)明，總是假定為定數(shù)，這時(shí)節(jié)點(diǎn)為.初值問(wèn)題（1.1），（1.2）的數(shù)值解法的基本特點(diǎn)是采取“步進(jìn)式”.即求解過(guò)程順著節(jié)點(diǎn)排列的次序一步一步地向前推進(jìn).（1.1）（1.2）3描述這類算法，只要給出用已知信息計(jì)算的遞推公式.一類是計(jì)算時(shí)只用到前一點(diǎn)的值，稱為單步法.另一類是用到前面點(diǎn)的值，稱為步法.其次，要研究公式的局部截?cái)嗾`差和階，數(shù)值解與精確解的誤差估計(jì)及收斂性，還有遞推公式的計(jì)算穩(wěn)定性等問(wèn)題.首先對(duì)方程離散化，建立求數(shù)值解的遞推公式.45.2

簡(jiǎn)單的數(shù)值方法與基本概念

5.2.1

Euler

Euler方法

積分曲線上一點(diǎn)的切線斜率等于函數(shù)的值.如果按函數(shù)在平面上建立一個(gè)方向場(chǎng)，那么，積分曲線上每一點(diǎn)的切線方向均與方向場(chǎng)在該點(diǎn)的方向一致.在平面上，微分方程的解稱作它的積分曲線.5基于上述幾何解釋，從初始點(diǎn)出發(fā)，先依方向場(chǎng)在該點(diǎn)的方向推進(jìn)到上一點(diǎn)，然后再?gòu)囊婪较驁?chǎng)的方向推進(jìn)到上一點(diǎn)，循此前進(jìn)做出一條折線（圖5-1).圖5-16一般地，設(shè)已作出該折線的頂點(diǎn)，過(guò)依方向場(chǎng)的方向再推進(jìn)到，顯然兩個(gè)頂點(diǎn)的坐標(biāo)有關(guān)系即（2.1）這就是著名的歐拉（Euler）公式.若初值已知，則依公式（2.1）可逐步算出7

例1（2.2）求解初值問(wèn)題

解取步長(zhǎng)，歐拉公式的具體形式為計(jì)算結(jié)果見(jiàn)表5-1.初值問(wèn)題（2.2）的解為，按這個(gè)解析式子算出的準(zhǔn)確值同近似值一起列在表5-1中，兩者相比較可以看出Euler方法的精度很差.8還可以通過(guò)幾何直觀來(lái)考察Euler方法的精度.假設(shè)，即頂點(diǎn)落在積分曲線上，那么，按Euler方法作出的折線便是過(guò)點(diǎn)的切線（圖5-2）.9圖5-2從圖形上看，這樣定出的頂點(diǎn)明顯地偏離了原來(lái)的積分曲線，可見(jiàn)Euler方法是相當(dāng)粗糙的.為了分析計(jì)算公式的精度，通?？捎肨aylor展開(kāi)將在處展開(kāi)，則有10在的前提下，稱為此方法的局部截?cái)嗾`差.于是可得Eulor法（2.1）的公式誤差（2.3）（2.4）（2.1）如果對(duì)方程從到積分，得11右端積分用左矩形公式近似.再以代替如果在（2.4）中右端積分用右矩形公式（2.5）稱為后退的Euler法.

Euler公式是關(guān)于的一個(gè)直接的計(jì)算公式，這類公式代替也得到（2.1），局部截?cái)嗾`差也是（2.3）.近似，則得另一個(gè)公式稱作是顯式的；（2.1）（2.3）（2.4）12公式（2.5）的右端含有未知的，它是關(guān)于的一個(gè)函數(shù)方程，隱式方程通常用迭代法求解，而迭代過(guò)程的實(shí)質(zhì)是逐步顯示化.設(shè)用歐拉公式給出迭代初值，用它代入（2.5）式的右端，使之轉(zhuǎn)化為顯式，直接計(jì)算得這類公式稱作是隱式的.（2.5）（2.5）13然后再用代入（2.5）式，又有如此反復(fù)進(jìn)行，得（2.6）由于對(duì)滿足利普希茨條件（1.3）.由（2.6）減（2.5）得由此可知，只要迭代法（2.6）就收斂到解.（2.5）（2.5）（1.3）14從積分公式可以看到后退Euler方法的公式誤差與Euler法是相似的.15

5.2.2

梯形方法

若用梯形求積公式近似等式（2.4）右端的積分并分別用代替則可得到比Euler法精度高的計(jì)算公式（2.7）稱為梯形方法.梯形方法是隱式單步法，可用迭代法求解.（2.4）16為了分析迭代過(guò)程的收斂性，將（2.7）與（2.8）式相減，得（2.8）同后退的Euler方法一樣，仍用Euler方法提供迭代初值，則梯形法的迭代公式為（2.7）17如果選取充分小，使得則當(dāng)時(shí)有這說(shuō)明迭代過(guò)程（2.8）是收斂的.于是有式中為關(guān)于的利普希茨常數(shù).18

5.2.3

單步法的局部截?cái)嗾`差與階

初值問(wèn)題(1.1)，(1.2)的單步法可用一般形式表示為（2.9）其中多元函數(shù)與有關(guān)，當(dāng)含有時(shí)，方法是隱式的，若不含則為顯式方法，（2.10）稱為增量函數(shù).所以顯式單步法可表示為例如對(duì)歐拉法（2.1）有它的局部截?cái)嗾`差已由（2.3）給出.（1.1）（1.2）（2.1）（2.3）19對(duì)一般顯式單步法則可如下定義.

定義1設(shè)是初值問(wèn)題(1.1)，(1.2)的準(zhǔn)確解，（2.11）為顯式單步法（2.10）的局部截?cái)嗾`差.

之所以稱為局部的，是假設(shè)在前各步?jīng)]有誤差.當(dāng)時(shí)，計(jì)算一步，則有稱（1.1）（1.2）（2.10）20在前一步精確的情況下用公式（2.10）計(jì)算產(chǎn)生的公式誤差.根據(jù)定義，Euler法的局部截?cái)嗾`差即為（2.3）的結(jié)果.這里稱為局部截?cái)嗾`差主項(xiàng).局部截?cái)嗾`差可理解為用方法（2.10）計(jì)算一步的誤差，即顯然（2.10）（2.10）（2.3）21

定義2設(shè)是初值問(wèn)題(1.1)，(1.2)的準(zhǔn)確解，若存在最大整數(shù)使顯式單步法(2.10)的局部截?cái)嗾`差滿足（2.12）則稱方法（2.10）具有階精度.若將（2.12）展開(kāi)式寫(xiě)成則稱為局部截?cái)嗾`差主項(xiàng).以上定義對(duì)隱式單步法（2.9）也是適用的.（1.1）（1.2）（2.10）（2.10）（2.9）22對(duì)后退Euler法（2.5）其局部截?cái)嗾`差為這里，是1階方法，局部截?cái)嗾`差主項(xiàng)為.（2.5）23對(duì)梯形法（2.7）有所以梯形方法是2階的，其局部誤差主項(xiàng)為（2.7）24

5.2.4

改進(jìn)的歐拉公式

梯形方法雖然提高了精度，但其算法復(fù)雜.在應(yīng)用迭代公式（2.9）進(jìn)行實(shí)際計(jì)算時(shí)，每迭代一次，都要重新計(jì)算函數(shù)的值.為了控制計(jì)算量，通常只迭代一兩次就轉(zhuǎn)入下一步的計(jì)算，這就簡(jiǎn)化了算法.具體地，先用Euler公式求得一個(gè)初步的近似值,而迭代又要反復(fù)進(jìn)行若干次，計(jì)算量很大，而且往往難以預(yù)測(cè).稱之為預(yù)測(cè)值，（2.9）25這樣建立的預(yù)測(cè)-校正系統(tǒng)通常稱為改進(jìn)的Euler公式：預(yù)測(cè)值的精度可能很差，再用梯形公式（2.7）將它校正一次，即按（2.8）式迭代一次得，這個(gè)結(jié)果稱校正值.預(yù)測(cè)校正（2.13）也可以表為下列平均化形式（2.7）（2.8）26

例2

解用改進(jìn)的Euler方法求解初值問(wèn)題（2.2）.（2.2）這里改進(jìn)的Euler公式為27仍取，計(jì)算結(jié)果見(jiàn)表5-2.同例1中Euler法的計(jì)算結(jié)果比較，改進(jìn)Euler法明顯改善了精度.285.3Runge-Kutta方法29

5.3.1

Runge-Kutta法的一般形式上節(jié)給出了顯式單步法的表達(dá)式其局部截?cái)嗾`差為對(duì)歐拉法，即方法為階.（3.1）若用改進(jìn)Euler法，它可表示為30此時(shí)增量函數(shù)（3.2）與歐拉法的相比，增加了計(jì)算一個(gè)右函數(shù)的值，可望.若要使得到的公式階數(shù)更大，就必須包含更多的值.（3.3）從方程等價(jià)的積分形式（2.4），即31若要使公式階數(shù)提高，就必須使右端積分的數(shù)值求積公式精度提高，必然要增加求積節(jié)點(diǎn).為此可將（3.3）的右端用求積公式表示為點(diǎn)數(shù)越多，精度越高，上式右端相當(dāng)于增量函數(shù)，為得到便于計(jì)算的顯式方法，可類似于改進(jìn)Euler法，將公式表示為（3.4）其中（3.3）32（3.5）這里均為常數(shù).

(3.4)和(3.5)稱為級(jí)顯式龍格-庫(kù)塔（Runge-Kutta）法，簡(jiǎn)稱R-K方法.當(dāng)時(shí)，就是Euler法，此時(shí)方法的階為.當(dāng)時(shí)，改進(jìn)Euler法(3.1)，(3.2)也是其中的一種.（3.4）33下面將證明階.要使公式(3.4),(3.5)具有更高的階,就要增加點(diǎn)數(shù).下面就推導(dǎo)R-K方法.（3.4）（3.5）34

5.3.2

二階R-K方法

對(duì)的R-K方法，計(jì)算公式如下（3.6）這里均為待定常數(shù).希望適當(dāng)選取這些系數(shù)，使公式階數(shù)盡量高.根據(jù)局部截?cái)嗾`差的定義，（3.6）的局部截?cái)嗾`差為（3.7）35這里.為得到的階，要將上式各項(xiàng)在處做Talor展開(kāi)，由于是二元函數(shù)，故要用到二元Talor展開(kāi)，其中各項(xiàng)展開(kāi)式為（3.8）36將以上結(jié)果代入局部截?cái)嗾`差公式則有要使公式（3.6）具有階，必須使（3.6）37即（3.9）的解是不惟一的.令，則得這樣得到的公式稱為二階R-K方法，如取，則這就是改進(jìn)Euler法（3.1）.（3.9）38若取則.稱為中點(diǎn)公式，（3.10）也可表示為得計(jì)算公式（3.10）相當(dāng)于數(shù)值積分的中矩形公式.的R-K公式(3.6)的局部誤差不可能提高到.（3.6）39把多展開(kāi)一項(xiàng)，從（3.8）的看到展開(kāi)式中的項(xiàng)是不能通過(guò)選擇參數(shù)消掉的.實(shí)際上要使的項(xiàng)為零，需增加3個(gè)方程，要確定4個(gè)參數(shù)，這是不可能的.故的顯式R-K方法的階只能是，而不能得到三階公式.40

5.3.3

三階與四階R-K方法

要得到三階R-K方法，必須.（3.11）其中及均為待定參數(shù).此時(shí)（3.4）,（3.5）的公式表示為公式（3.11）的局部截?cái)嗾`差為（3.4）（3.5）41只要將按二元函數(shù)Taylor展開(kāi)，使，可得待定參數(shù)滿足方程（3.12）42這是8個(gè)未知數(shù)6個(gè)方程的方程組，解也不是惟一的.所以這是一簇公式.滿足條件（3.12）的公式（3.11）統(tǒng)稱為三階R-K公式.一個(gè)常見(jiàn)的公式為此公式稱為庫(kù)塔三階方法.43繼續(xù)上述過(guò)程，經(jīng)過(guò)較復(fù)雜的數(shù)學(xué)演算，可以導(dǎo)出各種四階R-K公式，下列經(jīng)典公式是其中常用的一個(gè)：可以證明其截?cái)嗾`差為.四階Runge-Kutta方法的每一步需要計(jì)算4次函數(shù)值，（3.13）44

例3

解設(shè)取步長(zhǎng)，從到用四階龍格-庫(kù)塔方法求解初值問(wèn)題這里，經(jīng)典的四階龍格-庫(kù)塔公式為45表5-3列出了計(jì)算結(jié)果，同時(shí)列出了相應(yīng)的精確解.比較例3和例2的計(jì)算結(jié)果，顯然龍格-庫(kù)塔方法的精度高.46但由于這里放大了步長(zhǎng)，所以表5-3和表5-2所耗費(fèi)的計(jì)算量幾乎相同.龍格-庫(kù)塔方法的推導(dǎo)基于泰勒展開(kāi)方法，因而它要求所求的解具有較好的光滑性質(zhì).反之，如果解的光滑性差，那么，使用四階龍格-庫(kù)塔方法求得的數(shù)值解，其精度可能反而不如改進(jìn)的歐拉方法.四階龍格-庫(kù)塔方法每一步要4次計(jì)算函數(shù),其計(jì)算量比每一步只要2次計(jì)算函數(shù)的二階龍格-庫(kù)塔方法----改進(jìn)的歐拉方法大一倍.47

5.3.4

變步長(zhǎng)的龍格-庫(kù)塔方法