MATLAB第十二講-數(shù)值積分

數(shù)值積分21

數(shù)值積分概論

1.1數(shù)值積分的根本思想依據(jù)微積分根本定理，對于積分只要找到被積函數(shù)的原函數(shù)，便有下列牛頓-萊布尼茲(Newton-Leibniz)公式：

但對于以下情形：3〔1〕被積函數(shù)，諸如等，找不到用初等函數(shù)表示的原函數(shù)，或者即使能求得原函數(shù)但原函數(shù)的表達式非常復(fù)雜，計算困難；

（2）當是由測量或數(shù)值計算給出的一張數(shù)據(jù)表時，牛頓-萊布尼茲公式也不能直接運用.

因此有必要研究積分的數(shù)值計算問題.

由積分中值定理知，在積分區(qū)間內(nèi)存在一點ξ，成立4就是說，底為而高為的矩形面積恰等于所求曲邊梯形的面積(圖4-1).圖4-15

問題在于點ξ的具體位置一般是不知道的，因而難以準確算出的值.

將稱為區(qū)間上的平均高度.

這樣只要對平均高度提供一種算法，相應(yīng)地便獲得一種數(shù)值求積方法.

用兩端點“高度”與的算術(shù)平均作為平均高度的近似值，這樣導(dǎo)出的求積公式（1.1）是梯形公式(幾何意義參看圖4-2).

6圖4-2用區(qū)間中點的“高度”近似地取代平均高度，那么又可導(dǎo)出所謂中矩形公式(簡稱矩形公式)（1.2）左矩形公式梯形公式中矩形公式右矩形公式簡單的數(shù)值積分公式左/中/右矩形公式梯形公式9

一般地，可以在區(qū)間上適當選取某些節(jié)點，然后用加權(quán)平均得到平均高度的近似值，（1.3）式中稱為求積節(jié)點；稱為求積系數(shù)，亦稱伴隨節(jié)點

的權(quán).

權(quán)僅僅與節(jié)點的選取有關(guān)，而不依賴于被積函數(shù)的具體形式.這樣構(gòu)造出的求積公式具有以下形式：10

這類數(shù)值積分方法通常稱為機械求積，其特點是將積分求值問題歸結(jié)為函數(shù)值的計算，這就避開了牛頓-萊布尼茲公式需要尋求原函數(shù)的困難.

111.2

代數(shù)精度的概念

定義1如果某個求積公式對于次數(shù)不超過的多項式均能準確地成立，但對于次多項式就不準確成立，則稱該求積公式具有次代數(shù)精度.

梯形公式(1.1)和矩形公式(1.2)均具有一次代數(shù)

精度.

數(shù)值求積是近似方法，為保證精度，自然希望求積公式對盡可能多的函數(shù)準確成立.12

欲使求積公式(1.3)具有次代數(shù)精度，則只要令它對都準確成立，就得到（1.4）13

如果事先選定求積節(jié)點，譬如，以區(qū)間的等距分點作為節(jié)點，這時取，求解方程組(1.4)即可確定求積系數(shù)，而使求積公式(1.3)至少具有次代數(shù)精度.

構(gòu)造形如(1.3)的求積公式，原則上是一個確定參數(shù)和的代數(shù)問題.（1.3）14

例如時，取，求積公式為在線性方程組〔1.4〕中令，那么得解得于是得這就是梯形公式，說明利用線性方程組〔1.4〕推出的求積公式，與用通過兩點與的直線近似曲線得到的結(jié)果是一致的.15當時〔1.4〕式的第3個式子不成立，因為所以梯形公式〔1.1〕的代數(shù)精度為1.在〔1.4〕中如果節(jié)點和系數(shù)都不確定，那么〔1.4〕就是關(guān)于及的個參數(shù)的非線性方程組，該方程組在時求解是很困難的.但在和時還是可以通過求解〔1.4〕得到相應(yīng)的求積公式的.16

如

，此時求積公式為其中，及為待定參數(shù).根據(jù)代數(shù)精度的定義可令，由〔1.4〕知于是

所得到的就是〔1.2〕式的中矩形公式.17再令，代入〔1.4〕的第3式有說明公式〔1.2〕對不精確成立，故它的代數(shù)精度為1.方程組〔1.4〕是根據(jù)形如〔1.3〕式的求積公式得到的，按照代數(shù)精度的定義，如果求積公式中除了還有在某些節(jié)點上的值，也同樣可得到相應(yīng)的求積公式.18

例1

給定形如的求積公式，試確定系數(shù)，使公式具有盡可能高的代數(shù)精度.

解根據(jù)題意可令分別代入求積公式使它精確成立

當時，得

當時，得19

當時，得解得,于是得

當時，而上式右端為，故公式對不精確成立，其代數(shù)精度為2.201.3

插值型的求積公式

設(shè)給定一組節(jié)點且已知函數(shù)在這些節(jié)點上的值，作插值函數(shù)

.取

作為積分的近似值，（1.5）值，這樣構(gòu)造出的求積公式稱為是插值型的，式中求積系數(shù)通過插值基函數(shù)積分得出21（1.6）求積公式的余項式中ξ依賴于，〔1.7〕其中22

當是次數(shù)不超過的多項式時，插值多項式就是函數(shù)本身，余項為零，

反之，如果求積公式(1.5)至少具有次代數(shù)精度，則它必定是插值型的.

事實上，這時公式(1.5)對于插值基函數(shù)應(yīng)準確成立，即有至少具有次代數(shù)精度.所以這時插值型求積公式23注意到上式右端實際上即等于，因而成立.這樣，有

定理1形如(1.5)的求積公式至少有次代數(shù)精度的充分必要條件是，它是插值型的.

（1.5）24假設(shè)求積公式(1.3)的代數(shù)精度為，那么由求積公式余項的表達式(1.7)可以證明余項形如〔1.8〕其中為不依賴于的待定參數(shù),結(jié)果說明當是次數(shù)小于等于的多項式時，由于，故此時，即求積公式(1.3)精確成立.

而當時，(1.8)的右端故可求得1.4

求積公式的余項25〔1.9〕代入余項(1.8)中可以得到更細致的余項表達式.

梯形公式(1.1)的代數(shù)精度為1，可以證明它的余項表達式為其中于是得到梯形公式(1.1)的余項為〔1.10〕26

對中矩形公式(1.2)，其代數(shù)精度為1，可以證明其中于是得到中矩形公式(1.2)的余項為〔1.11〕27

例2

求例1中求積公式的余項

解由于此求積公式的代數(shù)精度為2，故余項表達式為

.令，得，于是有故得282牛頓-柯特斯〔Newton-Cotes〕公式

2.1

柯特斯系數(shù)與辛普森公式

設(shè)將積分區(qū)間劃分為等分，選取等距節(jié)點構(gòu)造出的插值型求積公式（2.1）稱為牛頓-柯特斯公式，式中稱為柯特斯系數(shù).

按(1.6)式，引進變換步長那么利用等距節(jié)點的插值公式，有29（2.2）

當時，這時的求積公式就是梯形公式30

當時，按(2.2)式，相應(yīng)的求積公式是辛普森(Simpson)公式

（2.3）柯特斯系數(shù)為

31

的牛頓-柯特斯公式稱為柯特斯公式，其形式是

（2.4）這里

按〔2.2〕式，可構(gòu)造柯特斯系數(shù)表.3233

從柯特斯系數(shù)表看到時，柯特斯系數(shù)出現(xiàn)負值，特別地，假定于是有且那么有34它表明初始數(shù)據(jù)誤差將會引起計算結(jié)果誤差增大，即計算不穩(wěn)定，故的牛頓-柯特斯公式是不用的.352.2

偶階求積公式的代數(shù)精度

由定理1，階的牛頓-柯特斯公式至少具有次的代數(shù)精度.

先看辛普森公式(2.3)，它是二階牛頓-柯特斯公式，因此至少具有二次代數(shù)精度.

用進行檢驗，本節(jié)討論代數(shù)精度的進一步提高問題.

按辛普森公式計算得36這時有，即辛普森公式對次數(shù)不超過三次的多項式均能準確成立，而它對通常是不準確的，因此，辛普森公式實際上具有三次代數(shù)精度.另一方面，直接求積得

定理3當階為偶數(shù)時，牛頓-柯特斯公式(2.1)至少有次代數(shù)精度.

（2.1）372.3

辛普森公式的余項

對牛頓-柯特斯求積公式，通常只使用時的三個公式，時為梯形公式，余項為

為辛普森公式代數(shù)精度為3，可以證明余項表達式為其中由(1.9)及(2.3)可得38從而可得辛普森公式的余項為〔2.5〕

為柯特斯公式代數(shù)精度為5，可以證明余項〔2.6〕393

復(fù)合求積公式復(fù)合求積的根本思想是把積分區(qū)間分成假設(shè)干子區(qū)間(通常是等分)，再在每個子區(qū)間上用低階求積公式，目的是提高精度.3.1

復(fù)合梯形公式

將區(qū)間劃分為等分，分點在每個子區(qū)間上采用梯形公式(1.1)，那么得40（3.1）記

稱為復(fù)合梯形公式.

（3.2）41

由(1.10)，余項由于,且

所以使

于是復(fù)合梯形公式余項為

42（3.3）誤差是階，且當時有

即復(fù)合梯形公式是收斂的.

事實上只要設(shè)，就可以得到收斂性，因為只要將改寫為

43

此外，的求積系數(shù)為正，由定理2知復(fù)合梯形公式是穩(wěn)定的.

當時，上式右端括號內(nèi)的兩個和式均收斂到積分所以復(fù)化梯形公式(3.2)收斂.443.2

復(fù)合辛普森求積公式將區(qū)間分為等分，在每個子區(qū)間上采用辛普森公式(2.3)，假設(shè)記，那么得記〔3.5〕（3.4）（2.3）稱為復(fù)合辛普森求積公式.45

由(2.5)，其余項于是當時，〔3.6〕誤差階為，顯然是收斂的.與復(fù)合梯形公式相似有

實際上，只要就有

此外，由于中求積系數(shù)均為正數(shù)，故知復(fù)合辛普森公式計算穩(wěn)定.

46

例3對于函數(shù)，給出的函數(shù)表(見表4-2)，試用復(fù)合梯形公式(3.2)及復(fù)合辛普森公式(3.5)計算積分并估計誤差.

解將積分區(qū)間劃分為8等分，應(yīng)用復(fù)化梯形法求得

47同積分的準確值比較，復(fù)合梯形法的結(jié)果只有兩位有效數(shù)字，而復(fù)合辛普森法的結(jié)果卻有6位有效數(shù)字.而如果將分為4等分，應(yīng)用復(fù)化辛普森法有

接下來看誤差估計，由于所以有

以上得到的兩個結(jié)果與，都需要提供9個點上的函數(shù)值，計算量根本相同，然而精度卻差異很大.48于是由(3.3)得復(fù)合梯形公式誤差對復(fù)合辛普森公式，由(3.6)得（3.3）（3.6）49例4計算積分，假設(shè)用復(fù)合梯形公，問區(qū)間應(yīng)分多少等份才能使誤差不超過，假設(shè)改用復(fù)合辛普森公式，要到達同樣的精度，區(qū)間應(yīng)分多少等份？解此題只要根據(jù)及的余項公式即可求得其截斷誤差應(yīng)滿足的精度.

由于，由復(fù)合梯形公式的余項公式得誤差的上界為50因此有，可取，即將區(qū)間213等份，即可使誤差不超過假設(shè)采用復(fù)合辛普森公式計算積分，那么由余項公式，要滿足精度要求，必須使由此得可取，即用的復(fù)合辛普森公式計算即可到達精度要求，此時區(qū)間實際上應(yīng)分為8等份.51從這個例子可以看出，為到達同樣的精度，復(fù)合辛普森公式只需計算9個函數(shù)值，而復(fù)合梯形公式那么需214個函數(shù)值，工作量相差近24倍.524自適應(yīng)積分方法

復(fù)合求積方法是用于被積函數(shù)變化不太大的積分.如果在求積區(qū)間中被積函數(shù)變化很大，有的局部函數(shù)值變化劇烈，另一局部變化平緩，這時統(tǒng)一將區(qū)間等份用復(fù)合求積公式計算工作量就會很大.要到達誤差要求對變化劇烈局部必須將區(qū)間細分，而平緩局部那么可用大步長，即針對被積函數(shù)在區(qū)間上不同情形采用不同的步長，使得在滿足精度前提下積分計算的工作量盡可能小.53

針對這類問題的算法技巧是在不同區(qū)間上預(yù)測被積函數(shù)變化的劇烈程度確定相應(yīng)的步長.

這種方法稱為自適應(yīng)積分方法.以常用的復(fù)合辛普森公式為例說明方法的根本思想.54

設(shè)給定精度要求，計算積分的近似值.先取步長，應(yīng)用辛普森公式有其中假設(shè)把區(qū)間對分，步長，在每個小區(qū)間上用辛普森公式，那么得〔5.1〕〔5.2〕55實際上(5.2)即為與(5.1)比較，假設(shè)在上變化不大，可假定其中〔5.2〕’從而可得56假設(shè)不等式(5.3)不成立，那么應(yīng)分別對子區(qū)間及再用辛普森公式，此時步長,得到這里.如果有那么可期望得到與(5.2)比較，那么得此時可取作為的近似，那么可到達給定的誤差精度.〔5.3〕57

及.

只要分別考察及是否成立.對滿足要求的區(qū)間不再細分，對不滿足要求的還要繼續(xù)上述過程，直到滿足要求為止，最后還要應(yīng)用龍貝格法那么求出相應(yīng)區(qū)間積分的近似值.58

例7計算積分假設(shè)用復(fù)合辛普森法(3.5)，計算結(jié)果見表4-6.〔此處即為公式中的，積分精確值為4〕計算到為止，此時的近似值，假設(shè)再用龍貝格法那么得到整個計算是將做32等分，即需要計算33個的值.59現(xiàn)在假設(shè)用自適應(yīng)積分法，當時有由于

大于允許誤差，故要對及兩區(qū)間再用做積分.

先計算的積分

由于60小于允許誤差0.01，故在區(qū)間的積分值為下面再計算子區(qū)間的積分，其中由于而對可求得61大于允許誤差0.01，因此還要分別計算及的積分.

當時可求得而小于允許誤差0.005，故可得的積分近似而對區(qū)間，其誤差不小于0.005，故還要分別計算及的積分，62且小于允許誤差0.0025，故有最后子區(qū)間的積分可檢驗出它的誤差小于0.0025，且可得其中，當可求得63將以上各區(qū)間的積分近似值相加可得它一共只需計算17個的值.646.1

一般理論

求積公式

含有個待定參數(shù)

當為等距節(jié)點時得到的插值求積公式其代數(shù)精度至少為次.

如果適當選取有可能使求積公式具有次代數(shù)精度.6

高斯求積公式65試確定節(jié)點及和系數(shù)，使其具有近可能高的代數(shù)精度.（6.2）

例8求積公式（6.1）

解令公式(6.1)對于準確成立，得66

用(6.2)式中的第3式減去乘(6.2)中的第2式有用第4式減去第2式乘，得由此得于是可取.用前一式代入那么得由此得出與異號，即，從而有67

再由(6.2)式的第1式得，于是有（6.3）當時，(6.3)式兩端分別為及，(6.3)式對不精確成立，故公式(6.3)的代數(shù)精度為3.實際上，形如(6.1)的求積公式其代數(shù)精度不可能超過3，因為當時，設(shè)這是4次多項式，代入(6.1)式左端有，而右端為0.說明兩個節(jié)點的求積公式的代數(shù)精度為3.

一般節(jié)點的求積公式的代數(shù)精度最高為次.68為求積節(jié)點，可適中選取及使〔6.4〕式具有次代數(shù)精度.

下面研究帶權(quán)積分這里為權(quán)函數(shù)，類似(1.3)，求積公式為（6.4）為不依賴于的求積系數(shù).（1.3）69

根據(jù)定義要使(6.4)具有次代數(shù)精度，只要?。?.5）當給定權(quán)函數(shù)，求出右端積分，則可由(6.5)解得令(6.4)精確成立，即定義4如果求積公式(6.4)具有次代數(shù)精度，那么稱其節(jié)點為高斯點，相應(yīng)公式(6.4)稱為高斯求積公式.70(6.4)是關(guān)于及的非線性方程組，當時求解非常困難.如果事先確定了節(jié)點，那么可以利用(6.5)求解.此時(6.5)是關(guān)于的線性方程組.

下面討論如何選取節(jié)點才能使求積公式(6.4)具有次代數(shù)精度.71

設(shè)上的個節(jié)點的拉格朗日插值多項式為其中那么72用乘上式并從到積分，那么得其中余項73顯然當取為時有，此時有即求積公式至少具有次代數(shù)精度.

現(xiàn)在考察如何選取節(jié)點才能使求積公式精度提高到次.

此時要求為次多項式時，而當時，為次多項式.74假設(shè)要求對，積分即相當于要求與每個帶權(quán)在上正交.

也就是以節(jié)點為零點的次多項式是上帶權(quán)正交的多項式，故有以下定理.75

定理5插值型求積公式(6.4)的節(jié)點是高斯點的充分必要條件是以這些節(jié)點為零點的多項式與任何次數(shù)不超過的多項式帶權(quán)正交，（6.7）即76

定理表明在上帶權(quán)的次正交多項式的零點就是求積公式(6.4)的高斯點.

有了求積節(jié)點，再利用對成立，解此方程則得的線性方程.則得到一組關(guān)于求積系數(shù)

也可直接由的插值多項式求出求積系數(shù)77

例9確定求積公式

解具有最高代數(shù)精度的求積公式是高斯型求積公式，其節(jié)點為關(guān)于權(quán)函數(shù)的正交多項式零點及,的系數(shù)及節(jié)點和，使它具有最高的代數(shù)精度.

設(shè)由正交性知與1及帶權(quán)正交，即得于是得78由此解得即令，那么得

由于兩個節(jié)點的高斯求積公式具有3次代數(shù)精度，故公式對精確成立，即

當時

當時由此解出79

下面討論高斯求積公式(6.4)的余項.

利用在節(jié)點的埃爾米特插值于是即80兩端乘，并由到積分，則得（6.9）其中右端第一項積分對次多項式精確成立，故由于（6.10）由積分中值定理得(6.4)的余項為

關(guān)于高斯求積公式的穩(wěn)定性與收斂性，有：（6.4）81

定理6高斯求積公式(6.4)的求積系數(shù)全是正的.

證明考察它是次多項式，因而是次多項式，注意到故高斯求積公式(6.4)對于它能準確成立，即有上式右端實際上即等于從而有（6.4）82由本定理及定理2，那么得

推論高斯求積公式(6.4)是穩(wěn)定的.定理7設(shè)那么高斯求積公式(6.4)收斂，即定理得證.（6.4）83利用余項公式(8.4)知，帶余項的兩點公式是（8.4）84用MATLAB作數(shù)值積分矩形公式Sum(x)輸入數(shù)組x(即fk),輸出x的和(數(shù))cumsum(x)輸入數(shù)組x,輸出x的依次累加和(數(shù)組)梯形公式trapz(x)輸入數(shù)組x,輸出按梯形公式x的積分(單位步長)trapz(x,y)輸入同長度數(shù)組x,y,輸出按梯形公式y(tǒng)對x的積分(步長不一定相等)85用MATLAB作數(shù)值積分辛普森公式quad(@fun,a,b,tol,trace)[I,fn]=quad(…)用自適應(yīng)辛普森公式計算tol為絕對誤差，缺省時為10-6Gauss-Lobatto公式quadl(@fun,a,b,tol,trace)[I,fn]=quadl(…)用自適應(yīng)Gauss-Lobatto公式計算tol為絕對誤差，缺省時為10-6注意：fun.m中應(yīng)以自變量為矩陣的形式輸入(點運算)86矩形域上計算二重積分的命令：dblquad(@fun,xmin,xmax,ymin,ymax,tol)廣義積分、二重和三重積分長方體上計算三重積分的命令：triplequad(@fun,xmin,xmax,ymin,ymax,zmin,zmax,tol)注：fun是被積函數(shù)，本身可以有自己的參數(shù)廣義積分：通過分析和控制誤差，轉(zhuǎn)換成普通積分quadv(@fun,a,b,tol,trace)向量值積分：87用MATLAB作數(shù)值積分例.計算1〕矩形公式和梯