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

數(shù)值積分21

數(shù)值積分概論

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

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

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

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

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

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

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

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

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

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

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

的權(quán).

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

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

111.2

代數(shù)精度的概念

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

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

精度.

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

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

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

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

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

如

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

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

例1

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

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

當(dāng)時(shí)，得

當(dāng)時(shí)，得19

當(dāng)時(shí)，得解得,于是得

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

插值型的求積公式

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

.取

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

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

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

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

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

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

而當(dāng)時(shí)，(1.8)的右端故可求得1.4

求積公式的余項(xiàng)25〔1.9〕代入余項(xiàng)(1.8)中可以得到更細(xì)致的余項(xiàng)表達(dá)式.

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

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

例2

求例1中求積公式的余項(xiàng)

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

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

2.1

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

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

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

當(dāng)時(shí)，這時(shí)的求積公式就是梯形公式30

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

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

31

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

（2.4）這里

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

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

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

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

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

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

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

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

（2.1）372.3

辛普森公式的余項(xiàng)

對(duì)牛頓-柯特斯求積公式，通常只使用時(shí)的三個(gè)公式，時(shí)為梯形公式，余項(xiàng)為

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

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

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

復(fù)合梯形公式

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

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

（3.2）41

由(1.10)，余項(xiàng)由于,且

所以使

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

42（3.3）誤差是階，且當(dāng)時(shí)有

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

事實(shí)上只要設(shè)，就可以得到收斂性，因?yàn)橹灰獙⒏膶憺?/p>

43

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

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

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

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

實(shí)際上，只要就有

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

46

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

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

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

接下來看誤差估計(jì)，由于所以有

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

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

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

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

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

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

及.

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

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

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

先計(jì)算的積分

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

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

一般理論

求積公式

含有個(gè)待定參數(shù)

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

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

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

例8求積公式（6.1）

解令公式(6.1)對(duì)于準(zhǔn)確成立，得66

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

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

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

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

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

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

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

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

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

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

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

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

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

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

例9確定求積公式

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

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

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

當(dāng)時(shí)

當(dāng)時(shí)由此解出79

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

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

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

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

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

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