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

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

2、種計算公式，如矩形求積公式、梯形求積公式和辛普森求積公式及相應(yīng)的MATLAB命令,并給出了用 MATLAB編程求數(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編程的各種求積公式對問題的求解.96總結(jié).13參考文獻(xiàn) .14附錄.151 引言 實際問題當(dāng)中常常需要計算積分。有些數(shù)值方法，如微分方程和積分方程的求解，也都和積分計

3、算相聯(lián)系。在一元微積分學(xué)中,對于積分 ，只要找到被積函數(shù)f（x）原函數(shù)為F( x) ,求f(x)在該區(qū)間上的定積分便可用牛頓 - 萊布尼茲公式求解,即 。用牛頓 - 萊布尼茲公式計算定積分的方法在理論上和解決實際問題中起到了很大的作用 ,但它并不能解決定積分計算的所有問題。在工程技術(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)雖然簡單且其原函數(shù)存在,但其原函數(shù)的結(jié)構(gòu)相對復(fù)雜。(4) 函數(shù)f(x)沒有具體的表達(dá)式,

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

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

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

7、的誤差為 （2.9） 其中。對上式兩邊從到作定積分，便可得出它的截斷誤差 （2.10）2.3常見的牛頓-科特斯求積公式2.3.1矩形求積公式 在牛頓-科特斯求積公式中，如果取，用零次多項式（即常數(shù)）代替被積函數(shù)，即用矩形面積代替曲邊梯形的面積，則有 （2.11）稱式（2.11）為矩形求積公式根據(jù)牛頓-科特斯求積公式的誤差理論式，矩形求積公式的誤差估計為2.3.2梯形求積公式 在牛頓-科特斯求積公式中，如果取，用一次多項式代替被積函數(shù)，即用梯形面積代替曲邊梯形的面積，則有其中，,查表可得代入上式得出 （2.12）稱式（2.12）為梯形求積公式由于用一次多項式近似代替被積函數(shù)，所以它的精度是1。也

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

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

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

11、 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從離散點圖可以看出應(yīng)該使用正弦曲線逼近。然后翻轉(zhuǎn)圖像 General 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,

12、 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 (-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

13、個矩形時： 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(1:499)*(20)/(500-1)I1 = 29.5219當(dāng)分為5000個矩形時： 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)*(

14、20)/(5000-1)I1 = 29.5221當(dāng)分為50000個矩形時： 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(1:49999)*(20)/(50000-1)I1 = 29.5222當(dāng)分為500000個矩形時： x=linspace(0,20,500000); y= 2.368*sin(0.159*x-0.00736) + 0.1294*sin(1.341*x-8.72) -0.1724

15、*sin(0.9055*x-0.5845); I1=sum( y(1:499999)*(20)/(500000-1)I1 = 29.5222所以面積為29.5222 ()。 2）右矩形公式求積當(dāng)分為500個矩形時： 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個矩形時： x=linspace(0,20,5000); y= 2

16、.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:5000)*(20)/(5000-1)I1 = 29.5222當(dāng)分為50000個矩形時： 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所以面

17、積為29.5222 ()。（2）梯形公式求積： 1）直接去離散點進(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;trapz(x,y)ans = 29.4000 2）運用擬合曲線進(jìn)行求積：當(dāng)分為500個梯形時： 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個梯

18、形時： 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);I3=trapz(x,y)I3 = 29.5222當(dāng)分為50000個梯形時： 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 ()

19、。（3）用辛普森公式求積：I4=quad(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é)本文主要討論了對實際問題的數(shù)值積分計算方法，并通過MATLAB軟件編程實現(xiàn)，通過前面的研究我們知道求數(shù)值積分近似值的計算方法很多，有矩形求積公式、梯形求積公式、辛普森求積公式和相應(yīng)的復(fù)合公式等等。 牛頓-科特斯方法是一種利用插值多項式來構(gòu)造數(shù)值積分的常用方法，

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

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

22、 Applications，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 fittin

23、g% window, 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 =

24、 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 488);legh_ = ; legt_ = ; % handles and text

25、 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. xx = x(:);y = y(:);% This dataset does not appear on the plot% Add it to the plot by removing the if/end statements

26、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),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,Paren

27、t,ax_,Color,0.333333 0.666667 0,. LineStyle,none, LineWidth,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 1fo_ = fitoptions(method,NonlinearLeastSquares,Lower,-Inf 0 -Inf -Inf 0 -Inf -Inf 0 -Inf );ok_ = (isnan(x) | isnan(y);st_