第7章數值積分§1插值型求積公式§2復化求積公式§3龍貝格(Romberg)求積方法2/6/20231

§1插值型求積公式在一元函數的積分學中,我們已經熟知,若函數f(x)在區間［a,b]上連續且其原函數為F(x),則可用牛頓―萊布尼茲公式(7―1)來求定積分。2/6/20232

公式(7―1)雖然在理論上或在解決實際問題中都起了很大的作用,但它并不能完全解決定積分的計算問題。因為定積分的計算常常會碰到以下三種情況:

(1)被積函數f(x)的原函數F(x)不易找到。許多很簡單的函數,例如

其原函數都不能用初等函數表示成有限形式。2/6/20233(2)被積函數f(x)沒有具體的解析表達式。其函數關系由表格或圖形表示,無法求出原函數。其被積函數的原函數就比較復雜,從數值計算角度來看,計算量太大。(3)盡管f(x)的原函數能表示成有限形式但其表達式相當復雜。例如定積分2/6/20234圖7.1如圖7.1,若用左矩形近似地代替曲邊梯形,則得到左矩形公式

(7―2)2/6/20235同樣可得到右矩形公式：(7―3)2/6/20236圖7.2如圖7.2,若用梯形的面積近似地代替曲邊梯形的面積,則得到計算定積分的梯形公式(7―4)2/6/20237如圖7.3,若用拋物線代替曲線f(x),則可得到拋物線公式(或辛普生公式)(7―5)圖7.32/6/20238此外,眾所周知的梯形公式:

I(f)≈(b-a)[f(a)+f(b)]/2和Simpson公式:

I(f)≈(b-a)[f(a)+4f((a+b)/2)+f(b)]/6則分別可以看作用a,b,c=(a+b)/2,

三點高度的加權平均值[f(a)+f(b)]/2和[f(a)+4f(c)+f(b)]/6作為平均高度f(ξ)的近似值.2/6/20239

更一般地,取區間[a,b]內n+1個點{xi},(i=0,1,2,…n)處的高度{f(xi)}(i=0,1,…,n)通過加權平均的方法近似地得出平均高度f(ξ),這類求積方法稱為機械求積:2/6/202310或寫成:數值積分公式求積系數

求積節點

(1)2/6/202311記稱(2)為數值求積公式,(3)為求積公式余項(誤差).構造或確定一個求積公式，要討論解決的問題有(i)

確定求積系數Ak和求積節點xk；(ii)

求積公式的誤差估計和收斂性為了構造形如式(2)的求積公式,需要提供一種判定求積方法精度高低準則2/6/202312求積公式的代數精度定義1

稱求積公式(2)具有m次代數精度,如果它滿足如下兩個條件:

(i)對所有次數≤m次的多項式,有

(ii)存在m+1次多項式,使得定義1中的條件(i),(ii)等價于:2/6/202313插值型求積公式

在積分區間[a,b]上取n+1個節點xi，i=0,1,2,…,n,作f(x)的n次代數插值多項式（拉格朗日插值公式）:則有

為插值余項于是有2/6/202314取稱(4)式為插值型求積公式,其中求積系數Ak由(5)式確定.(4)(5)Ak由節點決定，與f(x)無關。2/6/2023152/6/202316推論1

求積系數滿足:誤差定理1形如的求積公式至少有

n

次代數精度該公式為插值型（即：）2/6/202317現用第六章介紹的插值多項式Pn(x)來代替被積函數f(x),即有1.1牛頓―柯特斯公式(Newton―Cotes)

取節點為等距,即a=x0＜x1＜…＜xn=b

建立數值積分公式最基本的思想是選取一個既簡單又有足夠精度的函數φ(x),用φ(x)代替被積函數f(x),于是有

2/6/202318利用拉格朗日插值多項式(7―6)其中(7―7)2/6/202319這里yi=f(xi),對式(7―6)兩邊積分得

2/6/202320為牛頓―柯特斯(Newton-Cotes)求積公式,Rn(f)為牛頓―柯特斯求積公式的余項。我們稱2/6/202321令x=x0+sh,0≤s≤ndx=hds=(b-a)/nds(7―11)2/6/202322Newton-Cotes公式的誤差為:與x有關注意:由(7-11)式確定的Cotes系數只與i和n有關,與f(x)和積分區間[a,b]無關,且滿足:(7-9)2/6/202323稱Ci(n)為柯特斯求積系數。很顯然,當n=1時,可算得此時式(7―10)為(7―12)這是梯形公式。2/6/202324

當n=2時,可得于是(7―13)這是拋物線(Simpson)公式。2/6/202325當n=3時,代入(7―10)式得到求積公式

2/6/202326類似地可分別求出n=4,5,…時的柯特斯系數,從而建立相應的求積公式。具體結果見表7―1。從表中可以看出,當n≤7時,柯特斯系數為正;從n≥8開始,柯特斯系數有正有負。因此,當n≥8時,誤差有可能傳播擴大,牛頓―柯特斯求積公式不宜采用。柯特斯系數Ci

(n)

僅與n和i有關,與被積函數f(x)無關,且滿足(7―15)事實上,式(7―10)對f(x)=1是準確成立的。2/6/202327表7―12/6/202328定理當階數n為偶數時,Newton-Cotes公式(8)至少具有n+1次代數精度.證明

只需驗證當n為偶數時,Newton-Cotes公式對f(x)=xn+1的余項為零.由于f(x)=xn+1,所以f(n+1)(x)=(n+1)!.由式(7-9)得引進變換t=u+n/2,因為n為偶數,故n/2為整數,于是有據此可斷定R(f)=0,因為上述被積函數是個奇函數.2/6/202329Newton-Cotes公式的數值穩定性

現在討論舍入誤差對計算結果產生的影響.設用公式

近似計算積分時,其中計算函數值f(xj)有誤差εj

(j=0,1,2,…,n).設計算Cj(n)沒有誤差,中間計算過程中的舍入誤差也不考慮,則在式(10)的計算中,由εj引起的誤差為(10)2/6/202330如果Cj(n)都是正數,并設故en是有界的,即由εj引起的誤差受到控制,不超過ε的(b-a)倍,保證了數值計算的穩定性.而當n>7時,Cj(n)將出現負數,保證數值穩定性.因此高階公式不宜采用,有實用價值的僅僅是幾種低階的求積公式.將隨n增大,因而不能則有2/6/202331解利用梯形公式利用拋物線公式例1試分別用梯形公式和拋物線公式計算積分:原積分的準確值

2/6/202332現對牛頓―柯特斯求積公式所產生的誤差作一個分析。由式(7―9),牛頓―柯特斯求積公式的余項為

1.2誤差估計

易知,牛頓―柯特斯求積公式(7―10)對任何不高于n次的多項式是準確成立的。這是因為

f(n+1)(ξ)≡0

故Rn(f)≡02/6/202333牛頓―柯特斯求積公式的代數精確度至少為n。通常在基點個數相等的情況下,代數精確度愈高,求積公式就愈精確。一般說來,若某個求積公式對于次數不高于m的多項式都準確成立(即Rn(f)≡0),而對于某一次數為m+1的多項式并不準確成立,則稱這一求積公式的代數精確度為m。定理1(梯形公式的誤差)設f(x)在區間［a,b］上具有連續的二階導數,則梯形求積公式的誤差為2/6/202334由于ω1(x)=(x-a)(x-b)證由式(7―9)知,梯形公式的余項為ω1(x)在區間(a,b)內不變號,f″(ξ)是x的函數且在［a,b］上連續,故根據積分第二中值定理知,存在某一η∈(a,b)使2/6/202335

定理2(拋物線公式的誤差)設f(x)在［a,b］上有連續的四階導數,則拋物線公式的誤差為(7―17)2/6/202336n=1:TrapezoidalRule/*令x=a+th,h=ba,用中值定理*/代數精度=1n=2:Simpson’sRule代數精度=3n=4:CotesRule,代數精度=5,2/6/202337復合求積公式高次插值有Runge現象,高階Newton-Cotes公式會出現數值不穩定,低階Newton-Cotes公式有時又不能滿足精度要求.解決這個矛盾的辦法是將積分區間[a,b]分成若干小區間,在每個小區間上用低階求積公式計算,然后將它們加起來,這就是復合求積方法.2/6/202338§2復合求積公式

2.1復合梯形公式對于定積分(7―1),將積分區間［a,b］分成n個相等的子區間［xi,x

i+1］,這里步長在每一個子區間［xi,x

i+1］上使用梯形公式,則2/6/202339

復化梯形公式積分法2/6/202340相加后得

(7―18)(7―19)若f″(x)在［a,b］上連續,由連續函數的介值定理,存在某一ξ∈(a,b)使得2/6/202341因而

于是得到復合梯形公式(7―21)其余項為2/6/202342例2若用復合梯形公式計算積分

則當0＜x＜1時,有因為又解由余項(7―21)式問積分區間要等分多少才能保證有五位有效數字?2/6/202343由于原積分的準確值具有一位整數,因此要使近似積分值有五位有效數字,只需取n滿足

兩邊取對數得整理后得到

取n=68.2/6/202344類似復合梯形公式的做法,把區間［a,b］分成n個相等的子區間［x2i,x2i+2］(i=0,1,…,n-1),設每個子區間上的中點為x2i+1(i=0,1,…,n-1),且(7―22)2.2復合拋物線公式在每一個子區間［x2i,x2i+2

］上利用拋物線公式得2/6/202345

復化Simpson公式積分法2/6/202346相加后得(7―23)2/6/202347若f(4)(x)在［a,b］上連續,則從而得到復合拋物線公式(7―24)其余項為

(7―25)2/6/202348圖7.4復合拋物線公式框圖2/6/202349的數據表,分別用復合梯形公式和復合拋物線公式計算例3已知函數

x

f(x)011/80.99739782/80.98961583/80.97672674/80.95885105/80.93615566/80.90885167/80.877192510.84147092/6/202350解用復合梯形公式,這里2/6/202351用復合拋物線公式可得比較上面兩個結果T8和S4,它們都需要提供9個點上的函數值工作量基本相同,然而精度卻差別很大.同積分的準確值I(f)=0.9460831比較,復化梯形法的結果T8=0.9456909只有兩位有效數字,而復化Simpson法的結果S4=0.9460832卻有六位有效數字.2/6/202352

復化梯形公式：在每個上用梯形公式：=

Tn/*中值定理*/2/6/202353

復化Simpson公式：2/6/2023542.3變步長公式前面介紹的復合梯形公式和復合拋物線公式的步長都是預先確定的。它的主要缺點是事先很難估計出n的大小(或步長h的大小),使結果達到預先給定的精度。在實際計算中,我們常常借助于計算機來完成積分步長h的自動選擇,即采用變步長求積公式。具體地講,就是將步長逐次折半,反復利用復合求積公式,直到滿足精度要求為止。

2/6/202355下面介紹變步長復合拋物線公式(變步長復合梯形公式留給讀者作為練習)。

(7―26)其中再把每個子區間分成兩半,用逐次將區間［a,b］分成2,4,…,2m等分,并按復合拋物線公式逐次計算積分得到S1,S2,…,Sm,而2/6/202356作步長,按復合拋物線公式計算出積分的近似值S2m。對于相鄰兩次的積分近似值Sm、S2m,考察當｜S2m｜＜1當｜S2m｜≥1(7―27)設預先給定的精度為ε,若｜d｜＜ε則以S2m作為所要求的積分近似值,否則繼續將區間分半,利用復合拋物線公式求積分,直到滿足預給的精度為止。2/6/202357圖7.5變步長復合拋物線公式2/6/202358

圖7.5

變步長復合拋物線公式2/6/202359§3龍貝格(Romberg)積分方法我們已經知道,當被積函數f(x)在區間［a,b］上連續時,要使得復合梯形公式或復合拋物線公式比較精確地代替定積分可將分點(即基點)加密,也就是將區間［a,b］細分,然后利用復合梯形公式或復合拋物線公式求積。2/6/202360若用Tm表示把［a,b］作m等分并按復合梯形公式求積的結果,將每一小段再對分,令新的小段的長h′=h/2,則T2m與Tm之間有如下關系:(5―28)

其中

2/6/202361另外,若用Sm表示把［a,b］分成m(偶數)個小段按復合拋物線公式計算的結果,那么只要把Sm中的m改為2m,h改為h′就有從Tm的定義可得到關系式(5―29)2/6/202362我們再舉一個計算上半單位圓面積的例子(它的準確面積為π/2)。現用內接正多邊形的逼近方法來計算。如圖5.6,圖(a)、圖(b)是用同樣的內接正多邊形計算上半單位圓的面積。圖(a)是用梯形方法計算其面積,圖(b)是用三角形方法計算其面積。2/6/202363圖5.62/6/202364設正多邊形邊數為n=2k,則由圖(b)利用三角形公式算得面積為同理2/6/202365如果組合一下,就會得到更精確的結果,即同理2/6/202366再以類似方法組合得這樣繼續下去,其值越來越接近上半單位圓面積π/2。這種方法可以用到計算定積分2/6/202367為了推廣公式(5―29)和上述計算上半單位圓面積的組合方法,我們引進龍貝格求積算法。龍貝格求積算法本來是利用所謂外推法構造出的一種計算積分的方法。為了避免從外推引入而帶來理論上的麻煩,我們將直接從構造一個T數表開始。首先將［a,b］依次作20,21,22,…等分,記2/6/202368按復合梯形公式(5―20)算得的值相應地記為T(