数值积分及matlab实现.ppt

《数值积分及matlab实现.ppt》由会员分享，可在线阅读，更多相关《数值积分及matlab实现.ppt（43页珍藏版）》请在三一办公上搜索。

1、数值积分与微分,数值积分和数值微分,1 引言 我们知道,若函数f(x)在区间a,b上连续且其原函数为F(x),则可用Newton-Leibnitz公式,求得定积分,求定积分的值,Newton-Leibnitz公式 无论在理论上还是在解决实际问题上都起了很大作用，但它并不能完全解决定积分的计算问题，因为积分学涉及的实际问题极为广泛，而且极其复杂，在实际计算中经常遇到以下三种情况：,(1)被积函数f(x)并不一定能够找到用初等函数的 有限形式表示的原函数F(x)，例如：Newton-Leibnitz公式就无能为力了,(2)还有被积函数f(x)的原函数能用初等函数表示，但表达式太复杂，例如函数,并不

2、复杂，但积分后其表达式却很复杂，积分后其原函数F(x)为：,(3)被积函数f(x)没有具体的解析表达式,其函数 关系由表格或图形表示。对于这些情况,要计算积分的准确值都是十分困难的。由此可见,通过原函数来计算积分有它的局限性,因而研究一种新的积分方法来解决Newton-Leibniz公式所不能或很难解决的积分问题,这时需要用数值解法来建立积分的近似计算方法。将积分区间细分,在每一个小区间内用简单函数代替复杂函数进行积分，这就是数值积分的思想，用代数插值多项式去代替被积函数发f(x)进行积分是本章讨论数值积分的主要内容。,建立数值积分公式的途径比较多,其中最常用的有两种：（1）由积分中值定理可知

3、，对于连续函数f(x)，在积分区间a,b内存在一点，使得即所求的曲边梯形的面积恰好等于底为(b-a)，高为 的矩形面积。但是点的具体位置一般是未知的,因而 的值也是未知的,称 为f(x)在区间a,b上的平均高度。那么只要对平均高度 提供一种算法，相应地就获得一种数值求积方法,三个求积分公式,梯形公式,y=f(x),y,x,a,b,y=f(x),a,b,y,x,(a+b)/2,中矩形公式,按照这种思想，可构造出一些求积分值的近似公式。例如 分别取 和,则分别得到中矩形公式和梯形公式。,y=f(x),a,b,a,b,y=f(x),y,a,b,Simpson公式,(a+b)/2,f()的近似值而获得

4、的一种数值积分方法。中矩形公式把a,b 的中点处函数值 作为平均高度f()的近似值而获得的一种数值积分方法。,a,b,(a+b)/2,在这三个公式中,梯形公式把f(a),f(b)的加权平均值,作为平均高度,Simpson公式是以函数f(x)在a,b,(a+b)/2这三点的函数值f(a),f(b),的加权平均值 似值而获得的一种数值积分方法。,作为平均高度f()的近,（2）先用某个简单函数 近似逼近f(x),用 代替原被积函数f(x)，即,以此构造数值算法。从数值计算的角度考虑,函数 应对f(x)有充分的逼近程度,并且容易计算其积分。由于多项式能很好地逼近连续函数,且又容易计算积分,因此将 选取

5、为插值多项式,这样f(x)的积分就可以用其插值多项式的积分来近似代替,2.2 插值求积公式设已知f(x)在节点 有函数值,作n次拉格朗日插值多项式,式中,这里,多项式P(x)易于求积,所以可取 作为 的近似值，即,其中,称为求积系数。给出如下定义。,定义1 求积公式,其系数 时，则称求积公式为插值求积公式。,(4),设插值求积公式的余项为,由插值余项定理得,其中,当f(x)是次数不高于n的多项式时，有=0,求积公式(4)能成为准确的等式。由于闭区间a,b上的连续函数可用多项式逼近，所以一个求积公式能对多大次数的多项式f(x)成为准确等式，是衡量该公式的精确程度的重要指标，为此给出以下定义。,定

6、义2（代数精度）设求积公式（4）对于一 切次数小于等于m的多项式(,是准确的，而对于次数为m+1的多项式是不准确的，则称该求积公式具有m次代数精度（简称代数精度）,或,),定理1 n+1个节点的求积公式 为插值型求积公式的充要条件是公式 至少具有n次代数精度。,例1 设积分区间a,b为0,2，取时 时,分别用梯形和辛卜生公式,计算其积分结果并与准确值进行比较解:梯形公式和辛卜生的计算结果与准确值比 较如下表所示,f(x)1 x x2 x3 x4 ex 准确值 2 2 2.67 4 6.40 6.389 梯形公式计算值 2 2 4 8 16 8.389 辛卜生公式计算值 2 2 2.67 4 6

7、.67 6.421,从表中可以看出,当f(x)是 时,辛卜生公式比梯形公式更精确,一般说来，代数精度越高，求积公式越精确。梯形公式和中矩形公式具有1次代数精度，辛卜生公式有3次代数精度。下面以梯形公式为例进行验证,取f(x)=1时，,两端相等,取f(x)=x时,取f(x)=x2 时,两端不相等,所以梯形公式只有1次代数精度。,两端相等,构造插值求积公式有如下特点：复杂函数f(x)的积分转化为计算多项式的积分 求积系数Ak只与积分区间及节点xk有关，而与被积函数f(x)无关，可以不管f(x)如何，预先算出Ak的值 n+1个节点的插值求积公式至少具有n次代数精度 求积系数之和 可用此检验计算求积系

8、数的正确性,3 牛顿柯特斯(Newton-Cotes)求积公式 在插值求积公式,中,当所取节点是等距时称为牛顿-柯特斯公式其中 插值多项式 求积系数,这里 是插值基函数。即有,将积分区间a,b 划分为n等分,步长求积节点为 为了计算系数Ak,由于,所以,作变量代换 当 时,有,于是可得,(k=0,1,n),代入插值求积公式(4)有,称为牛顿-柯特斯求积公式,Ck称为柯特斯系数,引进记号,(k=0,1,n),则,容易验证,显然,Ck是不依赖于积分区间a,b以及被积函数f(x)的常数,只要给出n,就可以算出柯特斯系数,譬如当n=1时,当n=2时,4 几个低阶求积公式 在牛顿-柯特斯求积公式中n=1

9、,2,4时，就分别得到下面的梯形公式、辛卜生公式和柯特斯公式。(1)梯形公式 当n=1时，牛顿-柯特斯公式就是梯形公式,定理2（梯形公式的误差）设f(x)在a,b上具有连续的二阶导数，则梯形公式的误差（余项）为,（2）辛卜生公式 当n=2时，牛顿-柯特斯公式就是辛卜生公式（或 称抛物线公式）,定理3（辛卜生公式的误差）设在a,b上具有连续的四阶导数，则辛卜生求积公式的误差为,定理证明从略。,（3）柯特斯公式。当n=4时，牛顿-柯特斯公式为,定理4（柯特斯公式的误差）设在a,b上具有连续的6阶导数，则柯特斯求积公式的误差为,定理的证明从略。,例11 分别用梯形公式、辛卜生公式和柯特斯 公式计算定

10、积分 的近似值(计算结果取5位有效数字),(1)用梯形公式计算,(2)用辛卜生公式,(3)用柯特斯公式计算，系数为,积分的准确值为,可见，三个求积公式的精度逐渐提高。,例12 用辛卜生公式和柯特斯公式计算定积分,的近似值,并估计其误差(计算结果取5位小数),解:辛卜生公式,由于 由辛卜生公式余项,知其误差为,例12 用辛卜生公式和柯特斯公式计算定积分,的近似值,并估计其误差(计算结果取5位小数),解:柯特斯公式,知其误差为,例12 用辛卜生公式和柯特斯公式计算定积分,的近似值,并估计其误差(计算结果取5位小数),该定积分的准确值,这个例子告诉我们，对于同一个积分，当n2时，公式却是精确的，这是

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

12、积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0。返回参数I即定积分值，n为被积函数的调用次数。,例1 求定积分。(1)建立被积函数文件fesin.m。function f=fesin(x)f=exp(-0.5*x).*sin(x+pi/6);(2)调用数值积分函数quad求定积分。S,n=quad(fesin,0,3*pi)S=0.9008n=77,2牛顿柯特斯法基于牛顿柯特斯法，MATLAB给出了quad8函数来求定积分。该函数的调用格式为：I,n=qu

13、ad8(fname,a,b,tol,trace)其中参数的含义和quad函数相似，只是tol的缺省值取10-6。该函数可以更精确地求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。,例2 求定积分。(1)被积函数文件fx.m。function f=fx(x)f=x.*sin(x)./(1+cos(x).*cos(x);(2)调用函数quad8求定积分。I=quad8(fx,0,pi)I=2.4674,例3 分别用quad函数和quad8函数求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。调用函数quad求定积分：format

14、long;fx=inline(exp(-x);I,n=quad(fx,1,2.5,1e-10)I=0.28579444254766n=65,调用函数quad8求定积分：format long;fx=inline(exp(-x);I,n=quad8(fx,1,2.5,1e-10)I=0.28579444254754n=33,3被积函数由一个表格定义在MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y)函数。其中向量X,Y定义函数关系Y=f(X)。例4 用trapz函数计算定积分。命令如下：X=1:0.01:2.5;Y=exp(-X);%生成函数关系数据向量trapz(X

15、,Y)ans=0.28579682416393,1.3 二重定积分的数值求解使用MATLAB提供的dblquad函数就可以直接求出上述二重定积分的数值解。该函数的调用格式为：I=dblquad(f,a,b,c,d,tol,trace)该函数求f(x,y)在a,bc,d区域上的二重定积分。参数tol，trace的用法与函数quad完全相同。,例5 计算二重定积分(1)建立一个函数文件fxy.m：function f=fxy(x,y)global ki;ki=ki+1;%ki用于统计被积函数的调用次数f=exp(-x.2/2).*sin(x.2+y);(2)调用dblquad函数求解。global

16、 ki;ki=0;I=dblquad(fxy,-2,2,-1,1)kiI=1.57449318974494ki=1038,2 数值微分2.1 数值差分与差商2.2 数值微分的实现在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff，其调用格式为：DX=diff(X)：计算向量X的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,n-1。DX=diff(X,n)：计算X的n阶向前差分。例如，diff(X,2)=diff(diff(X)。DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(缺省状态)，按列计算差分；dim=2，按行计算差分。,例6

17、 生成以向量V=1,2,3,4,5,6为基础的范得蒙矩阵，按列进行差分运算。命令如下：V=vander(1:6)DV=diff(V)%计算V的一阶差分,例7 用不同的方法求函数f(x)的数值导数，并在同一个坐标系中做出f(x)的图像。程序如下：f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2);g=inline(3*x.2+4*x-1)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5);x=-3:0.01:3;p=polyfit(x,f(x),5);%用5次多项式p拟合f(x)dp=polyder(p);%对拟合多项式p求导数dpdpx=polyval(dp,x);%求dp在假设点的函数值dx=diff(f(x,3.01)/0.01;%直接对f(x)求数值导数gx=g(x);%求函数f的导函数g在假设点的导数plot(x,dpx,x,dx,.,x,gx,-);%作图,