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

第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)雖然在理論上或在解決實(shí)際問題中都起了很大的作用,但它并不能完全解決定積分的計(jì)算問題。因?yàn)槎ǚe分的計(jì)算常常會(huì)碰到以下三種情況:

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

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

(7―2)2/6/20235同樣可得到右矩形公式：(7―3)2/6/20236圖7.2如圖7.2,若用梯形的面積近似地代替曲邊梯形的面積,則得到計(jì)算定積分的梯形公式(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,

三點(diǎn)高度的加權(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個(gè)點(diǎn){xi},(i=0,1,2,…n)處的高度{f(xi)}(i=0,1,…,n)通過加權(quán)平均的方法近似地得出平均高度f(ξ),這類求積方法稱為機(jī)械求積:2/6/202310或?qū)懗?數(shù)值積分公式求積系數(shù)

求積節(jié)點(diǎn)

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

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

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

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

(i)對(duì)所有次數(shù)≤m次的多項(xiàng)式,有

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

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

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

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

n

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

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

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

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

2/6/202320為牛頓―柯特斯(Newton-Cotes)求積公式,Rn(f)為牛頓―柯特斯求積公式的余項(xiàng)。我們稱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時(shí),可算得此時(shí)式(7―10)為(7―12)這是梯形公式。2/6/202324

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

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

(n)

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

只需驗(yàn)證當(dāng)n為偶數(shù)時(shí),Newton-Cotes公式對(duì)f(x)=xn+1的余項(xiàng)為零.由于f(x)=xn+1,所以f(n+1)(x)=(n+1)!.由式(7-9)得引進(jìn)變換t=u+n/2,因?yàn)閚為偶數(shù),故n/2為整數(shù),于是有據(jù)此可斷定R(f)=0,因?yàn)樯鲜霰环e函數(shù)是個(gè)奇函數(shù).2/6/202329Newton-Cotes公式的數(shù)值穩(wěn)定性

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

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

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

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

1.2誤差估計(jì)

易知,牛頓―柯特斯求積公式(7―10)對(duì)任何不高于n次的多項(xiàng)式是準(zhǔn)確成立的。這是因?yàn)?/p>

f(n+1)(ξ)≡0

故Rn(f)≡02/6/202333牛頓―柯特斯求積公式的代數(shù)精確度至少為n。通常在基點(diǎn)個(gè)數(shù)相等的情況下,代數(shù)精確度愈高,求積公式就愈精確。一般說來,若某個(gè)求積公式對(duì)于次數(shù)不高于m的多項(xiàng)式都準(zhǔn)確成立(即Rn(f)≡0),而對(duì)于某一次數(shù)為m+1的多項(xiàng)式并不準(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)知,梯形公式的余項(xiàng)為ω1(x)在區(qū)間(a,b)內(nèi)不變號(hào),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公式會(huì)出現(xiàn)數(shù)值不穩(wěn)定,低階Newton-Cotes公式有時(shí)又不能滿足精度要求.解決這個(gè)矛盾的辦法是將積分區(qū)間[a,b]分成若干小區(qū)間,在每個(gè)小區(qū)間上用低階求積公式計(jì)算,然后將它們加起來,這就是復(fù)合求積方法.2/6/202338§2復(fù)合求積公式

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

i+1］,這里步長(zhǎng)在每一個(gè)子區(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)其余項(xiàng)為2/6/202342例2若用復(fù)合梯形公式計(jì)算積分

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

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

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

(7―25)2/6/202348圖7.4復(fù)合拋物線公式框圖2/6/202349的數(shù)據(jù)表,分別用復(fù)合梯形公式和復(fù)合拋物線公式計(jì)算例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ù)合拋物線公式可得比較上面兩個(gè)結(jié)果T8和S4,它們都需要提供9個(gè)點(diǎn)上的函數(shù)值工作量基本相同,然而精度卻差別很大.同積分的準(zhǔn)確值I(f)=0.9460831比較,復(fù)化梯形法的結(jié)果T8=0.9456909只有兩位有效數(shù)字,而復(fù)化Simpson法的結(jié)果S4=0.9460832卻有六位有效數(shù)字.2/6/202352

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

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

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

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

(7―26)其中再把每個(gè)子區(qū)間分成兩半,用逐次將區(qū)間［a,b］分成2,4,…,2m等分,并按復(fù)合拋物線公式逐次計(jì)算積分得到S1,S2,…,Sm,而2/6/202356作步長(zhǎng),按復(fù)合拋物線公式計(jì)算出積分的近似值S2m。對(duì)于相鄰兩次的積分近似值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變步長(zhǎng)復(fù)合拋物線公式2/6/202358

圖7.5

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

其中

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