數(shù)值計算_第7章數(shù)值微分和數(shù)值積分

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

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

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

4、60;得中心差商的截斷誤差：            差商的幾何意義微積分中的極限定義，表示在處切線的斜率，即圖7.1中直線的斜率；差商表示過和兩點直線 的斜率，是一條過的割線?？梢姅?shù)值微分是用近似值內(nèi)接弦的斜率代替準確值切線的斜率。圖7.1 微商與差商示意圖例7.1 給出下列數(shù)據(jù)，計算，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

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

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

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

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

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

10、造數(shù)值積分公式的過程。怎樣判斷數(shù)值積分公式的效果？代數(shù)精度是衡量數(shù)值積分公式優(yōu)劣的重要標準之一。定義7.1 (代數(shù)精度)記上以為積分節(jié)點的數(shù)值積分公式為若滿足而，則稱具有階代數(shù)精度。由此可知當具有階代數(shù)精度時，對任意的階多項式都有 。7.2.1 插值型數(shù)值積分對給定的被積函數(shù)在上的點列作拉格朗日插值多項式，以近似計算，即記，則有數(shù)值積分誤差，也就是對插值誤差的積分值 或?qū)σ话愕暮瘮?shù)，但若是一個不高于次的多項式，由于，而有。因此，階插值多項式型式的數(shù)值積分公式至少有階代數(shù)精度。例7.4 建立上節(jié)點為的數(shù)值積分公式。解：由得得以數(shù)值積分公式：7.2.2 牛頓-柯特斯（Newton-Cot

11、es）積分把積分區(qū)間分成等分，記步長為，取等分點作為數(shù)值積分節(jié)點，構(gòu)造拉格朗日插值多項式，取由此得到的數(shù)值積分稱為牛頓柯特斯積分。下面可以看到，牛頓柯特斯積分系數(shù)和積分節(jié)點以及積分區(qū)間無直接關(guān)系，系數(shù)固定而易于計算。梯形積分以和為插值節(jié)點構(gòu)造線性函數(shù)，有那么，提取公因子后，得到牛頓柯特斯積分的組合系數(shù)：，它們已與積分區(qū)間沒有任何關(guān)系了。記           （7.6）稱為梯形積分公式。它的幾何意義是用梯形面積近似代替積分值（圖7.2）。圖7.2 梯形積分怎樣確定梯形積分公式的代數(shù)精度？我們可以

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

13、。因此可設(shè)一個三次多項式滿足條件：，計算得到誤差為：  于是有但   故辛普森求積公式的截斷誤差：（7.9）牛頓-柯特斯積分系數(shù)等分區(qū)間，取等分點為積分節(jié)點，其中。以為插值節(jié)點構(gòu)造插值函數(shù)。 其中令，代入上式得（7.10）這里稱為牛頓-柯特斯系數(shù)可見在取等距節(jié)點時，積分系數(shù)與積分節(jié)點和積分區(qū)間無直接關(guān)系，只與插值的節(jié)點總數(shù)有關(guān)，而在例7.3中的積分系數(shù)是待定系數(shù)，這就簡化了數(shù)值積分公式，而不必對每一組插值節(jié)點xi都要計算一組相應的積分系數(shù)ai。在公式（7.10）中取=1，可算出梯形積分系數(shù)；取=2，可算出辛普森積分系數(shù)。在表7.1中列出從1到6的牛頓-柯

14、特斯系數(shù)。從表中可以看出牛頓-柯特斯系數(shù)具有對稱性。表7.11          2        3      4    5  67.2.3 求積公式的收斂性與穩(wěn)定性定義7.2 若，則稱求積公式是收斂的。定義中的包含了，通常都要求計算積分的求積公式是收斂的。穩(wěn)定性是研究計算公式當有誤差時，的誤差是否增長，現(xiàn)設(shè)，誤差記為。定義7.3 對任給，只要，就有則稱求積公式是穩(wěn)定的。定理7.1 若求積公式的系數(shù)，則求積公式是穩(wěn)定的。證明 由于,，故有于是對

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

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

17、。復化辛普森積分計算公式而，稱    （7.13）為復化辛普森積分公式，它是在上采用辛普森積分公式疊加而得。下面用圖7.5顯示復化辛普森積分計算公式中節(jié)點與系數(shù)的關(guān)系，取，在每個積分區(qū)間上提出因子后，三個節(jié)點的系數(shù)分別是1，4，1；將4個積分區(qū)間的系數(shù)按節(jié)點的位置累加，可以清楚地看到，首尾節(jié)點的系數(shù)是1，奇數(shù)點的系數(shù)是4，偶數(shù)點的系數(shù)是2。141141141141142424241圖7.5  復化辛普森積分系數(shù)復化辛普森公式的截斷誤差設(shè)，在上的誤差為因此，即        &

18、#160; （7.14）與復化梯形公式類似，誤差的截斷誤差按照或者的速度下降?？梢宰C明，只要在上有界并黎曼可積，當分點無限增多時，復化辛普森公式收斂到積分。記，則有對任給的誤差控制小量，只要    或    就有。例7.5 求，計算中要求有5位有效數(shù)字。用復化梯形和復化辛普森求積公式的分點應取多少？解：      由復化梯形誤差公式得到：計算出，復化梯形公式至少要在.00等分n = 68。由復化辛普森誤差公式，有在復化辛普森公式中取 或。   &

19、#160;   7.3.3 復化積分的自動控制誤差算法復化積分的誤差公式表明，截斷誤差隨分點的增大而減小，對于給定的誤差量，用估計函數(shù)導數(shù)的界的方法可計算出。用誤差公式計算滿足精度的分點數(shù)，像是在做一道計算導數(shù)上界的微積分習題（如例7.5所示）。但是在實際運算中，一般難以估計出函數(shù)的各階導數(shù)界，也就無法確定分點數(shù)。在計算中常用誤差的事后估計方法，即用估計誤差。T2n (f )的計算公式對定積分，取分點，計算得取分點，計算得這里，?？梢钥吹?，的值是與新增分點的組合。取分點，計算得這里，。同理，計算時只要在的基礎(chǔ)上計算新增分點，的值再做組合，如圖7.6所示。圖7.6 與一般地，

20、每次總對前一次的小區(qū)間分半，分點加密一倍，并可充分利用老分點上的函數(shù)值，每次只需計算新增分點的和。對上等分， ，則有記上的中點為則  （7.15）其中。或        其中。類似地，可得積分節(jié)點為，的辛普森求積公式的關(guān)系式：        （7.16）其中：由誤差公式：由于，分別為及個點上的均值，可視，于是 上式表明的誤差大約是誤差的4倍?；?#160;      

21、   （7.17）由此得到啟發(fā)，對任給的誤差控制量，要，只需即可，而用作為控制手段簡單直接，序列在計算機上也不難實現(xiàn)。復化積分的算法描述從數(shù)值積分的誤差公式可以看到，截斷誤差隨分點的增長而減少，控制計算的精度也就是確定分點數(shù)。在計算中不用數(shù)值積分的誤差公式確定分點數(shù)的理論模式，而用作為控制，通過增加分點自動滿足精度的方法稱為數(shù)值積分公式的自動積分法。即在計算中構(gòu)造序列，直到或時停止計算，由分點數(shù)自動控制積分值的誤差，并取。下面描述復化數(shù)值積分公式的自動控制誤差算法，詳細程序和算例請看本章7.6節(jié)。1輸入：誤差控制精度e = eps；初始分點值 。2計算分點的復化梯

22、形積分T1=T2+100      /迭代計算中T1和T2分別表示和3while | T1T2|eT1 = T2H=Hn           /計算新增節(jié)點的值T2= (T1+H )/2H = h/2，n =2n    /將區(qū)間一分為二end while4輸出積分值T2。在自動控制誤差算法中初始分點值不宜過小，以防假收斂。7.3.4 龍貝格（Romberg）積分由前面得到的關(guān)系式（7.17），可以將 作為的修正值補充

23、到，即 （7.18）其結(jié)果是將梯形求積公式組合成辛普森求積公式，截斷誤差由提高到。這種手段稱為外推算法。外推算法在不增加計算量的前題下提高了誤差的精度，是計算方法中一種常用手法。不妨對再做一次組合。由得到          （7.19）復化辛普森公式組成復化柯特斯公式，其截斷誤差是。同理對柯特斯公式進行組合： 得到具有7次代數(shù)精度和截斷誤差是 的龍貝格公式：還可以繼續(xù)對做上去。為了便于在計算機上實現(xiàn)龍貝格算法，將統(tǒng)一用表示，列標 分別表示梯形、辛普森、柯特斯積分，行標表示分點數(shù)或步長j。龍貝格計

24、算公式：對每一個從2做到，一直做到小于給定控制精度停止計算。龍貝格算法龍貝格算法按表7.2元素的行序進行運算， 在計算中每個元素只用到上一行和本行的元素。對上面的算法進一步優(yōu)化，對每k行可將計算定義在兩行元素之間，令；在每計算一行元素后，要將。表7.2 龍貝格算法計算元素順序表            1輸入?yún)^(qū)間端點 ，精度控制值，循環(huán)次數(shù)，定義函數(shù)，??；2；3for    =2  to          4輸出。7.4重積分計

25、算在微積分中計算二重積分是用化為累次積分的方法進行的。計算二重數(shù)值積分也是計算累次數(shù)值積分的過程。為了簡化問題，我們僅討論矩形域上的二重積分。有很多非矩形域上的二重積分可作變換將其轉(zhuǎn)換到矩形域上。                （7.20）其中：是常數(shù)，在上連續(xù)。像在微積分中一樣，將二重積分化為累次積分：      （7.21）或     

26、0;      二重積分的復化梯形公式對區(qū)間和分別選取正整數(shù)，在軸和軸上分別有步用復經(jīng)梯形公式計算，計算中將當作常數(shù)，有 （7.22）再將當作常數(shù)，在方向上計算式（7.23）中每一項的積分，有=則+=積分區(qū)域的4個角點的系數(shù)是1/4，4個邊界的系數(shù)是1/2，內(nèi)部節(jié)點的系數(shù)是1。誤差：和在積分區(qū)間內(nèi)。例子7.6 用復化梯形公式計算二重積分，取。解：如表7.3所示：表7.3 數(shù)值表的數(shù)值如表7.4所示：表7.4 數(shù)值表=0.873601的準確值是0.886176。二重復化辛普森求積公式對區(qū)間和分別等分和等分，在軸和軸上分別有步長 

27、;均為偶數(shù)類似于二重復化梯形公式推導，得記  誤差： 和在積分區(qū)間內(nèi)。按例7.6的化分區(qū)間，的值如表7.5所示：表7.57.5  高斯（Gauss）型積分公式介紹對插值型積分公式在牛頓-柯特斯積分公式中要求節(jié)點是等距的，其優(yōu)點是計算積分系數(shù)的公式規(guī)則相同，缺點是制約了求積公式的代數(shù)精度。可以證明：當節(jié)點個數(shù)為偶數(shù)時，求積公式具有階的代數(shù)精度；當節(jié)點個數(shù)為奇數(shù)時，求積公式具有階的代數(shù)精度。如果我們不預先指定求積節(jié)點的位置，和權(quán)系數(shù)都作為待定的常數(shù)，能否適當?shù)卮_定它們，以提高積分公式的代數(shù)精度？回答是肯定的。個待定常參數(shù)，需要個方程來確定，取一個函數(shù)組：，這一組函數(shù)構(gòu)

28、成了次多項式的基，任一小于等于次的多項式，都可以用這組函數(shù)的線性組合來表示。如果某一積分公式，對這組函數(shù)都能精確積分，則此積分公式就有次代數(shù)精度。例7.7 計算求積系數(shù)和求積節(jié)點，使得至少具有3階代數(shù)精度。解：按照求積公式的代數(shù)精度定義，分別令，得方程組：解方程組得：  ， 求積公式：按例7.7的方式，構(gòu)造更高階的代數(shù)精度的求積公式，生成求積系數(shù)和求積節(jié)點的方程組并無困難，而求解該方程組則無一定的章法可循。一般地，通過計算正交多項式的零點作為求積節(jié)點。當取積分節(jié)點為正交多項式的零點時，則節(jié)點個數(shù)是的求積公式具有階的代數(shù)精度。并稱積分節(jié)點為正交多項式的零點的數(shù)值積分公式為高斯型積分公式

29、。為了一般性，考慮積分其中稱為權(quán)函數(shù)。當時，即是普通的積分。對于不同的權(quán)函數(shù)選定的節(jié)點也不相同。如何構(gòu)造高斯型積分公式呢？對給定的及權(quán)函數(shù)，由施密特（Schmidt）正交化過程作出正交多項式；解出正交多項式的個零點，這個個零點就是積分節(jié)點；以這些節(jié)點構(gòu)造插值多項式，計算積分系數(shù)其中是拉格朗日插值基函數(shù)。高斯型求積公式為高斯型積分公式的優(yōu)點是它的代數(shù)精度高，特別是對無窮區(qū)間或瑕積分更有效，但計算正交多項式的零點即積分節(jié)點有一定的工作量，好在數(shù)值學家們已算出一些特定的函數(shù)的積分節(jié)點和積分系數(shù)，在計算中我們可以查表直接得到這些數(shù)值。本章并不構(gòu)造各種高斯型積分公式，有關(guān)的詳細內(nèi)容請參考有關(guān)的教材。下面

30、給出上，取權(quán)函數(shù)的高斯型積分。取1，1上權(quán)函數(shù)的正交多項式為勒讓德（Legendre）多項式：高斯-勒讓德相應的積分節(jié)點和積分系數(shù)表如下：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要計算一般區(qū)間上的積分，

31、只需作變量代換則有，其中，這樣 仍可用高斯積分求積，即 例7.8 應用兩點高斯-勒讓德積分公式計算。解：I = (0.5773503)2cos (0.5773503) + (0.5773503)2 cos (0.5773503)=0.558608例7.9 應用兩點高斯-勒讓德積分公式計算。解：令，得到積分例7.10 證明不存在使求積公式的代數(shù)精度超過次。證明：只要能找到一個次多項式，使求積公式兩邊不相等即可。用反證法，假定存在求積系數(shù)和節(jié)點及使求積公式對任何次多項式精確成立?，F(xiàn)取代入求積公式左端得而公式右端，故右端與左端不相等，與假設(shè)矛盾，說明不存在次代數(shù)精度的求積公式，故高斯型求積

32、公式是具有最高代數(shù)精度的求積公式。7.6程序示例程序7.1  復化梯形公式的自動控制誤差算法。算法描述輸入的值，定義；   ；T1=T2+100while |T1T2|T1 = T2；H = h        ！或用復化梯形公式計算T2T2 = (T1 + H) /2；h= h/2；n = 2n；end while輸出T2。程序源碼/  Purpose：梯形公式的自動控制誤差算法      /# include &l

33、t;stdio.h># include <math.h># define f (x)  (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；   dou

34、ble h = (ba) / n；   T_n = (f (a) + f(b) /2；   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

35、=0；i =1；in；i+          H_n+=f (a+hi+h/2)；      H_n= h；      T2 = (T1 + H_n) /2h = h /2；n = n2；        while (fabs (T1T2)epsilon)       

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

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