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

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

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

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

5.2.1

Euler

Euler方法

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

例1（2.2）求解初值問題

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

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

5.2.2

梯形方法

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

5.2.3

單步法的局部截斷誤差與階

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

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

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

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

5.2.4

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

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

例2

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

5.3.1

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

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

5.3.2

二階R-K方法

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

5.3.3

三階與四階R-K方法

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

例3

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

5.3.4

變步長的龍格-庫塔方法