《数值积分及matlab实现【稻谷书屋】》由会员分享，可在线阅读，更多相关《数值积分及matlab实现【稻谷书屋】（43页珍藏版）》请在金锄头文库上搜索。
1、数值积分与微分数值积分与微分2009.4.221实用内容数值积分和数值微分数值积分和数值微分1引言引言我们知道我们知道,若函数若函数f(x)在区间在区间a,b上连续且其原上连续且其原函数为函数为F(x),则可用则可用Newton-Leibnitz公式公式求得定积分求得定积分求定积分的值求定积分的值,Newton-Leibnitz公式公式无论在理论上无论在理论上还是在解决实际问题上都起了很大作用，但它并不还是在解决实际问题上都起了很大作用，但它并不能完全解决定积分的计算问题，因为积分学涉及的能完全解决定积分的计算问题，因为积分学涉及的实际问题极为广泛，而且极其复杂，在实际计算中实际问题极为广泛，
2、而且极其复杂，在实际计算中经常遇到以下三种情况：经常遇到以下三种情况：2实用内容(1)被积函数f(x)并不一定能够找到用初等函数的有限形式表示的原函数F(x)，例如：Newton-Leibnitz公式就无能为力了(2)还有被积函数f(x)的原函数能用初等函数表示，但表达式太复杂，例如函数并不复杂，但积分后其表达式却很复杂，积分后其原函数F(x)为：3实用内容(3)被积函数被积函数f(x)没有具体的解析表达式没有具体的解析表达式,其函数其函数关系由表格或图形表示。关系由表格或图形表示。对于这些情况对于这些情况,要计算积分的准确值都是十分要计算积分的准确值都是十分困难的。由此可见困难的。由此可见,
3、通过原函数来计算积分有它的通过原函数来计算积分有它的局限性局限性,因而研究一种新的积分方法来解决因而研究一种新的积分方法来解决Newton-Leibniz公式所不能或很难解决的积分问题公式所不能或很难解决的积分问题, , 这时需这时需要用数值解法来建立积分的近似计算方法。要用数值解法来建立积分的近似计算方法。 将积分区间细分将积分区间细分, ,在每一个小区间内用简单函在每一个小区间内用简单函数代替复杂函数进行积分，这就是数值积分的思想，数代替复杂函数进行积分，这就是数值积分的思想，用代数插值多项式去代替被积函数发用代数插值多项式去代替被积函数发f(x)f(x)进行积分进行积分是本章讨论数值积分
4、的主要内容。是本章讨论数值积分的主要内容。4实用内容建立数值积分公式的途径比较多建立数值积分公式的途径比较多,其中最常用的其中最常用的有两种：有两种：（1）由积分中值定理可知，对于连续函数由积分中值定理可知，对于连续函数f(x)，在积，在积分区间分区间a,b内存在一点内存在一点，使得，使得即所求的曲边梯形的面积恰好等于底为即所求的曲边梯形的面积恰好等于底为(b-a)，高为，高为的矩形面积。但是点的矩形面积。但是点的具体位置一般是未知的的具体位置一般是未知的,因而因而的值也是未知的的值也是未知的,称称为为f(x)在区间在区间a,b上上的平均高度。那么只要对平均高度的平均高度。那么只要对平均高度提
5、供一种算法，提供一种算法，相应地就获得一种数值求积方法相应地就获得一种数值求积方法5实用内容三个求积分公式三个求积分公式 梯形公式梯形公式y=f(x)yxaby=f(x)abyx(a+b)/2 中矩形公式中矩形公式按照这种思想，可构造出一些求积分值的近似公式。按照这种思想，可构造出一些求积分值的近似公式。例如例如 分别取分别取 和和则分别得到中矩形公式和梯则分别得到中矩形公式和梯形公式。形公式。y=f(x)abab6实用内容y=f(x)yabSimpson公式公式(a+b)/2f( )的近似值而获得的一种数值积分方法。的近似值而获得的一种数值积分方法。中矩形公式把中矩形公式把a,b的中点处函数
6、值的中点处函数值作为作为平均高度平均高度f( )的近似值而获得的一种数值的近似值而获得的一种数值积分方法。积分方法。ab(a+b)/2在这三个公式中在这三个公式中,梯形公式梯形公式把把f(a),f(b)的加权平均值的加权平均值作为平均高度作为平均高度7实用内容Simpson公式是以函数公式是以函数f(x)在在a,b,(a+b)/2这三点的函这三点的函数值数值f(a),f(b),的加权平均值的加权平均值似值而获得的一种数值积分方法。似值而获得的一种数值积分方法。作为平均高度作为平均高度f( )的的近近（2）先用某个简单函数先用某个简单函数近似逼近近似逼近f(x),用用代替原被积函数代替原被积函数
7、f(x)，即，即以此构造数值算法。从数值计算的角度考虑以此构造数值算法。从数值计算的角度考虑,函数函数应对应对f(x)有充分的逼近程度有充分的逼近程度,并且容易计算其积分。并且容易计算其积分。由于多项式能很好地逼近连续函数由于多项式能很好地逼近连续函数,且又容易计算积且又容易计算积分分,因此将因此将选取为插值多项式选取为插值多项式,这样这样f(x)的积分就的积分就可以用其插值多项式的积分来近似代替可以用其插值多项式的积分来近似代替8实用内容2.2 2.2 插值求积公式插值求积公式设已知设已知f(x)f(x)在节点在节点 有函数值有函数值, ,作作n n次拉格朗日插值多项式次拉格朗日插值多项式
8、式中式中 这里这里 多项式多项式P(x)P(x)易于求积易于求积, ,所以可取所以可取 作为作为 的近似值，即的近似值，即 9实用内容其中其中称为求积系数。给出如下定义称为求积系数。给出如下定义。 定义定义1 1 求积公式求积公式 其系数其系数 时，则称求积公式为插值时，则称求积公式为插值求积公式。求积公式。 (4)(4)10实用内容设插值求积公式的余项为设插值求积公式的余项为 , ,由插值余项定理得由插值余项定理得其中其中 当当f(x)f(x)是次数不高于是次数不高于n n的多项式时，有的多项式时，有 =0,=0,求积公式求积公式(4)(4)能成为准确的等式。由于闭区能成为准确的等式。由于闭
9、区间间a,ba,b上的连续函数可用多项式逼近，所以一个上的连续函数可用多项式逼近，所以一个求积公式能对多大次数的多项式求积公式能对多大次数的多项式f(x)f(x)成为准确等式，成为准确等式，是衡量该公式的精确程度的重要指标，为此给出以是衡量该公式的精确程度的重要指标，为此给出以下定义。下定义。 11实用内容定义定义2（代数精度）（代数精度）设求积公式（设求积公式（4）对于一）对于一 切次数小于等于切次数小于等于m的多项式的多项式( (是准确的，而对于次数为是准确的，而对于次数为m+1m+1的多项式是不准确的，的多项式是不准确的，则称该求积公式具有则称该求积公式具有m m次代数精度（简称代数精度
10、）次代数精度（简称代数精度） 或或)12实用内容定理定理1n+1个节点的求积公式个节点的求积公式 为插值型求积公式的充要条件是公式为插值型求积公式的充要条件是公式 至少具有至少具有n次代数精度。次代数精度。13实用内容例例1设积分区间设积分区间a,b为为0,2，取时，取时 时时, , 分别用梯形和辛卜生公式分别用梯形和辛卜生公式计算其积分结果并与准确值进行比较计算其积分结果并与准确值进行比较解解: :梯形公式和辛卜生的计算结果与准确值比梯形公式和辛卜生的计算结果与准确值比 较如下表所示较如下表所示 14实用内容 f(x)1xx2x3x4ex准确值准确值222.6746.406.389梯形公式计
11、算值梯形公式计算值2248168.389辛卜生公式计算值辛卜生公式计算值222.6746.676.421 从表中可以看出从表中可以看出, ,当当f(x)是是 时时, ,辛卜辛卜生公式比梯形公式更精确生公式比梯形公式更精确 一般说来，代数精度越高，求积公式越精确。一般说来，代数精度越高，求积公式越精确。梯形公式和中矩形公式具有梯形公式和中矩形公式具有1 1次代数精度，辛卜生公次代数精度，辛卜生公式有式有3 3次代数精度次代数精度。下面以梯形公式为例进行验证。下面以梯形公式为例进行验证 15实用内容取取f(x)f(x)=1时，时，两端相等两端相等 取取f(x)=xf(x)=x时时, , 取取f(x
12、)=xf(x)=x2 2 时时, , 两端不相等两端不相等 所以梯形公式只有所以梯形公式只有1 1次代数精度。次代数精度。 两端相等两端相等 16实用内容构造插值求积公式有如下特点：构造插值求积公式有如下特点：(1)复杂函数复杂函数f(x)的积分转化为计算多项式的积分的积分转化为计算多项式的积分(2)求积系数求积系数Ak只与积分区间及节点只与积分区间及节点xk有关，而与被有关，而与被积函数积函数f(x)无关，可以不管无关，可以不管f(x)如何，预先算出如何，预先算出Ak的值的值(3)n+1个节点的插值求积公式至少具有个节点的插值求积公式至少具有n次代数精度次代数精度(4)求积系数之和求积系数之
13、和可用此检验计算求积系数的正确性可用此检验计算求积系数的正确性17实用内容3牛顿牛顿柯特斯柯特斯(Newton-Cotes)求积公式求积公式在插值求积公式在插值求积公式中中,当所取节点是等距时称为牛顿当所取节点是等距时称为牛顿-柯特斯公式柯特斯公式其中其中插值多项式插值多项式求积系数求积系数这里这里是插值基函数。即有是插值基函数。即有18实用内容将积分区间将积分区间 a,b 划分为划分为n等分等分,步长步长求积节点为求积节点为为了计为了计算系数算系数Ak,由于由于,所以所以作变量代换作变量代换当当时时,有有,于是可得于是可得19实用内容(k=0,1,n)代入插值求积公式代入插值求积公式( (4
14、)有有称为牛顿称为牛顿- -柯特斯求积公式柯特斯求积公式,C,Ck k称为柯特斯系数称为柯特斯系数引进记号引进记号(k=0,1,n)则则20实用内容容易验证容易验证显然显然, , C Ck k是不依赖于积分区间是不依赖于积分区间a,ba,b以及被积函数以及被积函数f(x)f(x)的常数的常数, ,只要给出只要给出n,n,就可以算出柯特斯系数就可以算出柯特斯系数, ,譬譬如当如当n=1n=1时时 21实用内容当当n=2=2时时22实用内容4 4 几个低阶求积公式几个低阶求积公式 在牛顿在牛顿- -柯特斯求积公式中柯特斯求积公式中n=1,2,4=1,2,4时，就分别时，就分别得到下面的梯形公式、辛
15、卜生公式和柯特斯公式。得到下面的梯形公式、辛卜生公式和柯特斯公式。(1)(1)梯形公式梯形公式 当当n=1=1时，牛顿时，牛顿- -柯特斯公式就是梯形公式柯特斯公式就是梯形公式定理定理2（梯形公式的误差）设（梯形公式的误差）设f(x)在在 a,b 上具有连上具有连续的二阶导数，则梯形公式的误差（余项）为续的二阶导数，则梯形公式的误差（余项）为23实用内容（2 2） 辛卜生公式辛卜生公式 当当n=2=2时，牛顿时，牛顿- -柯特斯公式就是辛卜生公式（或柯特斯公式就是辛卜生公式（或 称抛物线公式）称抛物线公式）定理定理3 3（辛卜生公式的误差）设在（辛卜生公式的误差）设在a,ba,b上具有连续上具
16、有连续的四阶导数，则辛卜生求积公式的误差为的四阶导数，则辛卜生求积公式的误差为 定理证明从略。定理证明从略。24实用内容（3 3） 柯特斯公式。柯特斯公式。 当当n=4=4时，牛顿时，牛顿- -柯特斯公式为柯特斯公式为定理定理4 4（柯特斯公式的误差）设在（柯特斯公式的误差）设在a,ba,b上具有连上具有连续的续的6 6阶导数，则柯特斯求积公式的误差为阶导数，则柯特斯求积公式的误差为 定理的证明从略。定理的证明从略。 25实用内容例例11分别用梯形公式、辛卜生公式和柯特斯分别用梯形公式、辛卜生公式和柯特斯 公式计算定积分公式计算定积分 的近似值的近似值 ( (计算结果取计算结果取5 5位有效数
17、字位有效数字) ) (1) (1) 用梯形公式计算用梯形公式计算 (2) (2) 用辛卜生公式用辛卜生公式 26实用内容(3) (3) 用柯特斯公式计算，系数为用柯特斯公式计算，系数为积分的准确值为积分的准确值为可见，三个求积公式的精度逐渐提高。可见，三个求积公式的精度逐渐提高。 27实用内容例例12 12 用辛卜生公式和柯特斯公式计算定积分用辛卜生公式和柯特斯公式计算定积分的近似值的近似值, ,并估计其误差并估计其误差( (计算结果取计算结果取5 5位小数位小数) ) 解解: : 辛卜生公式辛卜生公式 由于由于 由辛卜生公式余项由辛卜生公式余项 知其误差为知其误差为 28实用内容例例12用辛
18、卜生公式和用辛卜生公式和柯特斯柯特斯公式计算定积分公式计算定积分的近似值的近似值,并估计其误差并估计其误差(计算结果取计算结果取5位小数位小数)解解:柯特斯公式柯特斯公式知其误差为知其误差为29实用内容例例12用辛卜生公式和用辛卜生公式和柯特斯柯特斯公式计算定积分公式计算定积分的近似值的近似值,并估计其误差并估计其误差(计算结果取计算结果取5位小数位小数) 该定积分的准确值该定积分的准确值,这个例子告诉我这个例子告诉我们，对于同一个积分，当们，对于同一个积分，当n2时，公式却是精确的，时，公式却是精确的，这是由于辛卜生公式具有三次代数精度，这是由于辛卜生公式具有三次代数精度，柯特斯公柯特斯公式
19、具有五次代数精度式具有五次代数精度，它们对被积函数为三次多项，它们对被积函数为三次多项式当然是精确成立的。式当然是精确成立的。30实用内容数值积分基本原理数值积分基本原理求解定积分的数值方法多种多样，如简单的梯形法、辛普生求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson) 法、牛顿柯特斯法、牛顿柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本法等都是经常采用的方法。它们的基本思想都是将整个积分区间思想都是将整个积分区间a,b分成分成n个子区间个子区间xi,xi+1，i=1,2,n，其中，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。这样求
20、定积分问题就分解为求和问题。31实用内容数值积分的实现方法数值积分的实现方法1变步长辛普生法变步长辛普生法基于变步长辛普生法，基于变步长辛普生法，MATLAB给出了给出了quad函数来求函数来求定积分。该函数的调用格式为：定积分。该函数的调用格式为：I,n=quad(fname,a,b,tol,trace)其中其中fname是被积函数名。是被积函数名。a和和b分别是定积分的下限和分别是定积分的下限和上限。上限。tol用来控制积分精度，缺省时取用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非控制是否展现积分过程，若取非0则展现积分过则展现积分过程，取程，取0则
21、不展现，缺省时取则不展现，缺省时取trace=0。返回参数。返回参数I即定即定积分值，积分值，n为被积函数的调用次数。为被积函数的调用次数。32实用内容例例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=7733实用内容2牛顿柯特斯法牛顿柯特斯法基于牛顿柯特斯法，基于牛顿柯特斯法，MATLAB给出了给出了quad8函数来求定积分。该函数的函数来求定积
22、分。该函数的调用格式为：调用格式为：I,n=quad8(fname,a,b,tol,trace)其中参数的含义和其中参数的含义和quad函数相似，只是函数相似，只是tol的缺省值取的缺省值取10-6。 该函数可以该函数可以更精确地求出定积分的值，且一般情况下函数调用的步数明显小于更精确地求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。函数，从而保证能以更高的效率求出所需的定积分值。34实用内容例2求定积分。(1)被积函数文件fx.m。functionf=fx(x)f=x.*sin(x)./(1+cos(x).*cos(x);(2)调用函
23、数quad8求定积分。I=quad8(fx,0,pi)I=2.467435实用内容例例3分别用分别用quad函数和函数和quad8函数求定积分的近似值，函数求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。并在相同的积分精度下，比较函数的调用次数。调用函数调用函数quad求定积分：求定积分：formatlong;fx=inline(exp(-x);I,n=quad(fx,1,2.5,1e-10)I=0.28579444254766n=6536实用内容调用函数调用函数quad8求定积分：求定积分：formatlong;fx=inline(exp(-x);I,n=quad8(fx,1,2
24、.5,1e-10)I=0.28579444254754n=3337实用内容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.2857968241639338实用内容1.3二重定积分的数值求解使用MATLAB提供的
25、dblquad函数就可以直接求出上述二重定积分的数值解。该函数的调用格式为：I=dblquad(f,a,b,c,d,tol,trace)该函数求f(x,y)在a,bc,d区域上的二重定积分。参数tol，trace的用法与函数quad完全相同。39实用内容例例5计算二重定积分计算二重定积分(1)建立一个函数文件建立一个函数文件fxy.m：functionf=fxy(x,y)globalki;ki=ki+1;%ki用于统计被积函数的调用次数用于统计被积函数的调用次数f=exp(-x.2/2).*sin(x.2+y);(2)调用调用dblquad函数求解。函数求解。globalki;ki=0;I=d
26、blquad(fxy,-2,2,-1,1)kiI=1.57449318974494ki=103840实用内容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(dif
27、f(X)。DX=diff(A,n,dim)：计算矩阵：计算矩阵A的的n阶差分，阶差分，dim=1时时(缺省状缺省状态态)，按列计算差分；，按列计算差分；dim=2，按行计算差分。，按行计算差分。41实用内容例例6生成以向量生成以向量V=1,2,3,4,5,6为基础的范得蒙矩为基础的范得蒙矩阵，按列进行差分运算。阵，按列进行差分运算。命令如下：命令如下：V=vander(1:6)DV=diff(V)%计算计算V的一阶差分的一阶差分42实用内容例例7用不同的方法求函数用不同的方法求函数f(x)的数值导数，并在同一个坐标系的数值导数，并在同一个坐标系中做出中做出f(x)的图像。的图像。程序如下：程序
28、如下：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,-);%作图作图43实用内容
