复化梯形算法求解数值积分

1、复化梯形算法求解数值积分复化梯形算法求解数值积分摘 要求某函数的定积分时，在多数情况下，被积函数的原函数很难用初等函数表达出来，因此能够借助微积分学的牛顿-莱布尼兹公式计算定积分的机会是不多的。另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解。由于以上原因，数值积分的理论与方法一直是计算数学研究的基本课题。构造数值积分公式最通常的方法是用积分区间上的n 次插值多项式代替被积函数，由此导出的求积公式称为插值型求积公式。特别在节点分布等距的情形称为牛顿-柯茨公式，例如梯形公式与抛物线公式就是最基本的近似公式。但是它们的精度较差。而且高

2、阶Newton-Cotes求积公式是不稳定的。因此，通常不用高阶求积公式得到比较精确的积分值，而是将整个积分区间分段，在每一小段上用低阶求积公式。这种方法称为复化求积方法。本文从三个积分实例出发，主要讨论复化梯形公式以及精确程度分析。关键词：数值积分；复化求积公式；复化梯形算法；MATLAB复化梯形算法求解数值积分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复化梯形算法求解数值积分目 录1 问题的提出.12 问题的分析.23 问题假设.24 符号说明.35 模型的建立及求解.35.1 模型的准备工作.35.1.1 复化梯形数值积分基本原理. .35.2 模型的建立及求解.46 模型验证及结果分析.8参考文献.9附录.10复化梯形算法求解数值积分1 问题提出有很多实际问题常常需要计算积分才能求解。有些数值方法，如微分方程和积分方程的求解，也都和积分计算有联系。依据人们所熟知的微积分基本原理，对于积分I=的原函数F（x），便有下列Newton-Liebniz公式：=F(b)-F(a) ，只要找到被积函数f(x)但实际使用这种求积方法往往有困难，

9、因为大量的被积函数，诸如，等等，找不到初等函数表示的原函数；另外，当f(x)是由测量或数值计算给出的一张数据表时，Newton-Leibniz公式也不能直接运用。因此有必要研究积分的数值计算问题。下面研究：1）数值积分的编程实现2）编写复合梯形算法的程序分别以三个例子 1) 2) 3)进行计算和讨论，注意：不允许用MATLAB内部函数；用复化梯形公式时与MATLAB中内部的quad、quadl命令函数比较结果。第 1 页 共 13页复化梯形算法求解数值积分2 问题分析在科学实验和生产实践中，经常要求函数S=在区间a,b上的定积分： 的原函数在高等数学中，计算积分依靠微积分基本定理，只要找到被积

10、函数F(x)，则可用牛顿-莱布尼茨（Newton-Leibniz）公式：来求出定积分。但是，在有些情况下，应用牛顿-莱布尼茨公式往往有困难，例如，当被积函数的原函数无法用初等函数表示或被积函数为仅知离散点处函数值的离散函数时。类似于积分问题，在微分学中函数的导数是用极限来定义的，如果一个函数是以数值给出的离散形式，那么它的导数是用极限来定义的，如果一个函数是以数值给出的离散形式，那么它的导数就无法用极限运算方法来求得，当然也就无法用求导方法去计算函数在某点处的导数。基于以上原因，在许多实际问题中要采取数值方法来求函数的积分或微分，下面我们介绍积分与微分数值方法的一种基本思想-复化梯形算法，并讨

11、论其在MATLAB中的实现方法。3 问题假设（1）设被积函数都是二次可积的。（2）所有要积分的问题都是定积分。第 2 页 共 13页复化梯形算法求解数值积分（3）积分过程中计算机内寸足够且计算机运转正常。4 符号说明5 模型的建立以及求解5.1 模型的准备工作5.1.1 复化梯形数值积分基本原理求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)法、牛顿柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间a,b分成n个子区间xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题：S=而在每一个小的区间

12、上定积分的值可以近似求得。采用复化梯形求积算法时，将积分区间a,b划分为n等份，步长h=,分点为=a+kh，k=0,1,n。所谓复化求积法，就是先用低阶的Newton-Cotes公式求得每个子区间,上的积分值，然后再求和，用第 3 页 共 13页作为所求积分I的近似值。复化梯形算法求解数值积分复化梯形公式的形式是：=其积分余项：=5.2 建立模型对于以下三个例子，我们把整个积分区间a,b分成n个子区间xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题： S=复化梯形公式的形式是=其积分余项= = = =对于一般的求定积分，MATLAB基于变步长辛普生法

13、给出了quad函数和quadl函数来求定积分。函数的调用格式为：I,n=quad(fname,a,b,tol,trace)I,n=quadl(fname,a,b,tol,trace)其中，fname是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，默认时取tol=。trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，默认时取trace=0。返回参数I即定积分值，n为被积函数的调用次数。 （一）1）数值积分的编程实现直接利用quad函数，MATLAB函数如下：format long第 4 页 共 13页复化梯形算法求解数值积分fx=inline('x

14、.*exp(x)');I,n=quad(fx ,1,2,1e-10)运行结果：I = 7.389056098930670 n = 129直接利用quadl函数，MATLAB函数如下：format longfx=inline('x.*exp(x)');I,n=quadl(fx ,1,2,1e-10)运行结果：I = 7.389056098934960 n = 182）编写复合梯形算法的程序求出f（x）的n阶导数，MATLAB函数如下：syms xf=x*exp(x) %定义函数f（x）n=input('输入所求导数阶数:')f2=diff(f,x,n) %

15、求f(x)的n阶导数输入： n = 2求得： f2 = 2*exp(x)+x*exp(x)复化梯形算法计算数值积分：具体程序请见附录（程序1）用复化梯形算法计算的结果：积分结果 Tn= 7.389056127230221 等分数 n= 7019已知值与计算值的误差 R= 2.829957068684053e-008 （二）1）数值积分的编程实现第 5 页 共 13页复化梯形算法求解数值积分直接利用quad函数，MATLAB函数如下：format longfx=inline('x.2');I,n=quad(fx ,0,1,1e-10)运行结果：I = 0.333333333333

16、333 n = 13直接利用quadl函数，MATLAB函数如下：format longfx=inline('x.2');I,n=quadl(fx ,0,1,1e-10)运行结果：I = 0.333333333333333 n = 182）编写复合梯形算法的程序求出f（x）的n阶导数，MATLAB函数如下：syms xf=x2 %定义函数f（x）n=input('输入所求导数阶数:')f2=diff(f,x,n) %求f(x)的n阶导数输入： n = 2求得： f2 = 2复化梯形算法计算数值积分：具体程序请见附录（程序2）用复化梯形算法计算的结果：积分结果 T

17、n= 0.333333383319197 等分数 n= 1826已知值与计算值的误差 R= 4.998586394799887e-008 （三）第 6 页 共 13页复化梯形算法求解数值积分1）数值积分的编程实现直接利用quad函数，MATLAB函数如下：format longfx=inline('sin(x)');I,n=quad(fx,0,1,1e-10)运行结果：I = 0.459697694131815 n = 61直接利用quadl函数，MATLAB函数如下：format longfx=inline('sin(x)');I,n=quadl(fx,0,1

18、,1e-10)运行结果：I = 0.459697694131821 n = 182）编写复合梯形算法的程序求出f（x）的n阶导数，MATLAB函数如下：syms xf=sin(x) %定义函数f（x）n=input('输入所求导数阶数:')f2=diff(f,x,n) %求f(x)的n阶导数输入： n = 2求得： f2 = -sin(x)复化梯形算法计算数值积分：具体程序请见附录（程序3）用复化梯形算法计算的结果：积分结果 Tn= 0.459697666851231 等分数 n= 1185已知值与计算值的误差 R= -2.728059006473771e-008第 7 页 共

19、 13页复化梯形算法求解数值积分6 模型验证及结果分析经验证运用复化梯形公式求解数值分析积分结果对比表如下：发现运用复化梯形公式求解积分效果较好，绝对误差较小。但是从理论上，我们分析，复合梯形公式的误差限的系数为-（b-a）/12，而且后面是h平方级还有f的二次导数，而辛普森的系数是-（b-a）/（180*16），后面h是四次方级的，f的导数为四次导数，显然辛普森的误差限更加小而高斯求积公式是对代数精度的方面有着更加好的结果。第 8 页 共 13页复化梯形算法求解数值积分参 考 文 献1 姜启源，谢金星，叶俊编数学模型（第三版）M北京：高等教育出版社，2005：1-2022 李庆阳，王能超，易

20、大义编数值分析（第四版）M武汉：华中科技大学出版社，20063 清华大学、北京大学计算方法编写组计算方法M北京：科学出版社，19804 冯康，等编数值计算方法M北京：国防工业出版社，19785 苏今明，阮沈勇编MATLAB实用教程（第二版）M北京：电子工业出版社，20086 刘卫国主编MATLAB程序设计教程（第二版）M . 北京：中国水利水电出版社，2010第 9 页 共 13页复化梯形算法求解数值积分附 录程序1：clcclearsyms x %定义自变量xf=inline('x*exp(x)','x') %定义函数f(x)=x*exp(x)，换函数时只需换

21、该函数表达式即可 f2=inline('(2*exp(x) + x*exp(x)','x') %定义f(x)的二阶导数,输入程序1里求出的f2即可。f3='-(2*exp(x) + x*exp(x)' %因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，以便求最大值e=5*10(-8) %精度要求值a=1 %积分下限b=2 %积分上限x1=fminbnd(f3,1,2) %求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值for n=2:1000000 %求等分数nRn=-(b-a)/12*(b-a)/n)2

22、*f2(x1) %计算余项if abs(Rn)<e %用余项进行判断break % 符合要求时结束endendh=(b-a)/n %求hTn1=0for k=1:n-1 %求连加和xk=a+k*hTn1=Tn1+f(xk)共 13页第 10 页复化梯形算法求解数值积分endTn=h/2*(f(a)+2*Tn1+f(b)z=exp(2)R=Tn-z %求已知值与计算值的差fprintf('用复化梯形算法计算的结果 Tn=')disp(Tn)fprintf('等分数 n=')disp(n) %输出等分数fprintf('已知值与计算值的误差 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') %定义函数f(x)=x2，换函数时只需换该函数表达式即可f2=inline('(2)','x') %定义f(x)的二阶导数,输入程序1里求出的f2即可。f3='-(2)' %因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，以便求最大值e=5*10(-8) %精度要

24、求值a=0 %积分下限b=1 %积分上限x1=fminbnd(f3,1,2) %求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值for n=2:1000000 %求等分数nRn=-(b-a)/12*(b-a)/n)2*f2(x1) %计算余项if abs(Rn)<e %用余项进行判断共 13页第 11 页复化梯形算法求解数值积分break % 符合要求时结束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 %求已知值与计算值的差fprintf('用复化梯形算法计算的结果 Tn=')disp(Tn)fprintf('等分