Numerical Differentiation

1、Chapter 6 Numerical Differentiation數(shù)值微分?jǐn)?shù)值導(dǎo)數(shù)的公式對于開發(fā)求解常微分方程（ODE）和偏微分方程（PDE）邊值問題的算法很重要（參見第九章和第十章）。數(shù)值微分中經(jīng)常采用已知函數(shù)為例，以便數(shù)值逼近的結(jié)果與準(zhǔn)確值相比較。為了說明問題，考慮貝塞爾函數(shù)（the Bessel function），它的值列表可以在標(biāo)準(zhǔn)參考書中查到。閉區(qū)間上八個等距點(diǎn)和。它所基于的原理是插值多項(xiàng)式的微分?？紤]尋求。插值多項(xiàng)式經(jīng)過三個結(jié)點(diǎn)和，并由此得到。這個二次多項(xiàng)式和它在點(diǎn)處的切線如圖所示。如果使用五個插值點(diǎn)，可以得到更好的逼近。經(jīng)過點(diǎn)和的多項(xiàng)式，并由此得到。這個四次多項(xiàng)式和它在點(diǎn)處

2、的切線如圖所示。這個導(dǎo)數(shù)的準(zhǔn)確值是，在使用和時的誤差分別為和。在這一章，我們研究數(shù)值微分精確度的相關(guān)理論。此處所舉例子是可查值列表的貝塞爾函數(shù)，以便進(jìn)行逼近值與真值的比較。采用的方法是用與函數(shù)相應(yīng)的插值多項(xiàng)式的導(dǎo)數(shù)值來近似其導(dǎo)數(shù)值。分別采用了二次插值多項(xiàng)式和四次插值多項(xiàng)式來近似，并計算了相應(yīng)誤差，顯然后者近似效果優(yōu)于前者。事實(shí)上，這里提供了一種近似函數(shù)導(dǎo)數(shù)值的方法，即用與函數(shù)相應(yīng)的插值多項(xiàng)式的導(dǎo)數(shù)值來近似其導(dǎo)數(shù)值。還要進(jìn)行的思考是，這種方法的近似效果如何？怎么提高近似精度？再就是，有沒有其他的近似方法？在前面的插值理論中已經(jīng)討論過，若函數(shù)的導(dǎo)數(shù)在定義閉區(qū)間上一致有界，且步長，則選擇較大的將得到

3、較小的（即有誤差界），從而高次逼近多項(xiàng)式將產(chǎn)生較小的誤差。也討論過高次逼近時在區(qū)間上發(fā)生的振蕩（節(jié)點(diǎn)數(shù)增加，振蕩更劇烈），即龍格現(xiàn)象，并提出了應(yīng)對策略分段低次插值。所以說，通過插值多項(xiàng)式的導(dǎo)數(shù)來近似函數(shù)的導(dǎo)數(shù)值是切實(shí)可行的，并且由于多項(xiàng)式求導(dǎo)簡單易行，這種方法就顯得更加方便?？紤]近似精度的話，相對而言，用高次插值多項(xiàng)式近似誤差要比低次的小，當(dāng)然這是在不致引起振蕩（次數(shù)不太高）的前提下。畢竟插值多項(xiàng)式的導(dǎo)出有時是很復(fù)雜的，比如三次樣條插值函數(shù)，需要數(shù)次求導(dǎo)積分并利用函數(shù)值、導(dǎo)數(shù)連續(xù)以及邊界條件求待定系數(shù)來推倒得到，計算復(fù)雜，況且我們所需的只是逼近導(dǎo)數(shù)值，沒有必要把插值多項(xiàng)式（函數(shù)的逼近）求出來，

4、所以本節(jié)探索僅用函數(shù)值的組合去尋求導(dǎo)數(shù)值的近似值，這是很自然的。6.1 Approximating The Derivative導(dǎo)數(shù)的近似值6.1.1 The Limit of the Difference Quotient差商的極限 考慮逼近函數(shù)的導(dǎo)數(shù)的數(shù)值過程：這個方法看起來很直接：選取序列，使得，計算下面序列的極限：讀者可能注意到我們只需要計算序列中的有限項(xiàng)，將作為答案。這樣的問題經(jīng)常被提出：為什么計算？同樣會問：如何取使得是的最佳逼近？為回答這個問題，我們再看一個例子，會發(fā)現(xiàn)為什么這不是件簡單的事。 例如，考慮函數(shù)，分別使用步長和去構(gòu)造點(diǎn)和點(diǎn)之間的割線。當(dāng)變小時，割線逼近于切線，如圖6

5、.2所示。盡管圖6.2給出了處理過程的可視化表示，我們必須計算時的值去獲得可以接受的數(shù)值解，而當(dāng)取這個值時割線與切線的圖像已不可分辨。上面處理為的過程，目的就是一個將連續(xù)問題離散化，把問題轉(zhuǎn)化為計算機(jī)可以處理的離散數(shù)據(jù)輸入。下面問題是計算不能無休止地進(jìn)行，只要達(dá)到一定的精度要求計算就可以停下來，也就是文中提到的有限項(xiàng)。從下面的例子我們還將看到，這不僅僅是考慮計算必須停止所要求的，也是因?yàn)椴⒎窃叫≡胶?，太小甚至?xí)U(kuò)大誤差。例6.1 設(shè)且。使用步長，其中，計算差商。精度為小數(shù)點(diǎn)后九位。 在計算中所用到的和的值如表6.1所示。 最大值不能得到好的近似，因?yàn)椴介L太大，使得兩點(diǎn)分割太遠(yuǎn)，差商是經(jīng)過這兩點(diǎn)

6、的割線的斜率，不能很好地近似切線。當(dāng)以小數(shù)點(diǎn)后9位的精度計算公式時，根據(jù)得到近似值，根據(jù)得到近似值。如果太小，那么算得的與就非常接近。作差時就會導(dǎo)致精度損失，因?yàn)閮蓚€量幾乎相等了。值太小了以至于存儲的值與是相等的，因此計算的差商是0。在例6.1中，極限值是?？梢钥吹?，取值可以給出最佳逼近，。從例6.1看出，找到方程的數(shù)值極限并不容易。這個差商序列一開始向收斂，是最接近的，然后開始遠(yuǎn)離。程序6.1中表明，當(dāng)時序列中的項(xiàng)停止計算。這是為了在項(xiàng)開始遠(yuǎn)離極限值之前，決定最佳逼近。將這個判定條件應(yīng)用到例6.1中，有，因此是我們選擇的答案?，F(xiàn)在我們開始討論如何根據(jù)較大的值得到合理精度的近似值公式。也就是說

7、，太大，割線不能很好地近似切線，但是也不能太小，否則由于精度損失，反而使誤差增大。上面的例子已經(jīng)顯示出，求的極限的數(shù)值近似解時，選擇適當(dāng)?shù)牡闹匾?。在表?.1中最后一列看到，在時最接近真實(shí)值，然后逐漸遠(yuǎn)離。其實(shí)，將相鄰的兩兩作差，可以發(fā)現(xiàn)，這個差值是由大變小，又由小變大的，所以必定存在，使得，這個就是所求。在實(shí)際的程序操作中，便可把上面的不等式，作為迭代運(yùn)算停止的判定條件。說完了存在性，那么為什么要以作為判定條件呢？這不禁讓人想到極限收斂的Cauchy判別方法，它的要義是要求當(dāng)自變量值在一定范圍內(nèi)時，函數(shù)值無限接近，也可以說函數(shù)值的差要任意小。我們在計算精度的限制下，不可能做到函數(shù)值之差任意

8、小，但起碼可以保證找到函數(shù)值差別最小的，這就是上面的，這樣這個就是最佳逼近所需的。6.1.2 The Central-difference Formulas中心差分公式如果函數(shù)在點(diǎn)的左邊和右邊的值可計算，則最佳二點(diǎn)公式包含兩邊的兩個對稱的橫坐標(biāo)。定理6.1（精度為的中心差分公式）設(shè)，且，則進(jìn)一步，存在數(shù)，滿足其中項(xiàng)稱為截斷誤差（the truncation error）。證明 從關(guān)于的二階泰勒展開式開始，寫出與：用式減去式，結(jié)果是因?yàn)槭沁B續(xù)的，根據(jù)中值定理可以找到值，使得將此式代入，得上式中右邊的第一項(xiàng)是中心差值公式，第二項(xiàng)是截斷誤差。定理證畢。假設(shè)三階導(dǎo)數(shù)的值變化不快，則式中的截斷誤差以與同

9、樣的方式趨近于0，表示為。當(dāng)用計算機(jī)進(jìn)行計算是，不宜將選得太小。因此，求解近似值的公式具有精度為的截斷誤差項(xiàng)是有用的。定理6.2（精度為的中心差分公式）設(shè)，且，則而且存在數(shù)，滿足其中。證明 一種導(dǎo)出公式式的方法如下。從關(guān)于四階泰勒展開式開始，作和的差：然后，使用步長代替，寫出以下近似：接下來用方程的8倍減去。包含的項(xiàng)就被消去了，得到如果的符號或是正或是負(fù)，而且它的值變化不快的話，在區(qū)間上可以找到一個值，使得把式代入，得到計算的結(jié)果式中右端的第一項(xiàng)是中心差分公式，第二項(xiàng)是截斷誤差。定理證畢。 假設(shè)在上有界，中的截斷誤差以與同樣的方式趨近于0，用符號表示。現(xiàn)在來比較一下公式和。假設(shè)有五階連續(xù)導(dǎo)數(shù)，

10、與大約相等。這樣四階公式的截斷誤差是的，比二階公式中的截斷誤差更快地趨近于0，這樣可以使用更大的步長。 回顧一下前面介紹的兩個定理以及得到的差分公式。在定理6.1中，用式減去式，消去了。同樣的，在定理6.2中，用方程的8倍減去。包含的項(xiàng)就被消去了。剩下的部分，根據(jù)中值定理并作必要的假設(shè)，合并整理出僅含的部分和誤差項(xiàng)部分。這樣得到的的表達(dá)式中就不包含導(dǎo)數(shù)信息，實(shí)現(xiàn)了僅用函數(shù)值來近似。定理6.2得到的四階公式的截斷誤差是的，比定理6.1得到的二階公式中的截斷誤差更快地趨近于0，這樣就可以使用較大的步長值以得到合理精度下的近似值。下面的例子，就在不同的步長下，分別采用公式和近似導(dǎo)數(shù)值，很好地說明了這

11、一點(diǎn)。例6.2 令步長取，使用公式和，計算的近似值。在所有的計算中均采用小數(shù)點(diǎn)后九位的精度。與準(zhǔn)確值比較。 取使用公式，得到取使用公式，得到由公式和公式得到的逼近值的誤差分別是-0.000011941和0.000000017。此例中，當(dāng)取時，公式比公式近似的效果要好。誤差分析將解釋這個例子并闡明為什么會這樣。其他的計算結(jié)果在表6.2中列出。再來分析一下由公式和公式得到的數(shù)值微分的結(jié)果。首先，用式進(jìn)行逼近時，在時，誤差最小。而用式進(jìn)行逼近時，在時，誤差就可以達(dá)到最小。換言之，使用式逼近導(dǎo)數(shù)值時，可以取較大步長實(shí)現(xiàn)要求精度。其次，可以看到，用式得到的最佳近似值的誤差即最小誤差是-0.0000000

12、91，而利用式得到的最佳近似值的誤差僅為0.000000017，所以說，即便是在不考慮步長的優(yōu)勢下，公式也比公式的近似效果要好。下面一個小節(jié)的工作就是將舍入誤差考慮進(jìn)去對誤差進(jìn)行分析，并通過考察總誤差取最小以找到步長的最優(yōu)值，討論的只是縱向最優(yōu)，沒有進(jìn)行公式和公式的橫向比較。6.1.3 Error Analysis and Optimum Step Size誤差分析和優(yōu)化步長關(guān)于數(shù)值微分的一個重要課題是研究計算機(jī)的舍入誤差。下面將對此進(jìn)行更深入的分析。假設(shè)使用計算機(jī)進(jìn)行數(shù)值計算，而且有和其中和是數(shù)值和的近似值，和分別是相關(guān)的舍入誤差。下面的結(jié)論說明了數(shù)值微分中誤差分析的復(fù)雜性。簡單的進(jìn)行一下符

13、號的說明，與中的的表示源于自變量，即，這里只是一種符號表示，對于往更高階推廣，使用更方便一些。推論6.1設(shè)函數(shù)滿足定理6.1中的假設(shè)，并利用計算公式誤差分析可通過如下的方程進(jìn)行解釋其中這里的總誤差項(xiàng)（the total error term）是舍入誤差（round-off error）與截斷誤差（truncation error）之和。推論6.1設(shè)函數(shù)滿足定理6.1中的假設(shè)，且進(jìn)行數(shù)值計算。如果，且，則式右邊最小時的值為 在式右邊關(guān)于求導(dǎo)，令導(dǎo)數(shù)為0，得到上式中的最優(yōu)值。當(dāng)較小時，式中包含的部分相對較大。在例6.2中，當(dāng)時，就會碰到這種情況。舍入誤差為截斷誤差為式中的誤差項(xiàng)為實(shí)際上，當(dāng)時的導(dǎo)數(shù)

14、近似值可用下式進(jìn)行計算顯然少了4位有效數(shù)字。誤差是-0.000003909，接近預(yù)計的誤差-0.000004001。在這一部分，使用例6.2，說明了上面兩個推論的結(jié)果，并表明當(dāng)時，舍入誤差相對較大。上面和 對應(yīng)的舍入誤差和，數(shù)量級都在，的數(shù)量級是，以致總的舍入誤差的數(shù)量級為，而截斷誤差項(xiàng)的數(shù)量級為，這樣舍入誤差相對于截斷誤差較大。而當(dāng)時，和 對應(yīng)的舍入誤差和，數(shù)量級仍是，的數(shù)量級是，這樣總的舍入誤差與截斷誤差項(xiàng)的數(shù)量級均為，在同一數(shù)量級上！此時的最優(yōu)，不只是減少舍入誤差或者截斷誤差，實(shí)際上保證了舍入誤差與截斷誤差和的最小，是一種博弈的結(jié)果。當(dāng)將式用于例6.2時，可用界和值計算舍入誤差。的優(yōu)化值

15、為0.001144714。步長時，最接近于優(yōu)化值0.001144714，而且在包含式的4個選項(xiàng)中，它給出了的最佳近似值（如表6.2和圖6.3所示）。 對于式的誤差分析跟上面的類似。設(shè)用計算機(jī)進(jìn)行數(shù)值計算，并有函數(shù)。推論6.2設(shè)函數(shù)滿足定理6.2的假設(shè)，使用計算公式誤差分析可用下列方程進(jìn)行解釋： 其中這里的總誤差項(xiàng)是舍入誤差與截斷誤差的和。推論6.2設(shè)函數(shù)滿足定理6.2的假設(shè)，且進(jìn)行數(shù)值計算。如果且，則式最小時的值為當(dāng)將式用于例6.2時，可用界和值計算舍入誤差。的優(yōu)化值為0.022388475，步長時最接近優(yōu)化值0.022388475，而且在包含式的4個選項(xiàng)中，給出了的最佳近似值（如表6.2和圖

16、6.4所示）。 這里可以將式和的結(jié)果再次橫向比較。圖6.3和6.4直觀地展示了式不僅可以去較大的步長，而且精度確實(shí)比式要高，即優(yōu)于。反映在圖像上就是圖線的最低點(diǎn)更低，離軸更近，而的標(biāo)度可以取更大。下面再次提到了篇首的做法，用插值多項(xiàng)式的微分近似函數(shù)的微分。通過另一種推導(dǎo)也可得到數(shù)值微分公式，即對插值多項(xiàng)式進(jìn)行微分。例如，經(jīng)過點(diǎn)，和的2次多項(xiàng)式的拉格朗日表達(dá)式為將它展開可得通過類似的計算可得到經(jīng)過點(diǎn)和4次多項(xiàng)式對這些多項(xiàng)式進(jìn)行微分，可得和，與表6.2中下面的值相符。和以及它們在點(diǎn)的切線分別如圖6.5和圖6.5所示。6.1.4 Richardsons Extrapolation Richardso

17、n外推法 這一節(jié)將重點(diǎn)研究式和式之間的關(guān)系。設(shè),且用和分別表示以和為步長根據(jù)式得到的的近似值，表示為和如果對式乘以4，并減去式，則可消去包含的項(xiàng)，結(jié)果為對式進(jìn)一步分解，可得式中最后一個表達(dá)式就是中心差分公式。例6.3 設(shè)，并利用式和式，取，說明式中的線性組合如何用來求出式給出的的近似值。精度為小數(shù)點(diǎn)后9位。利用式和式，取，可得和式中的線性組合為這與例6.2中直接用式所得的的近似值相同。這種從低階公式中推導(dǎo)出求解高階導(dǎo)數(shù)的方法稱為外推法。（這句話表述不夠清楚，意思是將誤差階由提高到的方法，稱為外推法）這里已經(jīng)看到了如何使用步長和消去包含的項(xiàng)。為了說明如何消去，用和分別表示以和為步長根據(jù)式得到的精度為的的近似值。則近似值可表示為和設(shè)只為正或負(fù)值，而且變化不快，則可用假設(shè)來消去式和式中的，結(jié)果為下面的結(jié)論描述了提高計算精度的一般形式。定理6.3（Richardson外推法）設(shè)的兩個精度為的近似值分別為和，而且它們滿足：和這樣可得到改進(jìn)的近似值表達(dá)式 本節(jié)由用插值多項(xiàng)式的微分來近似函數(shù)微分引入，探討了如何