数值积分及实现

演示文稿数值积分及实现当前1页，总共42页。（优选）数值积分及实现当前2页，总共42页。数值积分和数值微分1引言

我们知道,若函数f(x)在区间[a,b]上连续且其原函数为F(x),则可用Newton-Leibnitz公式求得定积分求定积分的值,Newton-Leibnitz公式无论在理论上还是在解决实际问题上都起了很大作用，但它并不能完全解决定积分的计算问题，因为积分学涉及的实际问题极为广泛，而且极其复杂，在实际计算中经常遇到以下三种情况：当前3页，总共42页。

(1)被积函数f(x)并不一定能够找到用初等函数的有限形式表示的原函数F(x)，例如：

Newton-Leibnitz公式就无能为力了(2)还有被积函数f(x)的原函数能用初等函数表示，但表达式太复杂，例如函数并不复杂，但积分后其表达式却很复杂，积分后其原函数F(x)为：当前4页，总共42页。(3)被积函数f(x)没有具体的解析表达式,其函数关系由表格或图形表示。对于这些情况,要计算积分的准确值都是十分困难的。由此可见,通过原函数来计算积分有它的局限性,因而研究一种新的积分方法来解决Newton-Leibniz公式所不能或很难解决的积分问题,这时需要用数值解法来建立积分的近似计算方法。将积分区间细分,在每一个小区间内用简单函数代替复杂函数进行积分，这就是数值积分的思想，用代数插值多项式去代替被积函数发f(x)进行积分是本章讨论数值积分的主要内容。

当前5页，总共42页。

建立数值积分公式的途径比较多,其中最常用的有两种：（1）由积分中值定理可知，对于连续函数f(x)，在积分区间[a,b]内存在一点ξ，使得即所求的曲边梯形的面积恰好等于底为(b-a)，高为的矩形面积。但是点ξ的具体位置一般是未知的,因而的值也是未知的,称为f(x)在区间[a,b]上的平均高度。那么只要对平均高度提供一种算法，相应地就获得一种数值求积方法当前6页，总共42页。三个求积分公式①梯形公式y=f(x)yxaby=f(x)abyx(a+b)/2②中矩形公式按照这种思想，可构造出一些求积分值的近似公式。例如分别取和则分别得到中矩形公式和梯形公式。y=f(x)abab当前7页，总共42页。y=f(x)yab③Simpson公式(a+b)/2f()的近似值而获得的一种数值积分方法。中矩形公式把[a,b]

的中点处函数值作为平均高度f()的近似值而获得的一种数值积分方法。

ab(a+b)/2

在这三个公式中,梯形公式把f(a),f(b)的加权平均值作为平均高度当前8页，总共42页。Simpson公式是以函数f(x)在a,b,(a+b)/2这三点的函数值f(a),f(b),的加权平均值似值而获得的一种数值积分方法。作为平均高度f()的近（2）先用某个简单函数近似逼近f(x),用代替原被积函数f(x)，即

以此构造数值算法。从数值计算的角度考虑,函数应对f(x)有充分的逼近程度,并且容易计算其积分。由于多项式能很好地逼近连续函数,且又容易计算积分,因此将选取为插值多项式,这样f(x)的积分就可以用其插值多项式的积分来近似代替当前9页，总共42页。2.2插值求积公式

设已知f(x)在节点有函数值,作n次拉格朗日插值多项式式中这里多项式P(x)易于求积,所以可取作为的近似值，即当前10页，总共42页。其中

称为求积系数。给出如下定义。定义1求积公式其系数时，则称求积公式为插值求积公式。(4)当前11页，总共42页。设插值求积公式的余项为,由插值余项定理得

其中

当f(x)是次数不高于n的多项式时，有

=0,求积公式(4)能成为准确的等式。由于闭区间[a,b]上的连续函数可用多项式逼近，所以一个求积公式能对多大次数的多项式f(x)成为准确等式，是衡量该公式的精确程度的重要指标，为此给出以下定义。当前12页，总共42页。定义2（代数精度）

设求积公式（4）对于一切次数小于等于m的多项式(是准确的，而对于次数为m+1的多项式是不准确的，则称该求积公式具有m次代数精度（简称代数精度）或)当前13页，总共42页。定理1n+1个节点的求积公式为插值型求积公式的充要条件是公式至少具有n次代数精度。

当前14页，总共42页。例1设积分区间[a,b]为[0,2]，取时时,

分别用梯形和辛卜生公式

计算其积分结果并与准确值进行比较解:梯形公式和辛卜生的计算结果与准确值比较如下表所示当前15页，总共42页。

f(x)1xx2x3x4ex

准确值222.6746.406.389

梯形公式计算值2248168.389

辛卜生公式计算值222.6746.676.421

从表中可以看出,当f(x)是时,辛卜生公式比梯形公式更精确

一般说来，代数精度越高，求积公式越精确。梯形公式和中矩形公式具有1次代数精度，辛卜生公式有3次代数精度。下面以梯形公式为例进行验证当前16页，总共42页。取f(x)=1时，

两端相等取f(x)=x时,取f(x)=x2时,两端不相等所以梯形公式只有1次代数精度。两端相等当前17页，总共42页。构造插值求积公式有如下特点：复杂函数f(x)的积分转化为计算多项式的积分求积系数Ak只与积分区间及节点xk有关，而与被积函数f(x)无关，可以不管f(x)如何，预先算出Ak的值

n+1个节点的插值求积公式至少具有n次代数精度求积系数之和可用此检验计算求积系数的正确性

当前18页，总共42页。3牛顿—柯特斯(Newton-Cotes)求积公式

在插值求积公式中,当所取节点是等距时称为牛顿-柯特斯公式其中插值多项式求积系数这里是插值基函数。即有当前19页，总共42页。将积分区间[a,b]

划分为n等分,步长求积节点为为了计算系数Ak,由于,所以作变量代换当时,有

,于是可得

当前20页，总共42页。

(k=0,1…,n)

代入插值求积公式(4)有

称为牛顿-柯特斯求积公式,Ck称为柯特斯系数引进记号

(k=0,1…,n)

则当前21页，总共42页。容易验证

∵

∴

显然,Ck是不依赖于积分区间[a,b]以及被积函数f(x)的常数,只要给出n,就可以算出柯特斯系数,譬如当n=1时当前22页，总共42页。当n=2时

当前23页，总共42页。4几个低阶求积公式

在牛顿-柯特斯求积公式中n=1,2,4时，就分别得到下面的梯形公式、辛卜生公式和柯特斯公式。(1)

梯形公式当n=1时，牛顿-柯特斯公式就是梯形公式

定理2（梯形公式的误差）设f(x)在[a,b]上具有连续的二阶导数，则梯形公式的误差（余项）为当前24页，总共42页。（2）辛卜生公式当n=2时，牛顿-柯特斯公式就是辛卜生公式（或称抛物线公式）

定理3（辛卜生公式的误差）设在[a,b]上具有连续的四阶导数，则辛卜生求积公式的误差为定理证明从略。

当前25页，总共42页。（3）柯特斯公式。当n=4时，牛顿-柯特斯公式为

定理4（柯特斯公式的误差）设在[a,b]上具有连续的6阶导数，则柯特斯求积公式的误差为

定理的证明从略。

当前26页，总共42页。例11分别用梯形公式、辛卜生公式和柯特斯公式计算定积分

的近似值(计算结果取5位有效数字)(1)用梯形公式计算(2)用辛卜生公式当前27页，总共42页。(3)用柯特斯公式计算，系数为

积分的准确值为

可见，三个求积公式的精度逐渐提高。当前28页，总共42页。例12用辛卜生公式和柯特斯公式计算定积分的近似值,并估计其误差(计算结果取5位小数)解:辛卜生公式由于由辛卜生公式余项知其误差为

当前29页，总共42页。例12用辛卜生公式和柯特斯公式计算定积分的近似值,并估计其误差(计算结果取5位小数)解:柯特斯公式

知其误差为当前30页，总共42页。例12用辛卜生公式和柯特斯公式计算定积分的近似值,并估计其误差(计算结果取5位小数)

该定积分的准确值,这个例子告诉我们，对于同一个积分，当n≥2时，公式却是精确的，这是由于辛卜生公式具有三次代数精度，柯特斯公式具有五次代数精度，它们对被积函数为三次多项式当然是精确成立的。当前31页，总共42页。数值积分基本原理求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)法、牛顿－柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间[a,b]分成n个子区间[xi,xi+1]，i=1,2,…,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。当前32页，总共42页。数值积分的实现方法1．变步长辛普生法基于变步长辛普生法，MATLAB给出了quad函数来求定积分。该函数的调用格式为：

[I,n]=quad('fname',a,b,tol,trace)其中fname是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0。返回参数I即定积分值，n为被积函数的调用次数。当前33页，总共42页。

例1求定积分。

(1)建立被积函数文件fesin.m。functionf=fesin(x)f=exp(-0.5*x).*sin(x+pi/6);(2)调用数值积分函数quad求定积分。[S,n]=quad('fesin',0,3*pi)S=0.9008n=77当前34页，总共42页。2．牛顿－柯特斯法基于牛顿－柯特斯法，MATLAB给出了quad8函数来求定积分。该函数的调用格式为：[I,n]=quad8('fname',a,b,tol,trace)其中参数的含义和quad函数相似，只是tol的缺省值取10-6。该函数可以更精确地求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。当前35页，总共42页。例2求定积分。(1)被积函数文件fx.m。functionf=fx(x)f=x.*sin(x)./(1+cos(x).*cos(x));(2)调用函数quad8求定积分。I=quad8('fx',0,pi)I=2.4674当前36页，总共42页。例3分别用quad函数和quad8函数求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。调用函数quad求定积分：formatlong;fx=inline('exp(-x)');[I,n]=quad(fx,1,2.5,1e-10)I=0.28579444254766n=65当前37页，总共42页。

调用函数quad8求定积分：formatlong;fx=inline('exp(-x)');[I,n]=quad8(fx,1,2.5,1e-10)I=0.28579444254754n=33当前38页，总共42页。3．被积函数由一个表格定义在MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y)函数。其中向量X,Y定义函数关系Y=f(X)。例4用trapz函数计算定积分。命令如下：X=1:0.01:2.5;Y=exp(-X);%生成函数关系数据向量trapz(X,Y)ans=0.28579682416393当前39页，总共42页。1.3二重定积分的数值