1、第7章  數值微分和數值積分7.1  數值微分 差商與數值微分當函數是以離散點列給出時，當函數的表達式過于復雜時，常用數值微分近似計算的導數。在微積分中，導數表示函數在某點上的瞬時變化率，它是平均變化率的極限；在幾何上可解釋為曲線的斜率；在物理上可解釋為物體變化的速率。以下是導數的三種定義形式：             （7.1）在微積分中，用差商的極限定義導數；在數值計算中返璞歸真，導數取用差商（平均變化率）作為其近似值。最簡單的計算數值微分的方法是用函

2、數的差商近似函數的導數，即取極限的近似值。下面是與式（7.1）相應的三種差商形式的數值微分公式以及相應的截斷誤差。向前差商用向前差商（平均變化率）近似導數有：                      （7.2）其中的位置在的前面，因此稱為向前差商。同理可得向后差商、中心差商的定義。由泰勒展開 得向前差商的截斷誤差： 向后差商用向后差商近似導數有： 

3、0;                   （7.3）與計算向前差商的方法類似，由泰勒展開得向后差商的截斷誤差：中心差商用中心差商（平均變化率）近似導數有：                    （7.4）由泰勒展開 得中

4、心差商的截斷誤差：            差商的幾何意義微積分中的極限定義，表示在處切線的斜率，即圖7.1中直線的斜率；差商表示過和兩點直線的斜率，是一條過的割線。可見數值微分是用近似值內接弦的斜率代替準確值切線的斜率。圖7.1 微商與差商示意圖例7.1 給出下列數據，計算，0.020.040.060.80.105.065.075.0655.055.055解：（5.075.06）/（0.040.02）= 0.5（5.055.07）/（0.080.04）= 0.5（5.055.055）/

5、（0.080.10）= 0.25（(0.10) (0.06)）/（0.100.06）= 18.75設定最佳步長在計算數值導數時，它的誤差由截斷誤差和舍入差兩部分組成。用差商或插值公式近似導數產生截斷誤差，由原始值的數值近似產生舍入誤差。在差商計算中，從截斷誤差的逼近值的角度看，越小，則誤差也越小；但是太小的會帶來較大的舍入誤差。怎樣選擇最佳步長，使截斷誤差與舍入誤差之和最小呢？一般對計算導數的近似公式進行分析可得到誤差的表示式，以中心差商為例，截斷誤差不超過而舍入誤差可用量估計（證明略），其中是函數的原始值的絕對誤差限，總誤差為當時，總誤差達到最小值，即    &

6、#160;                 （*）可以看到用誤差的表達式確定步長，難度較大，難以實際操作。通常用事后估計方法選取步長，例如，記為步長等于的差商計算公式，給定誤差界，當時，就是合適的步長。例7.2 對函數，取不同的步長計算，觀察誤差變化規律，從而確定最佳步長。解：誤差誤差0.100.090.080.070.063.16303.16223.16133.16073.16000.00480.00400.00310.00

7、250.00180.050.040.030.020.013.15903.15883.15833.15753.15500.00080.00060.00010.00070.0032表中數據顯示，當步長從0.10減少到0.03時，數值微分誤差的絕對值從0.0048減少到0.0001，而隨著的進一步減少，誤差的絕對值又有所反彈，表明當步長小于0.03時，舍入誤差起了主要作用。在實際計算中是無法得到誤差的準確數值的，這時以最小為標準確定步長，本例中取= 0.04。 插值型數值微分對于給定的的函數表，建立插值函數，用插值函數的導數近似函數的導數。設為上的節點，給定，以為插值點構造插值多項式，以的各階導數近

8、似的相應階的導數，即當時，            （7.5）誤差項為：例7.3 給定，并有，計算。解：作過的插值多項式：將代入 得三點端點公式和三點中點公式：    利用泰勒（Taylor）展開進行比較和分析，可得三點公式的截斷誤差是。類似地，可得到五點中點公式和五點端點公式：     樣條插值數值微分把離散點按大小排列成，用關系式構造插值點的樣條函數：當則當時，可用計算導數。7.2數值積分在微積分中用牛頓

9、萊布尼茲（Newton-Leibniz）公式計算連續函數的定積分： 但是，當被積函數是以點列的形式給出時，當被積函數的原函數難以得到時，例如，則無法用牛頓萊布尼茲積分公式計算。有時當被積函數的原函數過于復雜時，也不宜套用積分公式計算積分，而應采用數值積分公式計算定積分。在微積分中，定積分是黎曼（Rimann）和的極限，它是分割小區間長度趨于零時的極限，即在數值積分公式中，只能用有限項的和近似上面的極限，通常由函數在離散點函數值的線性組合形式給出。記，在本章中，用表示精確積分值，用表示近似積分值，稱為積分節點，稱為積分系數。確定中積分系數的過程就是構造數值積分公式的過程。怎樣判斷數值積

10、分公式的效果？代數精度是衡量數值積分公式優劣的重要標準之一。定義7.1 (代數精度)記上以為積分節點的數值積分公式為若滿足而，則稱具有階代數精度。由此可知當具有階代數精度時，對任意的階多項式都有 。 插值型數值積分對給定的被積函數在上的點列作拉格朗日插值多項式，以近似計算，即記，則有數值積分誤差，也就是對插值誤差的積分值 或對一般的函數，但若是一個不高于次的多項式，由于，而有。因此，階插值多項式型式的數值積分公式至少有階代數精度。例7.4 建立上節點為的數值積分公式。解：由得得以數值積分公式： 牛頓-柯特斯（Newton-Cotes）積分把積分區間分成等分，記步長為，取等分點作為數值

11、積分節點，構造拉格朗日插值多項式，取由此得到的數值積分稱為牛頓柯特斯積分。下面可以看到，牛頓柯特斯積分系數和積分節點以及積分區間無直接關系，系數固定而易于計算。梯形積分以和為插值節點構造線性函數，有那么，提取公因子后，得到牛頓柯特斯積分的組合系數：，它們已與積分區間沒有任何關系了。記           （7.6）稱為梯形積分公式。它的幾何意義是用梯形面積近似代替積分值（圖7.2）。圖7.2 梯形積分怎樣確定梯形積分公式的代數精度？我們可以取驗證取時，有即：取時，有即：取時，有得梯形求積公式具有

12、一階代數精度。梯形積分公式的誤差：由  得  因為在上不變號，由積分中值定理得到梯形求積公式的截斷誤差：   （7.7）辛普森（Simpson）積分對區間作二等分，記。以和為插值節點構造二次插值函數，那么，有計算得到積分組合系數：，。   （7.8）S (f)稱為辛普森或拋物線積分公式。它的幾何意義是用過三點的拋物線面積近似代替積分的曲邊面積（圖7.3）。圖7.3 拋物線積分面積分別將代入到和中，可以得到表明辛普森公式對于次數不超過三次的多項式準確成立，具有三階代數精度。因此可設一個三次多項式滿足條件：，計算得到誤差為：&#

13、160; 于是有但   故辛普森求積公式的截斷誤差：（7.9）牛頓-柯特斯積分系數等分區間，取等分點為積分節點，其中。以為插值節點構造插值函數。 其中令，代入上式得（7.10）這里稱為牛頓-柯特斯系數可見在取等距節點時，積分系數與積分節點和積分區間無直接關系，只與插值的節點總數有關，而在例7.3中的積分系數是待定系數，這就簡化了數值積分公式，而不必對每一組插值節點xi都要計算一組相應的積分系數ai。在公式（7.10）中取=1，可算出梯形積分系數；取=2，可算出辛普森積分系數。在表7.1中列出從1到6的牛頓-柯特斯系數。從表中可以看出牛頓-柯特斯系數具有對稱性。表7

14、.11          2        3      4    5  6 求積公式的收斂性與穩定性定義7.2 若，則稱求積公式是收斂的。定義中的包含了，通常都要求計算積分的求積公式是收斂的。穩定性是研究計算公式當有誤差時，的誤差是否增長，現設，誤差記為。定義7.3 對任給，只要，就有則稱求積公式是穩定的。定理7.1 若求積公式的系數，則求積公式是穩定的。證明 由于,，故有于是對，只要，就有故求積公式是穩定的。7.3復化數值積分由插值的龍格現象

15、可知，高階牛頓-柯特斯積分不能保證等距數值積分系列的收斂性，同時可證（略）高階牛頓-柯特斯積分的計算是不穩定的。因此，實際計算中常用低階復化梯形等積分公式。 復化梯形積分把積分區間分割成若干小區間，在每個小區間上用梯形積分公式，再將這些小區間上的數值積分累加起來，稱為復化梯形公式。復化梯形公式用若干個小梯形面積逼近積分比用一個大梯形公式效果顯然更好，如圖7.4所示。這種作法使我們想起定積分定義，即它為被積函數無限分割的代數和。這也正是計算定積分最樸素的算法。圖7.4 復化梯形公式積分視圖復化梯形積分計算公式對作等距分割，有，于是在上，則有 記等分的復化梯形公式為，有 （7.

16、11）復化梯形公式截斷誤差由，根據均值定理，當時，存在，有，于是   （7.12） 由此看到復化梯形公式的截斷誤差按照或者的速度下降，事實上，可以證明，只要在上有界并黎曼可積，當分點無限增多時，復化梯形公式收斂到積分。記 ，則有對于任給的誤差控制小量，有   或   就有，式中表示取其最大整數。 復化辛普森積分把積分區間分成偶數等分，記，其中是節點總數，是積分子區間的總數。記，在每個區間上用辛普森數值積分公式計算，則得到復化辛普森公式，記為。復化辛普森積分計算公式而，稱    （7.1

17、3）為復化辛普森積分公式，它是在上采用辛普森積分公式疊加而得。下面用圖7.5顯示復化辛普森積分計算公式中節點與系數的關系，取，在每個積分區間上提出因子后，三個節點的系數分別是1，4，1；將4個積分區間的系數按節點的位置累加，可以清楚地看到，首尾節點的系數是1，奇數點的系數是4，偶數點的系數是2。141141141141142424241圖7.5  復化辛普森積分系數復化辛普森公式的截斷誤差設，在上的誤差為因此，即          （7.14）與復化梯形公式類似，誤差的截斷誤差按照或者的速度下降

18、。可以證明，只要在上有界并黎曼可積，當分點無限增多時，復化辛普森公式收斂到積分。記，則有對任給的誤差控制小量，只要    或    就有。例7.5 求，計算中要求有5位有效數字。用復化梯形和復化辛普森求積公式的分點應取多少？解：      由復化梯形誤差公式得到：計算出，復化梯形公式至少要在.00等分n = 68。由復化辛普森誤差公式，有在復化辛普森公式中取 或。       復化積分的自動控制誤差算法復化積分的誤差公式表明，

19、截斷誤差隨分點的增大而減小，對于給定的誤差量，用估計函數導數的界的方法可計算出。用誤差公式計算滿足精度的分點數，像是在做一道計算導數上界的微積分習題（如例7.5所示）。但是在實際運算中，一般難以估計出函數的各階導數界，也就無法確定分點數。在計算中常用誤差的事后估計方法，即用估計誤差。T2n (f )的計算公式對定積分，取分點，計算得取分點，計算得這里，。可以看到，的值是與新增分點的組合。取分點，計算得這里，。同理，計算時只要在的基礎上計算新增分點，的值再做組合，如圖7.6所示。圖7.6 與一般地，每次總對前一次的小區間分半，分點加密一倍，并可充分利用老分點上的函數值，每次只需計算新增分點的和。

20、對上等分， ，則有記上的中點為則  （7.15）其中。或        其中。類似地，可得積分節點為，的辛普森求積公式的關系式：        （7.16）其中：由誤差公式：由于，分別為及個點上的均值，可視，于是 上式表明的誤差大約是誤差的4倍。或          （7.17）由此得到啟發，對任給的誤差控制量，要，只需即可，而

21、用作為控制手段簡單直接，序列在計算機上也不難實現。復化積分的算法描述從數值積分的誤差公式可以看到，截斷誤差隨分點的增長而減少，控制計算的精度也就是確定分點數。在計算中不用數值積分的誤差公式確定分點數的理論模式，而用作為控制，通過增加分點自動滿足精度的方法稱為數值積分公式的自動積分法。即在計算中構造序列，直到或時停止計算，由分點數自動控制積分值的誤差，并取。下面描述復化數值積分公式的自動控制誤差算法，詳細程序和算例請看本章7.6節。1輸入：誤差控制精度e = eps；初始分點值。2計算分點的復化梯形積分T1=T2+100      /迭代計算中T

22、1和T2分別表示和3while | T1T2|eT1 = T2H=Hn           /計算新增節點的值T2= (T1+H )/2H = h/2，n =2n    /將區間一分為二end while4輸出積分值T2。在自動控制誤差算法中初始分點值不宜過小，以防假收斂。 龍貝格（Romberg）積分由前面得到的關系式（7.17），可以將作為的修正值補充到，即 （7.18）其結果是將梯形求積公式組合成辛普森求積公式，截斷誤差由提高到。這種手段稱為外推算法

23、。外推算法在不增加計算量的前題下提高了誤差的精度，是計算方法中一種常用手法。不妨對再做一次組合。由得到          （7.19）復化辛普森公式組成復化柯特斯公式，其截斷誤差是。同理對柯特斯公式進行組合： 得到具有7次代數精度和截斷誤差是的龍貝格公式：還可以繼續對做上去。為了便于在計算機上實現龍貝格算法，將統一用表示，列標 分別表示梯形、辛普森、柯特斯積分，行標表示分點數或步長j。龍貝格計算公式：對每一個從2做到，一直做到小于給定控制精度停止計算。龍貝格算法龍貝格算法按表7.2元素的行序進行運算，在計

24、算中每個元素只用到上一行和本行的元素。對上面的算法進一步優化，對每k行可將計算定義在兩行元素之間，令；在每計算一行元素后，要將。表7.2 龍貝格算法計算元素順序表            1輸入區間端點 ，精度控制值，循環次數，定義函數，取；2；3for    =2  to          4輸出。7.4重積分計算在微積分中計算二重積分是用化為累次積分的方法進行的。計算二重數值積分也是計算累次數值積分的過程。為了簡化問題，我們

25、僅討論矩形域上的二重積分。有很多非矩形域上的二重積分可作變換將其轉換到矩形域上。                （7.20）其中：是常數，在上連續。像在微積分中一樣，將二重積分化為累次積分：      （7.21）或            二重積分的復化梯形公式對區間和分別選取正整數，在軸

26、和軸上分別有步用復經梯形公式計算，計算中將當作常數，有 （7.22）再將當作常數，在方向上計算式（7.23）中每一項的積分，有=則+=積分區域的4個角點的系數是1/4，4個邊界的系數是1/2，內部節點的系數是1。誤差：和在積分區間內。例子7.6 用復化梯形公式計算二重積分，取。解：如表7.3所示：表7.3 數值表的數值如表7.4所示：表7.4 數值表=0.873601的準確值是0.886176。二重復化辛普森求積公式對區間和分別等分和等分，在軸和軸上分別有步長 均為偶數類似于二重復化梯形公式推導，得記  誤差： 和在積分區間內。按例7.6的化分區間，的值

27、如表7.5所示：表7.57.5  高斯（Gauss）型積分公式介紹對插值型積分公式在牛頓-柯特斯積分公式中要求節點是等距的，其優點是計算積分系數的公式規則相同，缺點是制約了求積公式的代數精度。可以證明：當節點個數為偶數時，求積公式具有階的代數精度；當節點個數為奇數時，求積公式具有階的代數精度。如果我們不預先指定求積節點的位置，和權系數都作為待定的常數，能否適當地確定它們，以提高積分公式的代數精度？回答是肯定的。個待定常參數，需要個方程來確定，取一個函數組：，這一組函數構成了次多項式的基，任一小于等于次的多項式，都可以用這組函數的線性組合來表示。如果某一積分公式，對這組函數都能精確積分

28、，則此積分公式就有次代數精度。例7.7 計算求積系數和求積節點，使得至少具有3階代數精度。解：按照求積公式的代數精度定義，分別令，得方程組：解方程組得：  ， 求積公式：按例7.7的方式，構造更高階的代數精度的求積公式，生成求積系數和求積節點的方程組并無困難，而求解該方程組則無一定的章法可循。一般地，通過計算正交多項式的零點作為求積節點。當取積分節點為正交多項式的零點時，則節點個數是的求積公式具有階的代數精度。并稱積分節點為正交多項式的零點的數值積分公式為高斯型積分公式。為了一般性，考慮積分其中稱為權函數。當時，即是普通的積分。對于不同的權函數選定的節點也不相同。如何構造高斯型積分公

29、式呢？對給定的及權函數，由施密特（Schmidt）正交化過程作出正交多項式；解出正交多項式的個零點，這個個零點就是積分節點；以這些節點構造插值多項式，計算積分系數其中是拉格朗日插值基函數。高斯型求積公式為高斯型積分公式的優點是它的代數精度高，特別是對無窮區間或瑕積分更有效，但計算正交多項式的零點即積分節點有一定的工作量，好在數值學家們已算出一些特定的函數的積分節點和積分系數，在計算中我們可以查表直接得到這些數值。本章并不構造各種高斯型積分公式，有關的詳細內容請參考有關的教材。下面給出上，取權函數的高斯型積分。取1，1上權函數的正交多項式為勒讓德（Legendre）多項式：高斯-勒讓德相應的積分

30、節點和積分系數表如下：2 =0.5773503， =0.57735031.0000000，1.00000004=0.8611363，=0.3399810 =0.3399810， =0.86113630.3478548，0.65214520.6521452，0.34785485x1= x5=0.9061798x2= x4=0.538469=0.00.2369269，0.23692690.4786287，0.47862870.5688889要計算一般區間上的積分，只需作變量代換則有，其中，這樣 仍可用高斯積分求積，即 例7.8 應用兩點高斯-勒讓德積分公式計算。解：I

31、 = (0.5773503)2cos (0.5773503) + (0.5773503)2 cos (0.5773503)=0.558608例7.9 應用兩點高斯-勒讓德積分公式計算。解：令，得到積分例7.10 證明不存在使求積公式的代數精度超過次。證明：只要能找到一個次多項式，使求積公式兩邊不相等即可。用反證法，假定存在求積系數和節點及使求積公式對任何次多項式精確成立。現取代入求積公式左端得而公式右端，故右端與左端不相等，與假設矛盾，說明不存在次代數精度的求積公式，故高斯型求積公式是具有最高代數精度的求積公式。7.6程序示例程序7.1  復化梯形公式的自動控制誤差算法。算法描述輸入

32、的值，定義；   ；T1=T2+100while |T1T2|T1 = T2；H = h        ！或用復化梯形公式計算T2T2 = (T1 + H) /2；h= h/2；n = 2n；end while輸出T2。程序源碼/  Purpose：梯形公式的自動控制誤差算法      /# include <stdio.h># include <math.h># define f (x) 

33、(sin (x)         / f(x)由define定義# define m  100# define a 1.0# define b 2.0# define epsilon 0.00001void main( )  int n = m；   int i；   double T n，H_n，T1，T2；   double h = (ba) / n；   T_n = (f (a) + f(b) /2； 

34、;  for (i =1；in；i+       T_n + = f(a+h i)；       T_n= h；   T2 = T_n；   T1 = T2+100；   do     T1 = T2；    for (i=0，H_n=0；i =1；in；i+       

35、60;  H_n+=f (a+hi+h/2)；      H_n= h；      T2 = (T1 + H_n) /2h = h /2；n = n2；        while (fabs (T1T2)epsilon)        printf (“T=% 1f n”，T2)；      

36、60;   / End of File計算實例對于，在區間上驗證梯形公式的自動控制誤差公式。程序輸入輸出對于不同，修改程序# define f(x)項，本例，初始，區間。）        T = 0.956447程序7.2  龍貝格積分算法計算公式和算法描述第節所述。程序源碼/       Purpose：龍貝格算法        /# include<

37、;stdio.h># include<math.h># define f (x)             (sin (x)# define N_H             20# define MAXREPT        10# define a                 1.0# define b