計算方法 第5章 數(shù)值積分與數(shù)值微分

第五章 數(shù)值積分與數(shù)值微分在高等數(shù)學中我們學過定積分的計算方法，若找到被積函數(shù)在區(qū)間上的一個原函數(shù)，利用Newton-Leibniz公式可以輕易得計算出積分值，但在實際問題中，往往會遇到一些困難。有些函數(shù)雖然能找到原函數(shù),但表達式過于復(fù)雜，例如的原函數(shù)為有些函數(shù)找不到初等函數(shù)形式的原函數(shù)，例如積分有些情況下，函數(shù)值是用表格形式給出的，例如：對于以上這些積分問題，解決的方法就是使用數(shù)值積分方法。其實數(shù)值積分方法不僅可以解決上述問題，最為重要的優(yōu)點是對任意被積函數(shù)任意積分區(qū)間的積分問題都可以采用統(tǒng)一的數(shù)值積分公式，非常便于計算機編程實現(xiàn)。 對于微分問題，雖然對每一個初等函數(shù)都可以求出其導(dǎo)數(shù)，但是不同函數(shù)其求導(dǎo)方法依賴于各自不同的求導(dǎo)公式，沒有簡單、統(tǒng)一的處理方法，而數(shù)值微分法卻可以對不同的函數(shù)使用統(tǒng)一的數(shù)值微分公式或數(shù)值微分算法。 本章首先介紹一些數(shù)值積分公式，最后再簡單的介紹數(shù)值微分問題。5.1 數(shù)值積分公式 1.數(shù)值積分的基本思想 我們知道定積分的幾何意義就是所圍成的曲邊形面積，而數(shù)值積分的基本思想是利用函數(shù)在區(qū)間上某些點處函數(shù)值的線性組合來計算其定積分的近似值，把計算定積分這一復(fù)雜問題轉(zhuǎn)換為僅僅涉及到函數(shù)值的計算問題，而無需考慮函數(shù)本身的結(jié)構(gòu)以及函數(shù)值的真實來源，這樣就很便于計算機編程實現(xiàn)。 為求的數(shù)值解，我們考慮在積分區(qū)間上用的插值函數(shù)近似，則有 (5.1)稱(5.1)式為函數(shù)在區(qū)間上的求積公式，為求積系數(shù)，插值函數(shù)的插值節(jié)點在這里稱為求積節(jié)點，如果記 (5.2)稱(5.2)式為求積公式(5.1)的截斷誤差，實質(zhì)上求積公式的截斷誤差就是插值余項所致，求積公式的精確程度可由截斷誤差的絕對值來度量。從(5.1)式還容易看出求積系數(shù)僅僅與節(jié)點的選取有關(guān)，而不依賴被積函數(shù)，因此求積公式適用于所有函數(shù)，當然其截斷誤差的大小與節(jié)點的選擇和函數(shù)本身的特性均有關(guān)系。 2.求積公式的代數(shù)精確度 定義5.1如果求積公式，對于是精確成立的，而對于不再精確成立，則稱該求積公式的代數(shù)精確度為。根據(jù)定義5.1，如果求積公式的代數(shù)精確度為，則求積系數(shù)應(yīng)滿足 ……………. 如果將求積公式中的系數(shù)作為未知量，則我們已經(jīng)列出了確定求積系數(shù)的線性方程組，選出節(jié)點，取，求解此線性方程組，即可以通過求解此方程組確定出求積系數(shù)，得到其求積公式且求積公式的代數(shù)精確度至少為。 例5.1試確定中的待定參數(shù)，使該積分公式的代數(shù)精確度盡可能的高，并指出所求積分公式的代數(shù)精確度。 解：1)求待定系數(shù)，利用求積公式，得到關(guān)于的方程組 解之得 2)確定求積公式的代數(shù)精確度，按此方程組確定的待定系數(shù)，求積公式的代數(shù)精確度至少為2，現(xiàn)在依次驗證時的求積公式是否精確成立，直到出現(xiàn)首個不成立為止。 則直到均精確成立，但對于不再精確成立，表明所求的求積公式的代數(shù)精確度為5。用解線性方程組的方法構(gòu)造求積公式很不方便，通常是利用以為插值節(jié)點的插值多項式來構(gòu)造出常用的求積公式，在編程中直接使用這些常用求積公式，以下我們將構(gòu)造這些求積公式。3．梯形求積公式對于，取，用近似，則有稱積分公式 (5.3)為梯形求積公式。梯形求積公式具有直觀的幾何意義：圖5.1梯形求積公式的幾何意義關(guān)于梯形求積公式的截斷誤差，有如下結(jié)論。定理5.1若函數(shù)在上有連續(xù)的二階導(dǎo)數(shù)，則梯形求積公式的截斷誤差為 (5.4) 證明：由于梯形求積公式的誤差是由插值余項而引起的，其余項公式為其中與有關(guān)，再注意到在上有連續(xù)的二階導(dǎo)數(shù)，記，則在區(qū)間內(nèi)且不變符號，由廣義積分中值定理有 ，對上式右端求積分可得例5.2驗證梯形求積公式的代數(shù)精確度為1。 解： 對于精確成立，但對于不成立，所以代數(shù)精度為1，從梯形公式的代數(shù)精確度太小也可以看出梯形求積公式的精度不高。 4.Simpson求積公式 對于，將區(qū)間劃分為二等分,令，取，用近似，則有 稱積分公式 (5.5)為Simpson求積公式，Simpson求積公式的幾何意義：圖5.2Simpson求積公式的幾何意義關(guān)于Simpson求積公式的截斷誤差，仿照前邊定理5.1的證明過程可證明如下定理。定理5.2若函數(shù)在區(qū)間上有連續(xù)的4階導(dǎo)數(shù)，則Simpson求積公式的截斷誤差為 (5.6)例5.3求Simpson求積公式的代數(shù)精確度。解：對于，考察 所以Simpson求積公式的代數(shù)精確度為3，從代數(shù)精確度、截斷誤差和實際計算效果(稍后的例5.4比較了3個求積公式的計算結(jié)果)看，Simpson求積公式都優(yōu)于梯形求積公式。 5.Cotes求積公式 對于，將區(qū)間劃分為4等分，令，取，用近似，則有稱積分公式 (5.7)為Cotes求積公式，其截斷誤差為 (5.8)可以驗證Cotes求積公式的代數(shù)精度為5。 6．Newton-Cotes求積公式 現(xiàn)在我們討論一般情況下的求積公式，先將積分區(qū)間劃分為等分，令，積分節(jié)點為，用近似，則有 (5.9)稱此公式為Newton-Cotes求積公式，其中 (5.10)令，則對于有稱為Cotes系數(shù)?？梢宰C明當n為偶數(shù)時，Newton-Cotes型求積公式的代數(shù)精確度為n+1，當n為奇數(shù)時，其代數(shù)精確度為n。Cotes系數(shù)表從Cotes系數(shù)表可以看出，梯形求積公式、Simpson求積公式均是Newton-Cotes積分公式的特例。另一方面從Cotes系數(shù)的表達式還可以看出，Cotes系數(shù)與積分區(qū)間無關(guān)，與被積函數(shù)也無關(guān)，僅與區(qū)間的等份個數(shù)有關(guān)。Cotes系數(shù)具有如下性質(zhì)。由于Cotes系數(shù)是由插值公式推導(dǎo)而來，而插值公式在高階情況下具有不穩(wěn)定性，因此Cotes公式必然也是不穩(wěn)定的，從Cotes系數(shù)表中的出現(xiàn)負值就可以看出問題，通常情況下，Newton-Cotes求積公式中最常用的是梯形求積公式、Simpson求積公式和4階Cotes求積公式，通常不用高階的Newton-Cotes求積公式。 例5.4分別用梯形求積公式、Simpson求積公式和4階Cotes求積公式計算定積分的近似值。 解：梯形求積公式： Simpson求積公式： Cotes求積公式： 和精確值比較，梯形求積公式、Simpson求積公式和4階Cotes求積公式所求結(jié)果分別有2、4、6位有效數(shù)字。5.2 復(fù)化求積公式 Newton-Cotes求積公式的兩個問題高階求積公式穩(wěn)定性難以保證；從截斷誤差公式可看出，低階求積公式對于積分區(qū)間較大時同樣精度也難以保證。解決方法： 針對前述的兩個問題，注意到前節(jié)中各截斷誤差公式中均有因子，如果將區(qū)間縮小定會提高積分的精度，復(fù)化求積公式就是將大的積分區(qū)間分解為小的積分區(qū)間，在小的積分區(qū)間上使用低階求積公式。復(fù)化梯形公式將積分區(qū)間劃分為等分，則每個子區(qū)間長度均為，稱為步長，此時分點為，在每個子區(qū)間上施行梯形公式，得到稱 (5.11)為復(fù)化梯形公式。由定積分的定義，顯然有此極限也說明，當N充分大時，復(fù)化梯形公式會充分接近積分的精確值。復(fù)化梯形公式的截斷誤差由以下定理給出。定理5.3若函數(shù)在上有連續(xù)的二階導(dǎo)數(shù)，則復(fù)化梯形公式截斷誤差為 (5.12)其截斷誤差的上限為 (5.13) 證明：復(fù)化梯形公式的截斷誤差就是每個小區(qū)間上的截斷誤差之和，根據(jù)(5.4)式有，由于在上有連續(xù)的二階導(dǎo)數(shù)，根據(jù)高等數(shù)學中的介值定理，存在使從而有，截斷誤差限估計式(5.13)式是顯然成立的。復(fù)化Simpson公式類似于復(fù)化梯形公式的構(gòu)造方法，再構(gòu)造復(fù)化Simpson公式，將積分區(qū)間劃分為等分，則每個子區(qū)間長度均為，在每個子區(qū)間上施行Simpson公式(此時有個分點)，得到稱 (5.14)為復(fù)化Simpson公式，對于復(fù)化Simpson公式同樣也有這也表明當N充分大時，復(fù)化Simpson公式也能很好的近似定積分的精確值。關(guān)于復(fù)化Simpson公式截斷誤差，有比復(fù)化梯形公式更好的結(jié)論。 定理5.4若函數(shù)在上有連續(xù)的4階導(dǎo)數(shù)，則復(fù)化Simpson公式截斷誤差為 (5.15)其截斷誤差上限 (5.16)復(fù)化Cotes公式類似于復(fù)化梯形公式和復(fù)化Simpson公式的構(gòu)造過程可以得到復(fù)化Cotes公式 (5.17)其中，稱該式為復(fù)化Cotes公式，容易看出復(fù)化Cotes公式共有4N+1個分點。同樣根據(jù)定積分的定義也有 定理5.5若函數(shù)在上有連續(xù)的6階導(dǎo)數(shù)，則復(fù)化Cotes公式截斷誤差為 (5.18)且截斷誤差的上限為 (5.19) 例5.5試取9個節(jié)點分別用復(fù)化梯形公式、復(fù)化Simpson公式、復(fù)化Cotes公式計算定積分。 解：首先建立如下表格根據(jù)此表可以計算出的結(jié)果如下： 所求積分的精確值為，則的誤差分別為，，，則有效數(shù)字分別為2、7、6位。此例也說明同樣計算9次函數(shù)值，計算量相同，使用不同的公式所得結(jié)果的精度相差很多。 在復(fù)化積分公式中，編者推薦使用復(fù)化Simpson公式，此公式簡單易記且精度高。 例5.6如果計算積分的精度要求為，試問分別用復(fù)化梯形公式、復(fù)化Simpson公式、復(fù)化Cotes公式，積分區(qū)間至少應(yīng)劃分多少等分？ 解：，則 復(fù)化梯形公式： 復(fù)化Simpson公式： 復(fù)化Cotes公式： 因此可以得到如下結(jié)論，對于所給的定積分，到達指定的精度，復(fù)化梯形公式需要674等分積分區(qū)間，也就是將要計算675個分點的函數(shù)值；復(fù)化Simpson公式需要7等分積分區(qū)間，總共需要計算個分點的函數(shù)值；復(fù)化Cotes公式需要2等分積分區(qū)間，總共需要計算個分點的函數(shù)值。5.3 Romberg求積法復(fù)化公式的一個缺點是必須先給出步長，從截斷誤差公式可以看出，復(fù)化求積公式的計算精度嚴重依賴于步長，而如何選取步長也沒有給出一個可行的方法。本節(jié)所給出的求積公式采用變步長的方法，把步長逐次分半，直至達到要求的精度為止，這就是本節(jié)所要討論的內(nèi)容—Romberg求積法。變步長梯形公式變步長梯形公式的思想是，首先取步長是整個積分區(qū)間，使用梯形公式計算積分值，逐次將步長減半，使用復(fù)化梯形公式計算積分值，當積分值的計算精度到達要求時停止計算，否則為步長繼續(xù)減半。在算法實施過程中，我們要解決三個問題。1)在計算出后，如何計算，即導(dǎo)出與的遞推公式。2)如何避免重復(fù)計算相同節(jié)點函數(shù)值？3)在計算出后，如何估計其誤差，即算法的終止準則是什么？首先推導(dǎo)變步長梯形公式的遞推公式，在計算時，需要計算個點上的函數(shù)值。在求得后再計算時，需將每個子區(qū)間再二等分，新增加個分點。為了避免重復(fù)計算，計算時，已經(jīng)計算出的個點上的函數(shù)值應(yīng)保留下來，而只計算新增加的個分點上的函數(shù)值，為此把表示成兩部分之和。即這樣就可以得到積分值遞推公式 (5.20)有了遞推公式同時也就解決了重復(fù)計算問題。以下我們根據(jù)和的截斷誤差，推導(dǎo)出算法的終止準則，按照(5.12)式有這里，如果在積分區(qū)間內(nèi)的變化不大，可以假設(shè)，綜合以上兩式則有該式表明，用作為定積分的近似值，其誤差大約為，因此終止準則可以用和的差來度量定積分計算值的精度。即對于給定的精度，算法終止準則為 (5.21)例5.7用變步長梯形公式求定積分，要求終止準則為。 解：首先用梯形公式得到，再利用遞推公式得到、、、…直到滿足精度要求為止，計算結(jié)果如下，，，，精確值為，。Romberg(龍貝格)求積公式我們先利用前邊推導(dǎo)終止準則時得到的近似公式 (5.22)對例5.7中的和數(shù)據(jù)進行計算該值與精確值的誤差為，達到的四位有效數(shù)字，實際上好于的精度。 基于這一事實，我們猜想把(5.22)這一近似公式作為新的積分近似值應(yīng)該優(yōu)于，正是這一猜想產(chǎn)生了非常有效的Romberg積分法。首先我們可以驗證 (5.23)仿照對梯形公式的推導(dǎo)過程，對做同樣的處理，則有如果在積分區(qū)間內(nèi)的變化不大，可以假設(shè)，則有通過驗算還可以得出 (5.24)繼續(xù)對做同樣處理可推出還可以得到 (5.25)可以想象要優(yōu)于，用計算積分其精確將可能更高，如此遞推所得的遞推公式稱為Romberg遞推公式。Romberg遞推公式的算法設(shè)計通常將該計算步驟設(shè)計成以下表格來遞推計算積分值便于編程的Romberg遞推公式為 ，注意在此遞推格式中內(nèi)循環(huán)的是同步變化的，直到完成為止。在編程時可以將遞推公式中的作以下等價變換其中后一項就是本次遞推的修正項，也可以作為終止精度的判斷項。如同第4章計算差商表一樣，最新一行近似積分值的存放方式按倒序更方便，即每次迭代對角線元素為首個元素，算法結(jié)束時首元素就是所求積分值，設(shè)計算法如下：算法5.1Romberg求積法[初始化]算法中將最后一行按反序存入向量，令初始步長，給定終止精度，，計算首個梯形值；令，；[更新梯形值]；[遞推]令，然后對循環(huán)3.1)3.2)[判斷終止準則]如果，算法結(jié)束，為滿足精度要求的積分近似值。，轉(zhuǎn)1)。[算法結(jié)束]例5.8用Romberg遞推公式求定積分，終止精度為。 解：同樣先用梯形公式得到，再利用Romberg遞推公式得到以下結(jié)果20.78333330.03333330.002064840.78552940.00219610.000131280.78539640.00013300.0000017建議：對于已知函數(shù)表達式的情況計算定積分推薦使用Romberg算法，對于給出的僅僅是離散函數(shù)值的情況推薦使用復(fù)化Simpson公式。Richardson(理查森)外推加速法前面根據(jù)猜想的方式得到了Romberg求積公式并驗證了有效性，現(xiàn)在我們探討其理論依據(jù)，Romberg求積法本質(zhì)上是用誤差補償?shù)募铀偎惴?，由收斂較慢的序列構(gòu)造成收斂很快的序列，此方法即是Richardson外推加速法。記是以為步長的梯形值，可以證明當被積函數(shù)有任意階導(dǎo)數(shù)時，有 (5.26)其中為積分的理論值，與無關(guān)，該式說明，在被積函數(shù)和積分區(qū)間確定的情況下，Romberg遞推公式中的每一步積分值可以看作是步長的函數(shù)。 從(5.26)式可以看出，，即用復(fù)化梯形公式計算積分，其計算誤差階為。如果在式中用來代替，則得到 (5.27)用4乘(5.27)式減去(5.26)式再除則消除了項，記為，則得 (5.28)顯然對理論積分值的精度達到。前邊我們已經(jīng)驗證過是以為步長的復(fù)化Simpson公式，如果在中再用來代替又可以得到此時我們再消除其中的項，用乘，減去(5.28)式再除將其記為，則得到下式顯然對理論積分值的精度達到，就是步長為的復(fù)化Cotes公式。依次類推就得到Romberg求積法。我們從遞推過程可以看出，每次遞推其步長將折半，對積分理論值的精度又提高兩階，這就是Romberg求積法高效的原因所在。其實Richardson外推加速法不僅可以對梯形積分值序列加速，只要是滿足以上特征的所有序列均可加速，在下節(jié)數(shù)值微分中還會用到Richardson外推加速策略，現(xiàn)在我們看另外一個例子，體會一下Richardson外推加速的用途。例5.9已知有如下恒等式試依據(jù)此式的值按，用Richardson外推法求的近似值。 解：設(shè)，則有使用Richardson外推法外推一次得到再外推一次可得到將計算結(jié)果列于如下表中精確值為，誤差為。5.4 數(shù)值微分 所謂數(shù)值微分就是僅利用函數(shù)在一些點處的函數(shù)值求出其導(dǎo)數(shù)的近似值，其優(yōu)點是可以對不同的函數(shù)使用同一算法求各自對應(yīng)的導(dǎo)數(shù)，便于計算機編程實現(xiàn)。本節(jié)將討論數(shù)值微分的差商法、Richardson外推法、插值法等。差商代替導(dǎo)數(shù)為計算的導(dǎo)數(shù)，取一較小的步長，考慮的一階Taylor展開從該式可以解出 (5.29)同樣如果考慮的一階Taylor展開則有 (5.30)如果再考慮和的Taylor展開兩式相減則有 (5.31)以上(5.29)、(5.30)、(5.31)三個公式分別稱為向前差商、向后差商、中心差商公式，從無窮小項的階數(shù)可以看出中心差商近似導(dǎo)數(shù)更為精確。如果計算是精確的，顯然愈小差商愈能逼近導(dǎo)數(shù)值，但從誤差分析的角度討論，太小將意味著兩個相近的數(shù)相減，勢必導(dǎo)致有效數(shù)字的損失，這一點我們在第一章的例1.2已經(jīng)看到過，現(xiàn)在再看一個例子。例5.10用中心差商計算在處的導(dǎo)數(shù)。解：使用中心差商公式(5.31)，有取以下不同的步長，按5位有效數(shù)字進行計算，結(jié)果見如下表格(精確值0.3535533906)按5位有效數(shù)字進行計算在時誤差最小，達3位有效數(shù)字，誤差也不算大，反而是時僅1位有效數(shù)字。讀者可以自己編程試算一下，如果用雙精度進行計算，取，計算將達到8位有效數(shù)字，當后誤差將逐步增大；如果用單精度進行計算，時精度最高計算達到6位有效數(shù)字，當后將無一位有效數(shù)字。從誤差分析的角度，本例的計算應(yīng)該按如下方式用計算已達4為有效數(shù)字，用5位有效數(shù)字計算，達到4位有效數(shù)字已算相當精確。Richardson外推法計算微分使用差商近似導(dǎo)數(shù)并不是一個讓人滿意的選擇，但是對中心差商進行外推卻可以得到很好的算法，我們考察中心差商公式：根據(jù)Taylor展開公式，有兩式相減并除以h得這里常數(shù)與無關(guān)，顯然符合Richardson外推形式，從而有如下外推公式 (5.32)外推格式表如下所示注意到對的計算，當太大也就是步長太小時將導(dǎo)致兩個非常相近的數(shù)相減，造成有效數(shù)字的損失，因此使用外推要適可而止，一般控制步長，從以下例5.11也可看出外推法中的步長不須太小就可以得到較為精確的近似解。例5.11運用外推法計算在處的導(dǎo)數(shù)。 解：計算結(jié)果見下表(計算保留5位數(shù)字)精確值為e=2.718281828…。我們可以仿照外推法計算積分的算法5.1設(shè)計如下算法，算法中將外推表的最后一行按反序存入向量，因此算法結(jié)束時就是導(dǎo)數(shù)值。 算法5.2外推法計算導(dǎo)數(shù)[初始化]取初始步長，給定終止精度，，計算首個導(dǎo)數(shù)值；令，；[更新導(dǎo)數(shù)]；[遞推]令，然后對循環(huán)3.1)3.2)[判斷終止準則]如果，算法結(jié)束，為滿足精度要求的近似導(dǎo)數(shù)值。，轉(zhuǎn)1)。[算法結(jié)束]插值函數(shù)計算微分利用插值函數(shù)計算導(dǎo)數(shù)也是數(shù)值微分的一種方法。其思想很簡單，就是對于給定的的函數(shù)值表，建立的插值函數(shù)，令，在某些情況下，用這種方法計算導(dǎo)數(shù)是可以接受的選擇?？紤]函數(shù)在插值區(qū)間上有個插值節(jié)點，其插值函數(shù)為，則對此等式兩邊求導(dǎo)令則有其中介于之間，顯然用插值函數(shù)的導(dǎo)數(shù)近似表示在處的導(dǎo)數(shù)其誤差由上式中的余項決