1、数值积分与微分,2009.4.22,数值积分和数值微分,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)进行积分是本章讨论数值积分的主要内容。,建立数值积分公式的途径比较多, 其中最
3、常用的有两种： （1）由积分中值定理可知，对于连续函数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, Si
4、mpson公式,(a+b)/2,f()的近似值而获得的一种数值积分方法。中矩形公式把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
7、4 8 16 8.389辛卜生公式计算值 2 2 2.67 4 6.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个节点的插值
8、求积公式至少具有n次代数精度求积系数之和 可用此检验计算求积系数的正确性,3 牛顿柯特斯(Newton-Cotes)求积公式在插值求积公式,中,当所取节点是等距时称为牛顿-柯特斯公式 其中 插值多项式 求积系数,这里 是插值基函数。即有,将积分区间a,b 划分为n等分, 步长 求积节点为 为了计 算系数Ak, 由于 , 所以,作变量代换 当 时,有 ,于是可得,( k=0,1,n ),代入插值求积公式(4)有,称为牛顿-柯特斯求积公式,Ck称为柯特斯系数,引进记号,( k=0,1,n ),则,容易验证,显然, Ck是不依赖于积分区间a,b以及被积函数f(x)的常数,只要给出n,就可以算出柯特斯
9、系数,譬如当n=1时,当n=2时,4 几个低阶求积公式在牛顿-柯特斯求积公式中n=1,2,4时，就分别 得到下面的梯形公式、辛卜生公式和柯特斯公式。(1) 梯形公式当n=1时，牛顿-柯特斯公式就是梯形公式,定理2 （梯形公式的误差）设f(x)在a,b上具有连续的二阶导数，则梯形公式的误差（余项）为,（2） 辛卜生公式当n=2时，牛顿-柯特斯公式就是辛卜生公式（或称抛物线公式）,定理3（辛卜生公式的误差）设在a,b上具有连续的四阶导数，则辛卜生求积公式的误差为,定理证明从略。,（3） 柯特斯公式。当n=4时，牛顿-柯特斯公式为,定理4（柯特斯公式的误差）设在a,b上具有连续的6阶导数，则柯特斯求
10、积公式的误差为,定理的证明从略。,例11 分别用梯形公式、辛卜生公式和柯特斯公式计算定积分的近似值 (计算结果取5位有效数字),(1) 用梯形公式计算,(2) 用辛卜生公式,(3) 用柯特斯公式计算，系数为,积分的准确值为,可见，三个求积公式的精度逐渐提高。,例12 用辛卜生公式和柯特斯公式计算定积分,的近似值,并估计其误差(计算结果取5位小数),解: 辛卜生公式,由于 由辛卜生公式余项,知其误差为,例12 用辛卜生公式和柯特斯公式计算定积分,的近似值,并估计其误差(计算结果取5位小数),解:柯特斯公式,知其误差为,例12 用辛卜生公式和柯特斯公式计算定积分,的近似值,并估计其误差(计算结果取
11、5位小数),该定积分的准确值 ,这个例子告诉我们，对于同一个积分，当n2时，公式却是精确的，这是由于辛卜生公式具有三次代数精度，柯特斯公式具有五次代数精度，它们对被积函数为三次多项式当然是精确成立的。,数值积分基本原理求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)法、牛顿柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间a,b分成n个子区间xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。,数值积分的实现方法 1变步长辛普生法 基于变步长辛普生法，MATLAB给出了quad函数来求定积分
12、。该函数的调用格式为：I,n=quad(fname,a,b,tol,trace) 其中fname是被积函数名。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.9008 n
13、 =77,2牛顿柯特斯法 基于牛顿柯特斯法，MATLAB给出了quad8函数来求定积分。该函数的调用格式为： I,n=quad8(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
14、,例3 分别用quad函数和quad8函数求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。 调用函数quad求定积分： format long; fx=inline(exp(-x); I,n=quad(fx,1,2.5,1e-10) I =0.28579444254766 n =65,调用函数quad8求定积分： format long; fx=inline(exp(-x); I,n=quad8(fx,1,2.5,1e-10) I =0.28579444254754 n =33,3被积函数由一个表格定义 在MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y
15、)函数。其中向量X,Y定义函数关系Y=f(X)。 例4 用trapz函数计算定积分。 命令如下： X=1:0.01:2.5; Y=exp(-X); %生成函数关系数据向量 trapz(X,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： f
16、unction f=fxy(x,y) global ki; ki=ki+1; %ki用于统计被积函数的调用次数 f=exp(-x.2/2).*sin(x.2+y); (2) 调用dblquad函数求解。 global ki;ki=0; I=dblquad(fxy,-2,2,-1,1) ki I =1.57449318974494 ki =1038,2 数值微分 2.1 数值差分与差商 2.2 数值微分的实现 在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff，其调用格式为： DX=diff(X)：计算向量X的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,
17、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 生成以向量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求导数dp dpx=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,-); %作图,
