數(shù)值積分課件 (《計算方法》)

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

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

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

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

其原函數(shù)都不能用初等函數(shù)表示成有限形式。2/6/20233(2)被積函數(shù)f(x)沒有具體的解析表達式。其函數(shù)關(guān)系由表格或圖形表示,無法求出原函數(shù)。其被積函數(shù)的原函數(shù)就比較復(fù)雜,從數(shù)值計算角度來看,計算量太大。(3)盡管f(x)的原函數(shù)能表示成有限形式但其表達式相當(dāng)復(fù)雜。例如定積分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,

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

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

求積節(jié)點

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

確定求積系數(shù)Ak和求積節(jié)點xk；(ii)

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

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

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

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

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

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

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

n

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

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

建立數(shù)值積分公式最基本的思想是選取一個既簡單又有足夠精度的函數(shù)φ(x),用φ(x)代替被積函數(shù)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有關(guān)注意:由(7-11)式確定的Cotes系數(shù)只與i和n有關(guān),與f(x)和積分區(qū)間[a,b]無關(guān),且滿足:(7-9)2/6/202323稱Ci(n)為柯特斯求積系數(shù)。很顯然,當(dāng)n=1時,可算得此時式(7―10)為(7―12)這是梯形公式。2/6/202324

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

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

(n)

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

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

現(xiàn)在討論舍入誤差對計算結(jié)果產(chǎn)生的影響.設(shè)用公式

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

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

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

1.2誤差估計

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

f(n+1)(ξ)≡0

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

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

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

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

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

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

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

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

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

兩邊取對數(shù)得整理后得到

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

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

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

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

x

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

復(fù)化梯形公式：在每個上用梯形公式：=

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

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

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

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

圖7.5

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

其中

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