基于Matlab數(shù)值積分公式問題

1、 數(shù)值分析學(xué) 號(hào)：130080402015 學(xué)生所在學(xué)院：測(cè)試與光電工程學(xué)院學(xué) 生 姓 名 ：張翀任 課 教 師 ：鄭華盛教師所在學(xué)院：數(shù)信學(xué)院 基于Matlab的數(shù)值積分公式問題 張翀，測(cè)試與光電工程學(xué)院 測(cè)試計(jì)量技術(shù)及儀器，130080402015 摘 要：在求一些函數(shù)的定積分時(shí)，由于原函數(shù)十分復(fù)雜難以求出或用初等函數(shù)表達(dá)，導(dǎo)致積分很難精確求出，只能設(shè)法求其近似值，因此能夠直接借助牛頓-萊布尼茲公式計(jì)算定積分的情形是不多的。數(shù)值積分就是解決此類問題的一種行之有效的方法。積分的數(shù)值計(jì)算是數(shù)值分析的一個(gè)重要分支；因此，探討近似計(jì)算的數(shù)值積分方法是有著明顯的實(shí)際意義的。 本文介紹了數(shù)值積分法的幾

2、種計(jì)算公式，如矩形求積公式、梯形求積公式和辛普森求積公式及相應(yīng)的MATLAB命令,并給出了用 MATLAB編程求數(shù)值積分的實(shí)例。關(guān)鍵詞: MATLAB；數(shù)值積分；矩形求積公式；梯形求積公式；辛普森求積公式 目錄1引言.12數(shù)值積分算法介紹.1 2.1數(shù)值求積公式的構(gòu)造.1 2.2求積公式的推導(dǎo).2 2.3常見的牛頓-科特斯求積公式.5 2.4復(fù)合求積公式.73關(guān)于河流橫斷面積的數(shù)值積分問題.84問題的求解過程.95基于MATLAB編程的各種求積公式對(duì)問題的求解.96總結(jié).13參考文獻(xiàn) .14附錄.151 引言 實(shí)際問題當(dāng)中常常需要計(jì)算積分。有些數(shù)值方法，如微分方程和積分方程的求解，也都和積分計(jì)

3、算相聯(lián)系。在一元微積分學(xué)中,對(duì)于積分 ，只要找到被積函數(shù)f（x）原函數(shù)為F( x) ,求f(x)在該區(qū)間上的定積分便可用牛頓 - 萊布尼茲公式求解,即 。用牛頓 - 萊布尼茲公式計(jì)算定積分的方法在理論上和解決實(shí)際問題中起到了很大的作用 ,但它并不能解決定積分計(jì)算的所有問題。在工程技術(shù)領(lǐng)域常遇到十分復(fù)雜的情況而無法用牛頓 - 萊布尼茲公式求解.其可能出現(xiàn)的情況有:(1) 某些被積函數(shù)f(x),其原函數(shù)無法用初等函數(shù)表示 ,如 ， 等。(2) 函數(shù)f(x)結(jié)構(gòu)復(fù)雜,求其原函數(shù)非常困難。(3) 函數(shù)f(x)的結(jié)構(gòu)雖然簡(jiǎn)單且其原函數(shù)存在,但其原函數(shù)的結(jié)構(gòu)相對(duì)復(fù)雜。(4) 函數(shù)f(x)沒有具體的表達(dá)式,

4、只有一些由試驗(yàn)測(cè)試數(shù)據(jù)形成的表格或圖形。而在這些情況下 ,可采用 “數(shù)值積分”的方法求出定積分(近似值) 。2 數(shù)值積分算法介紹2.1數(shù)值求積公式的構(gòu)造大多數(shù)實(shí)際問題的積分是需要用數(shù)值積分方法求出近似結(jié)果的。數(shù)值積分原則上可以用于計(jì)算各種被積函數(shù)的定積分，無論被積函數(shù)是解析形式還是數(shù)表形式，其基本原理都是用多項(xiàng)式函數(shù)近似代替被積函數(shù)，用多項(xiàng)式的積分結(jié)果近似代替對(duì)被積函數(shù)的積分。由于所選多項(xiàng)式形式的不同，可以有許多種數(shù)值積分方法。而利用插值多項(xiàng)式來構(gòu)造數(shù)值求積公式是最常用的一種方法。對(duì)于積分，用一個(gè)容易積分的函數(shù)去代替被積函數(shù)，這樣的自然以多項(xiàng)式為最佳，因?yàn)槎囗?xiàng)式能很好的逼近任何連續(xù)函數(shù)，而且容

5、易求出其原函數(shù)。2.2求積公式的推導(dǎo)在積分區(qū)間上取有限個(gè)點(diǎn)，作的次插值多項(xiàng)式，其中，為次插值基函數(shù)。用近似代替被積函數(shù)，則得 （2.1）若記 （2.2）則得數(shù)值求積公式 （2.3）其中稱為求積系數(shù)，稱為求積節(jié)點(diǎn)。則稱該求積公式為插值型求積公式。為了便于計(jì)算與應(yīng)用，常將積分區(qū)間的等分點(diǎn)作為求積節(jié)點(diǎn)，這樣構(gòu)造出來的插值型求積公式就稱為牛頓-科特斯求積公式。在積分區(qū)間上取個(gè)等距節(jié)點(diǎn)，其中，做次拉格朗日插值多項(xiàng)式，因?yàn)?，所?記 （2.4） （2.5）截去第二項(xiàng)得 顯然與無關(guān)，只與節(jié)點(diǎn)有關(guān)。令，則當(dāng)時(shí)，于是 （2.6）而 從而得記 （2.7）則 故求積公式（2.3）可寫成 （2.8）這就是牛頓-科特斯

6、求積公式，其中稱為科特斯系數(shù)。部分科特斯系數(shù)取值如下表2.1科特斯系數(shù)具有以下特點(diǎn)（1） （2） （3）當(dāng) ³ 8 時(shí)，出現(xiàn)負(fù)數(shù)，穩(wěn)定性得不到保證。而且當(dāng) 較大時(shí)，由于Runge現(xiàn)象，收斂性也無法保證。故一般不采用高階的牛頓-科特斯求積公式。 （4）當(dāng) £ 7 時(shí)，牛頓-科特斯公式是穩(wěn)定的。表2.1 部分科特斯系數(shù)表知道了什么是牛頓-科特斯求積公式，下面我們來看它的誤差估計(jì)，首先來看看牛頓-科特斯求積公式的截?cái)嗾`差。我們知道牛頓-科特斯求積公式是一個(gè)插值型數(shù)值求積公式，當(dāng)用插值多項(xiàng)式代替進(jìn)行積分時(shí)，其截?cái)嗾`差即積分真值和近似值之差，推導(dǎo)如下，由插值多項(xiàng)式的誤差估計(jì)可知，用次

7、拉格朗日多項(xiàng)式逼近函數(shù)時(shí)產(chǎn)生的誤差為 （2.9） 其中。對(duì)上式兩邊從到作定積分，便可得出它的截?cái)嗾`差 （2.10）2.3常見的牛頓-科特斯求積公式2.3.1矩形求積公式 在牛頓-科特斯求積公式中，如果取，用零次多項(xiàng)式（即常數(shù)）代替被積函數(shù)，即用矩形面積代替曲邊梯形的面積，則有 （2.11）稱式（2.11）為矩形求積公式根據(jù)牛頓-科特斯求積公式的誤差理論式，矩形求積公式的誤差估計(jì)為2.3.2梯形求積公式 在牛頓-科特斯求積公式中，如果取，用一次多項(xiàng)式代替被積函數(shù)，即用梯形面積代替曲邊梯形的面積，則有其中，,查表可得代入上式得出 （2.12）稱式（2.12）為梯形求積公式由于用一次多項(xiàng)式近似代替被

8、積函數(shù)，所以它的精度是1。也就是說，只有當(dāng)被積函數(shù)是一次多項(xiàng)式時(shí)，梯形求積公式才是準(zhǔn)確的。根據(jù)牛頓-科特斯求積公式的誤差理論式（2.10），梯形求積公式的誤差估計(jì)為是被積函數(shù)二階導(dǎo)數(shù)在點(diǎn)的取值，2.3.3辛普森求積公式在牛頓-科特斯求積公式中，如果取，用二次多項(xiàng)式代替被積函數(shù)，即曲邊用拋物線代替，則有其中，,查表可得，代入上式得出 （2.13）稱式（2.13）為辛普森求積公式，也稱拋物線求積公式。它的幾何意義是：用過3個(gè)點(diǎn)，的拋物線和，構(gòu)成的曲邊梯形面積，近似地代替了被積函數(shù)形成的曲邊和，構(gòu)成的曲線梯形面積。下面對(duì)辛普森求積公式的誤差進(jìn)行估計(jì)。由于辛普森求積公式是用二次多項(xiàng)式逼近被積函數(shù)推得的

9、，原則上它的代數(shù)精度為2.但因多項(xiàng)式次數(shù)是偶數(shù)，根據(jù)定理1.1可知，它的代數(shù)精度為3過，和3個(gè)點(diǎn)，構(gòu)造一個(gè)的三次Lagrange插值多項(xiàng)式,且使。根據(jù)Lagrange插值余項(xiàng)定理得 對(duì)上式兩邊從到進(jìn)行積分，即可得到 （2.14）根據(jù)定積分中值定理可知，在上總有一點(diǎn)滿足下述關(guān)系：通過變量代換,，很容易求得把這個(gè)結(jié)果代入式（2.14），便得出辛普森求積公式的誤差估計(jì)式 （2.15）2.4復(fù)合求積公式前面導(dǎo)出的誤差估計(jì)式表明，用牛頓-科特斯公式計(jì)算積分近似值時(shí)，步長(zhǎng)越小，截?cái)嗾`差越小。但縮小步長(zhǎng)等于增加節(jié)點(diǎn)，亦即提高插值多項(xiàng)式的次數(shù)。龍格現(xiàn)象表明，這樣做并不一定能提高精度。理論上已經(jīng)證明，當(dāng)時(shí)，牛頓

10、-科特斯公式所求得的近似值不一定收斂于積分的準(zhǔn)確值，而且隨著的增大，牛頓-科特斯公式是不穩(wěn)定的。因此，實(shí)際中不常用高階牛頓-科特斯公式。為了提高計(jì)算精度，可考慮對(duì)被積函數(shù)用分段低次多項(xiàng)式插值，由此導(dǎo)出復(fù)合求積公式。用數(shù)值積分的方法求一個(gè)函數(shù)在區(qū)間a,b上的定積分,可利用定積分的定義來求解:，設(shè)，則。此時(shí)稱In為數(shù)值積分。顯然數(shù)值積分In就是I的近似值,并且當(dāng)n越大,In就越接近于精確值I.由于k取值不同,數(shù)值積分In的結(jié)果會(huì)有所不同。數(shù)值積分的計(jì)算公式也有多種:(1) 復(fù)合矩形公式將積分區(qū)間 a ,b n 等分 ,每個(gè)小區(qū)間寬度均為h = (b - a) / n ,h 稱為積分步長(zhǎng)。記a =

11、x0 < x1 < < xk < xn = b ,在小區(qū)間上用小矩形面積近似小曲邊梯形的面積 ,若分別取左端點(diǎn)和右端點(diǎn)的函數(shù)值為小矩形的高,則分別得到兩個(gè)曲邊梯形面積的近似計(jì)算公式: （2.16） （2.17）稱公式（2.16）,（2.17）分別為左、右矩形公式,兩個(gè)矩形面積分別小于和大于所求曲邊梯形的面積。(2) 復(fù)合梯形公式如果將二者求平均值 ,則每個(gè)小區(qū)間上的小矩形變?yōu)樾√菪?,整個(gè)區(qū)間上的值變?yōu)? （2.18）其中。（2.18）式稱為梯形公式.(3) 辛普森公式為了提高計(jì)算結(jié)果的精度 ,用分段二次插值函數(shù)代替f (x) .由于每段要用到相鄰兩個(gè)小區(qū)間端點(diǎn)的三個(gè)函

12、數(shù)值 ,所以 ,小區(qū)間的數(shù)目必須是偶數(shù) ,記n = 2m ,k = 0 ,1 , ,m - 1 ,在第k個(gè)小區(qū)間上用三個(gè)節(jié)點(diǎn)(x2k ,f2k) , (x2k+1 ,f2k+1) , (x2k+2 ,f 2k+2) 作二次插值函數(shù) Sk (x) ,然后積分可得: （2.19）（2.19）式稱為辛普森公式。3 關(guān)于河流橫斷面積的數(shù)值積分問題設(shè)河面寬20 m，從河的一岸向另一岸每隔2 m測(cè)得的水深如下（單位：m） 表3.1 河流水深數(shù)據(jù)表河寬 02468水深 00.61.42.02.3河寬 101214161820水深 2.12.51.91.20.70其中試求河流的橫斷面積。4 問題的求解過程首先

13、可通過梯形積分公式直接對(duì)離散點(diǎn)進(jìn)行求解，但是精度不一定好。通過對(duì)離散點(diǎn)進(jìn)行觀察可發(fā)現(xiàn)用正弦曲線逼近擬合效果較好，這樣可得到擬合曲線，通過對(duì)擬合曲線進(jìn)行復(fù)合矩形公式和復(fù)合梯形公式以及辛普森公式進(jìn)行求解。通過對(duì)各種方法得到的結(jié)果進(jìn)行對(duì)比可以得到最優(yōu)的結(jié)果。5 基于MATLAB編程的各種求積公式對(duì)問題的求解首先對(duì)離散點(diǎn)進(jìn)行擬合>> x=0 2 4 6 8 10 12 14 16 18 20;>> y=0 0.6 1.4 2.0 2.3 2.1 2.5 1.9 1.2 0.7 0;>> cftool從離散點(diǎn)圖可以看出應(yīng)該使用正弦曲線逼近。然后翻轉(zhuǎn)圖像 General

14、 model Sin3: f(x) = a1*sin(b1*x+c1) + a2*sin(b2*x+c2) + a3*sin(b3*x+c3)Coefficients (with 95% confidence bounds): a1 = 2.368 (2.326, 2.411) b1 = 0.159 (0.1566, 0.1614) c1 = -0.00736 (-0.0367, 0.02198) a2 = 0.1294 (0.08601, 0.1728) b2 = 1.341 (1.299, 1.383) c2 = -8.72 (-9.197, -8.243) a3 = -0.1724 (-

15、0.2099, -0.1348) b3 = 0.9055 (0.8626, 0.9483) c3 = -0.5845 (-1.044, -0.1247)Goodness of fit: SSE: 0.0007643 R-square: 0.9999 Adjusted R-square: 0.9995 RMSE: 0.01955（1）矩形公式求積： 1）左矩形公式求積 當(dāng)分為500個(gè)矩形時(shí)：>> x=linspace(0,20,500); y= 2.368*sin(0.159*x-0.00736) + 0.1294*sin(1.341*x-8.72) -0.1724 *sin(0.9

16、055*x-0.5845); I1=sum( y(1:499)*(20)/(500-1)I1 = 29.5219當(dāng)分為5000個(gè)矩形時(shí)：>> x=linspace(0,20,5000); y= 2.368*sin(0.159*x-0.00736) + 0.1294*sin(1.341*x-8.72) -0.1724 *sin(0.9055*x-0.5845); I1=sum( y(1:4999)*(20)/(5000-1)I1 = 29.5221當(dāng)分為50000個(gè)矩形時(shí)：>> x=linspace(0,20,50000); y= 2.368*sin(0.159*x-0.

17、00736) + 0.1294*sin(1.341*x-8.72) -0.1724 *sin(0.9055*x-0.5845); I1=sum( y(1:49999)*(20)/(50000-1)I1 = 29.5222當(dāng)分為500000個(gè)矩形時(shí)：>> x=linspace(0,20,500000); y= 2.368*sin(0.159*x-0.00736) + 0.1294*sin(1.341*x-8.72) -0.1724 *sin(0.9055*x-0.5845); I1=sum( y(1:499999)*(20)/(500000-1)I1 = 29.5222所以面積為29

18、.5222 ()。 2）右矩形公式求積當(dāng)分為500個(gè)矩形時(shí)：>> x=linspace(0,20,500); y= 2.368*sin(0.159*x-0.00736) + 0.1294*sin(1.341*x-8.72) -0.1724 *sin(0.9055*x-0.5845); I1=sum( y(2:500)*(20)/(500-1) I1 = 29.5223當(dāng)分為5000個(gè)矩形時(shí)：>> x=linspace(0,20,5000); y= 2.368*sin(0.159*x-0.00736) + 0.1294*sin(1.341*x-8.72) -0.1724

19、*sin(0.9055*x-0.5845);I1=sum( y(2:5000)*(20)/(5000-1)I1 = 29.5222當(dāng)分為50000個(gè)矩形時(shí)：>> x=linspace(0,20,50000); y= 2.368*sin(0.159*x-0.00736) + 0.1294*sin(1.341*x-8.72) -0.1724 *sin(0.9055*x-0.5845); I1=sum( y(2:50000)*(20)/(50000-1)I1 = 29.5222所以面積為29.5222 ()。（2）梯形公式求積： 1）直接去離散點(diǎn)進(jìn)行求積：>> x=0,2,4

20、,6,8,10,12,14,16,18,20;y=0,0.6,1.4,2.0,2.3,2.1,2.5,1.9,1.2,0.7,0;trapz(x,y)ans = 29.4000 2）運(yùn)用擬合曲線進(jìn)行求積：當(dāng)分為500個(gè)梯形時(shí)：>> x=linspace(0,20,500);y=2.368*sin(0.159*x-0.00736) + 0.1294*sin(1.341*x-8.72) -0.1724 *sin(0.9055*x-0.5845);I3=trapz(x,y)I3 = 29.5221當(dāng)分為5000個(gè)梯形時(shí)：>> x=linspace(0,20,5000);y=2

21、.368*sin(0.159*x-0.00736) + 0.1294*sin(1.341*x-8.72) -0.1724 *sin(0.9055*x-0.5845);I3=trapz(x,y)I3 = 29.5222當(dāng)分為50000個(gè)梯形時(shí)：>> x=linspace(0,20,50000);y=2.368*sin(0.159*x-0.00736) + 0.1294*sin(1.341*x-8.72) -0.1724 *sin(0.9055*x-0.5845);I3=trapz(x,y)I3 = 29.5222所以面積為29.5222 ()。（3）用辛普森公式求積：I4=quad(

22、'2.368*sin(0.159*x-0.00736)+0.1294*sin(1.341*x-8.72)-0.1724*sin(0.9055*x-0.5845)',0 ,20)I4 = 29.5222所以面積為29.5222 ()。比較上述方法的結(jié)果，認(rèn)定河流的橫切面積為29.5222 ()。6 總結(jié)本文主要討論了對(duì)實(shí)際問題的數(shù)值積分計(jì)算方法，并通過MATLAB軟件編程實(shí)現(xiàn)，通過前面的研究我們知道求數(shù)值積分近似值的計(jì)算方法很多，有矩形求積公式、梯形求積公式、辛普森求積公式和相應(yīng)的復(fù)合公式等等。 牛頓-科特斯方法是一種利用插值多項(xiàng)式來構(gòu)造數(shù)值積分的常用方法，這其中梯形積分方法的誤

23、差最大，近似效果最差，辛普森方法的精度比梯形積分高了一個(gè)數(shù)量級(jí)，它的代數(shù)精度比梯形積分的代數(shù)精度高，能更好地近似積分值。復(fù)合梯形積分方法比單獨(dú)的梯形積分精度高，它的積分精度和被積函數(shù)有關(guān)，還和復(fù)合積分時(shí)的步長(zhǎng)有關(guān)。 一般來說，牛頓-科特斯方法的代數(shù)精度越高，數(shù)值積分的效果越好、越精確。當(dāng)積分區(qū)間比較大的時(shí)候，可以采用復(fù)合積分方法可以得到更好的效果。 參考文獻(xiàn)1 樊守芳Newton-Cotes數(shù)值求積公式的注記J棗莊學(xué)院學(xué)報(bào)，2011，28(2)：186-1902 劉小偉基于MATLAB的復(fù)合梯形數(shù)值積分法的研究與實(shí)驗(yàn)J甘肅聯(lián)合大學(xué)學(xué)報(bào)(自然科學(xué)版)，2010，24(4)：20-233 3陳佩寧

24、,劉競(jìng)用. MATLAB求數(shù)值積分的方法J，石家莊職業(yè)技術(shù)學(xué)院學(xué)報(bào)2008，12(6)：58-604 王建強(qiáng) 多種數(shù)值積分方法比較分析J 武漢大學(xué)測(cè)繪軒轅學(xué)報(bào)，2010，2(1)：104-1065 伍麗華，周玲麗數(shù)學(xué)軟件教程M廣州：中山大學(xué)出版社，20086 張德豐等MATLAB數(shù)值計(jì)算方法M北京：機(jī)械工業(yè)出版社，20107 Li F，Li XThe neighbor-scattering number can be computed in polynomial time for interval graphsJComputers and Mathematics with Applicatio

25、ns，2007，54：679-6868 GUO RuA Note on Newton-CoteS Numerical Integral FormulaJ. Journal of Xinjiang University(Natural Science Edition)，2010，27(2)：186-190附錄擬合M文件：function nihe(x,y)%NIHE Create plot of datasets and fits% NIHE(X,Y)% Creates a plot, similar to the plot in the main curve fitting% window,

26、using the data that you provide as input. You can% apply this function to the same data you used with cftool% or with different data. You may want to edit the function to% customize the code and this help message.% Number of datasets: 2% Number of fits: 1 % Data from dataset "y vs. x":% X

27、= x:% Y = y:% Unweighted % Data from dataset "y vs. x (2 )":% X = x:% Y = y:% Unweighted% This function was automatically generated on 16-Dec-2013 21:12:43 % Set up figure to receive datasets and fitsf_ = clf;figure(f_);set(f_,'Units','Pixels','Position',719 71 688

28、488);legh_ = ; legt_ = ; % handles and text for legendxlim_ = Inf -Inf; % limits of x axisax_ = axes;set(ax_,'Units','normalized','OuterPosition',0 0 1 1);set(ax_,'Box','on');axes(ax_); hold on; % - Plot data originally in dataset "y vs. x"x = x(:);y

29、 = y(:);% This dataset does not appear on the plot% Add it to the plot by removing the if/end statements that follow% and by selecting the desired color and markerif 0 h_ = line(x,y,'Color','r','Marker','.','LineStyle','none'); xlim_(1) = min(xlim_(1),

30、min(x); xlim_(2) = max(xlim_(2),max(x); legh_(end+1) = h_; legt_end+1 = 'y vs. x'end % end of "if 0" % - Plot data originally in dataset "y vs. x (2 )"h_ = line(x,y,'Parent',ax_,'Color',0.333333 0.666667 0,. 'LineStyle','none', 'LineWid

31、th',1,. 'Marker','.', 'MarkerSize',12);xlim_(1) = min(xlim_(1),min(x);xlim_(2) = max(xlim_(2),max(x);legh_(end+1) = h_;legt_end+1 = 'y vs. x (2 )' % Nudge axis limits beyond data limitsif all(isfinite(xlim_) xlim_ = xlim_ + -1 1 * 0.01 * diff(xlim_); set(ax_,'XLim',xlim_)end % - Create fit "fit 1"fo_ = fitoptions('method','NonlinearLeastSquares','Lower',-Inf 0 -Inf -Inf 0 -Inf -Inf 0 -Inf );ok_ = (isnan(x) | isnan(y);st_ = 0.6945672404255 0.6213101307954 0.7948210802