復(fù)化梯形算法求解數(shù)值積分

1、復(fù)化梯形算法求解數(shù)值積分復(fù)化梯形算法求解數(shù)值積分摘 要求某函數(shù)的定積分時(shí)，在多數(shù)情況下，被積函數(shù)的原函數(shù)很難用初等函數(shù)表達(dá)出來(lái)，因此能夠借助微積分學(xué)的牛頓-萊布尼茲公式計(jì)算定積分的機(jī)會(huì)是不多的。另外，許多實(shí)際問(wèn)題中的被積函數(shù)往往是列表函數(shù)或其他形式的非連續(xù)函數(shù)，對(duì)這類函數(shù)的定積分，也不能用不定積分方法求解。由于以上原因，數(shù)值積分的理論與方法一直是計(jì)算數(shù)學(xué)研究的基本課題。構(gòu)造數(shù)值積分公式最通常的方法是用積分區(qū)間上的n 次插值多項(xiàng)式代替被積函數(shù)，由此導(dǎo)出的求積公式稱為插值型求積公式。特別在節(jié)點(diǎn)分布等距的情形稱為牛頓-柯茨公式，例如梯形公式與拋物線公式就是最基本的近似公式。但是它們的精度較差。而且高

2、階Newton-Cotes求積公式是不穩(wěn)定的。因此，通常不用高階求積公式得到比較精確的積分值，而是將整個(gè)積分區(qū)間分段，在每一小段上用低階求積公式。這種方法稱為復(fù)化求積方法。本文從三個(gè)積分實(shí)例出發(fā)，主要討論復(fù)化梯形公式以及精確程度分析。關(guān)鍵詞：數(shù)值積分；復(fù)化求積公式；復(fù)化梯形算法；MATLAB復(fù)化梯形算法求解數(shù)值積分THE REHABILITATION OF TRAPEZOID FORMULA TOSOLVE THE NUMERICAL INTEGRATIONABSTRACTFind the definite integral of a function, in most cases, the

3、original integrand function is difficult toexpress the elementary functions, it can use calculus of Newton - Leibniz formula to calculate thedefinite integral of the few opportunities . In addition, many practical problems in the integrand is often a list of functions or other forms of non-continuou

4、s function, the definite integral of suchfunctions, indefinite integral method can not solve. For these reasons, the numerical integration oftheory and method has been the subject of calculation of the basic mathematical research.Structural formula for numerical integration method is used most often

5、 on the n-th integration interval polynomial interpolation instead of the integrand, thus derived is called interpolation-typequadrature formula quadrature formula. Especially in the case of equidistant distribution of nodesis called Newton - Keci formula, such as trapezoidal formula and the formula

6、 is the most basicparabolic approximation formula. But their accuracy is poor. And high-level Newton-Cotesquadrature formula is unstable. So it is usually not higher-order quadrature formula to be moreprecise integral values, but the whole range of sub-points, with each short on low-level quadrature

7、 formula. This method is called complex method of quadrature.This example from three points of departure, the main complex of the trapezoid formula anddiscuss the accuracy of the analysis.Key words: Numerical integration；Rehabilitation of numerical integration；Rehabilitation of trapezoid formula；MAT

8、LAB復(fù)化梯形算法求解數(shù)值積分目 錄1 問(wèn)題的提出.12 問(wèn)題的分析.23 問(wèn)題假設(shè).24 符號(hào)說(shuō)明.35 模型的建立及求解.35.1 模型的準(zhǔn)備工作.35.1.1 復(fù)化梯形數(shù)值積分基本原理. .35.2 模型的建立及求解.46 模型驗(yàn)證及結(jié)果分析.8參考文獻(xiàn).9附錄.10復(fù)化梯形算法求解數(shù)值積分1 問(wèn)題提出有很多實(shí)際問(wèn)題常常需要計(jì)算積分才能求解。有些數(shù)值方法，如微分方程和積分方程的求解，也都和積分計(jì)算有聯(lián)系。依據(jù)人們所熟知的微積分基本原理，對(duì)于積分I=的原函數(shù)F（x），便有下列Newton-Liebniz公式：=F(b)-F(a) ，只要找到被積函數(shù)f(x)但實(shí)際使用這種求積方法往往有困難，

9、因?yàn)榇罅康谋环e函數(shù)，諸如，等等，找不到初等函數(shù)表示的原函數(shù)；另外，當(dāng)f(x)是由測(cè)量或數(shù)值計(jì)算給出的一張數(shù)據(jù)表時(shí)，Newton-Leibniz公式也不能直接運(yùn)用。因此有必要研究積分的數(shù)值計(jì)算問(wèn)題。下面研究：1）數(shù)值積分的編程實(shí)現(xiàn)2）編寫(xiě)復(fù)合梯形算法的程序分別以三個(gè)例子 1) 2) 3)進(jìn)行計(jì)算和討論，注意：不允許用MATLAB內(nèi)部函數(shù)；用復(fù)化梯形公式時(shí)與MATLAB中內(nèi)部的quad、quadl命令函數(shù)比較結(jié)果。第 1 頁(yè) 共 13頁(yè)復(fù)化梯形算法求解數(shù)值積分2 問(wèn)題分析在科學(xué)實(shí)驗(yàn)和生產(chǎn)實(shí)踐中，經(jīng)常要求函數(shù)S=在區(qū)間a,b上的定積分： 的原函數(shù)在高等數(shù)學(xué)中，計(jì)算積分依靠微積分基本定理，只要找到被積

10、函數(shù)F(x)，則可用牛頓-萊布尼茨（Newton-Leibniz）公式：來(lái)求出定積分。但是，在有些情況下，應(yīng)用牛頓-萊布尼茨公式往往有困難，例如，當(dāng)被積函數(shù)的原函數(shù)無(wú)法用初等函數(shù)表示或被積函數(shù)為僅知離散點(diǎn)處函數(shù)值的離散函數(shù)時(shí)。類似于積分問(wèn)題，在微分學(xué)中函數(shù)的導(dǎo)數(shù)是用極限來(lái)定義的，如果一個(gè)函數(shù)是以數(shù)值給出的離散形式，那么它的導(dǎo)數(shù)是用極限來(lái)定義的，如果一個(gè)函數(shù)是以數(shù)值給出的離散形式，那么它的導(dǎo)數(shù)就無(wú)法用極限運(yùn)算方法來(lái)求得，當(dāng)然也就無(wú)法用求導(dǎo)方法去計(jì)算函數(shù)在某點(diǎn)處的導(dǎo)數(shù)。基于以上原因，在許多實(shí)際問(wèn)題中要采取數(shù)值方法來(lái)求函數(shù)的積分或微分，下面我們介紹積分與微分?jǐn)?shù)值方法的一種基本思想-復(fù)化梯形算法，并討

11、論其在MATLAB中的實(shí)現(xiàn)方法。3 問(wèn)題假設(shè)（1）設(shè)被積函數(shù)都是二次可積的。（2）所有要積分的問(wèn)題都是定積分。第 2 頁(yè) 共 13頁(yè)復(fù)化梯形算法求解數(shù)值積分（3）積分過(guò)程中計(jì)算機(jī)內(nèi)寸足夠且計(jì)算機(jī)運(yùn)轉(zhuǎn)正常。4 符號(hào)說(shuō)明5 模型的建立以及求解5.1 模型的準(zhǔn)備工作5.1.1 復(fù)化梯形數(shù)值積分基本原理求解定積分的數(shù)值方法多種多樣，如簡(jiǎn)單的梯形法、辛普生(Simpson)法、牛頓柯特斯(Newton-Cotes)法等都是經(jīng)常采用的方法。它們的基本思想都是將整個(gè)積分區(qū)間a,b分成n個(gè)子區(qū)間xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。這樣求定積分問(wèn)題就分解為求和問(wèn)題：S=而在每一個(gè)小的區(qū)間

12、上定積分的值可以近似求得。采用復(fù)化梯形求積算法時(shí)，將積分區(qū)間a,b劃分為n等份，步長(zhǎng)h=,分點(diǎn)為=a+kh，k=0,1,n。所謂復(fù)化求積法，就是先用低階的Newton-Cotes公式求得每個(gè)子區(qū)間,上的積分值，然后再求和，用第 3 頁(yè) 共 13頁(yè)作為所求積分I的近似值。復(fù)化梯形算法求解數(shù)值積分復(fù)化梯形公式的形式是：=其積分余項(xiàng)：=5.2 建立模型對(duì)于以下三個(gè)例子，我們把整個(gè)積分區(qū)間a,b分成n個(gè)子區(qū)間xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。這樣求定積分問(wèn)題就分解為求和問(wèn)題： S=復(fù)化梯形公式的形式是=其積分余項(xiàng)= = = =對(duì)于一般的求定積分，MATLAB基于變步長(zhǎng)辛普生法

13、給出了quad函數(shù)和quadl函數(shù)來(lái)求定積分。函數(shù)的調(diào)用格式為：I,n=quad(fname,a,b,tol,trace)I,n=quadl(fname,a,b,tol,trace)其中，fname是被積函數(shù)名。a和b分別是定積分的下限和上限。tol用來(lái)控制積分精度，默認(rèn)時(shí)取tol=。trace控制是否展現(xiàn)積分過(guò)程，若取非0則展現(xiàn)積分過(guò)程，取0則不展現(xiàn)，默認(rèn)時(shí)取trace=0。返回參數(shù)I即定積分值，n為被積函數(shù)的調(diào)用次數(shù)。 （一）1）數(shù)值積分的編程實(shí)現(xiàn)直接利用quad函數(shù)，MATLAB函數(shù)如下：format long第 4 頁(yè) 共 13頁(yè)復(fù)化梯形算法求解數(shù)值積分fx=inline('x

14、.*exp(x)');I,n=quad(fx ,1,2,1e-10)運(yùn)行結(jié)果：I = 7.389056098930670 n = 129直接利用quadl函數(shù)，MATLAB函數(shù)如下：format longfx=inline('x.*exp(x)');I,n=quadl(fx ,1,2,1e-10)運(yùn)行結(jié)果：I = 7.389056098934960 n = 182）編寫(xiě)復(fù)合梯形算法的程序求出f（x）的n階導(dǎo)數(shù)，MATLAB函數(shù)如下：syms xf=x*exp(x) %定義函數(shù)f（x）n=input('輸入所求導(dǎo)數(shù)階數(shù):')f2=diff(f,x,n) %

15、求f(x)的n階導(dǎo)數(shù)輸入： n = 2求得： f2 = 2*exp(x)+x*exp(x)復(fù)化梯形算法計(jì)算數(shù)值積分：具體程序請(qǐng)見(jiàn)附錄（程序1）用復(fù)化梯形算法計(jì)算的結(jié)果：積分結(jié)果 Tn= 7.389056127230221 等分?jǐn)?shù) n= 7019已知值與計(jì)算值的誤差 R= 2.829957068684053e-008 （二）1）數(shù)值積分的編程實(shí)現(xiàn)第 5 頁(yè) 共 13頁(yè)復(fù)化梯形算法求解數(shù)值積分直接利用quad函數(shù)，MATLAB函數(shù)如下：format longfx=inline('x.2');I,n=quad(fx ,0,1,1e-10)運(yùn)行結(jié)果：I = 0.333333333333

16、333 n = 13直接利用quadl函數(shù)，MATLAB函數(shù)如下：format longfx=inline('x.2');I,n=quadl(fx ,0,1,1e-10)運(yùn)行結(jié)果：I = 0.333333333333333 n = 182）編寫(xiě)復(fù)合梯形算法的程序求出f（x）的n階導(dǎo)數(shù)，MATLAB函數(shù)如下：syms xf=x2 %定義函數(shù)f（x）n=input('輸入所求導(dǎo)數(shù)階數(shù):')f2=diff(f,x,n) %求f(x)的n階導(dǎo)數(shù)輸入： n = 2求得： f2 = 2復(fù)化梯形算法計(jì)算數(shù)值積分：具體程序請(qǐng)見(jiàn)附錄（程序2）用復(fù)化梯形算法計(jì)算的結(jié)果：積分結(jié)果 T

17、n= 0.333333383319197 等分?jǐn)?shù) n= 1826已知值與計(jì)算值的誤差 R= 4.998586394799887e-008 （三）第 6 頁(yè) 共 13頁(yè)復(fù)化梯形算法求解數(shù)值積分1）數(shù)值積分的編程實(shí)現(xiàn)直接利用quad函數(shù)，MATLAB函數(shù)如下：format longfx=inline('sin(x)');I,n=quad(fx,0,1,1e-10)運(yùn)行結(jié)果：I = 0.459697694131815 n = 61直接利用quadl函數(shù)，MATLAB函數(shù)如下：format longfx=inline('sin(x)');I,n=quadl(fx,0,1

18、,1e-10)運(yùn)行結(jié)果：I = 0.459697694131821 n = 182）編寫(xiě)復(fù)合梯形算法的程序求出f（x）的n階導(dǎo)數(shù)，MATLAB函數(shù)如下：syms xf=sin(x) %定義函數(shù)f（x）n=input('輸入所求導(dǎo)數(shù)階數(shù):')f2=diff(f,x,n) %求f(x)的n階導(dǎo)數(shù)輸入： n = 2求得： f2 = -sin(x)復(fù)化梯形算法計(jì)算數(shù)值積分：具體程序請(qǐng)見(jiàn)附錄（程序3）用復(fù)化梯形算法計(jì)算的結(jié)果：積分結(jié)果 Tn= 0.459697666851231 等分?jǐn)?shù) n= 1185已知值與計(jì)算值的誤差 R= -2.728059006473771e-008第 7 頁(yè) 共

19、 13頁(yè)復(fù)化梯形算法求解數(shù)值積分6 模型驗(yàn)證及結(jié)果分析經(jīng)驗(yàn)證運(yùn)用復(fù)化梯形公式求解數(shù)值分析積分結(jié)果對(duì)比表如下：發(fā)現(xiàn)運(yùn)用復(fù)化梯形公式求解積分效果較好，絕對(duì)誤差較小。但是從理論上，我們分析，復(fù)合梯形公式的誤差限的系數(shù)為-（b-a）/12，而且后面是h平方級(jí)還有f的二次導(dǎo)數(shù)，而辛普森的系數(shù)是-（b-a）/（180*16），后面h是四次方級(jí)的，f的導(dǎo)數(shù)為四次導(dǎo)數(shù)，顯然辛普森的誤差限更加小而高斯求積公式是對(duì)代數(shù)精度的方面有著更加好的結(jié)果。第 8 頁(yè) 共 13頁(yè)復(fù)化梯形算法求解數(shù)值積分參 考 文 獻(xiàn)1 姜啟源，謝金星，葉俊編數(shù)學(xué)模型（第三版）M北京：高等教育出版社，2005：1-2022 李慶陽(yáng)，王能超，易

20、大義編數(shù)值分析（第四版）M武漢：華中科技大學(xué)出版社，20063 清華大學(xué)、北京大學(xué)計(jì)算方法編寫(xiě)組計(jì)算方法M北京：科學(xué)出版社，19804 馮康，等編數(shù)值計(jì)算方法M北京：國(guó)防工業(yè)出版社，19785 蘇今明，阮沈勇編MATLAB實(shí)用教程（第二版）M北京：電子工業(yè)出版社，20086 劉衛(wèi)國(guó)主編MATLAB程序設(shè)計(jì)教程（第二版）M . 北京：中國(guó)水利水電出版社，2010第 9 頁(yè) 共 13頁(yè)復(fù)化梯形算法求解數(shù)值積分附 錄程序1：clcclearsyms x %定義自變量xf=inline('x*exp(x)','x') %定義函數(shù)f(x)=x*exp(x)，換函數(shù)時(shí)只需換

21、該函數(shù)表達(dá)式即可 f2=inline('(2*exp(x) + x*exp(x)','x') %定義f(x)的二階導(dǎo)數(shù),輸入程序1里求出的f2即可。f3='-(2*exp(x) + x*exp(x)' %因fminbnd（）函數(shù)求的是表達(dá)式的最小值，且要求表達(dá)式帶引號(hào)，故取負(fù)號(hào)，以便求最大值e=5*10(-8) %精度要求值a=1 %積分下限b=2 %積分上限x1=fminbnd(f3,1,2) %求負(fù)的二階導(dǎo)數(shù)的最小值點(diǎn)，也就是求二階導(dǎo)數(shù)的最大值點(diǎn)對(duì)應(yīng)的x值for n=2:1000000 %求等分?jǐn)?shù)nRn=-(b-a)/12*(b-a)/n)2

22、*f2(x1) %計(jì)算余項(xiàng)if abs(Rn)<e %用余項(xiàng)進(jìn)行判斷break % 符合要求時(shí)結(jié)束endendh=(b-a)/n %求hTn1=0for k=1:n-1 %求連加和xk=a+k*hTn1=Tn1+f(xk)共 13頁(yè)第 10 頁(yè)復(fù)化梯形算法求解數(shù)值積分endTn=h/2*(f(a)+2*Tn1+f(b)z=exp(2)R=Tn-z %求已知值與計(jì)算值的差fprintf('用復(fù)化梯形算法計(jì)算的結(jié)果 Tn=')disp(Tn)fprintf('等分?jǐn)?shù) n=')disp(n) %輸出等分?jǐn)?shù)fprintf('已知值與計(jì)算值的誤差 R=

23、9;)disp(R)程序2：clcclearformat longfx=inline('x.2');I,n=quadl(fx,0,1,1e-10)syms x %定義自變量xf=inline('x2','x') %定義函數(shù)f(x)=x2，換函數(shù)時(shí)只需換該函數(shù)表達(dá)式即可f2=inline('(2)','x') %定義f(x)的二階導(dǎo)數(shù),輸入程序1里求出的f2即可。f3='-(2)' %因fminbnd（）函數(shù)求的是表達(dá)式的最小值，且要求表達(dá)式帶引號(hào)，故取負(fù)號(hào)，以便求最大值e=5*10(-8) %精度要

24、求值a=0 %積分下限b=1 %積分上限x1=fminbnd(f3,1,2) %求負(fù)的二階導(dǎo)數(shù)的最小值點(diǎn)，也就是求二階導(dǎo)數(shù)的最大值點(diǎn)對(duì)應(yīng)的x值for n=2:1000000 %求等分?jǐn)?shù)nRn=-(b-a)/12*(b-a)/n)2*f2(x1) %計(jì)算余項(xiàng)if abs(Rn)<e %用余項(xiàng)進(jìn)行判斷共 13頁(yè)第 11 頁(yè)復(fù)化梯形算法求解數(shù)值積分break % 符合要求時(shí)結(jié)束endendh=(b-a)/n %求hTn1=0for k=1:n-1 %求連加和xk=a+k*hTn1=Tn1+f(xk)endTn=h/2*(f(a)+2*Tn1+f(b)z=IR=Tn-z %求已知值與計(jì)算值的差fprintf('用復(fù)化梯形算法計(jì)算的結(jié)果 Tn=')disp(Tn)fprintf('等分