数值积分及matlab实现

1、数值积分及matlab实现数值积分与微分数值积分与微分2009.4.22数值积分及matlab实现数值积分和数值微分数值积分和数值微分1 引言引言 我们知道我们知道,若函数若函数f(x)在区间在区间a,b上连续且其原上连续且其原函数为函数为F(x),则可用则可用Newton-Leibnitz公式公式baaFbFdxxf)()()(求得定积分求得定积分求定积分的值求定积分的值 , Newton-Leibnitz公式公式 无论在理论上无论在理论上还是在解决实际问题上都起了很大作用，但它并不还是在解决实际问题上都起了很大作用，但它并不能完全解决定积分的计算问题，因为积分学涉及的能完全解决定积分的计算

2、问题，因为积分学涉及的实际问题极为广泛，而且极其复杂，在实际计算中实际问题极为广泛，而且极其复杂，在实际计算中经常遇到以下三种情况：经常遇到以下三种情况：数值积分及matlab实现 (1) 被积函数f(x)并不一定能够找到用初等函数的 有限形式表示的原函数F(x)，例如： Newton-Leibnitz公式就无能为力了dxedxxxx10102sin和(2) 还有被积函数f(x)的原函数能用初等函数表示，但表达式太复杂，例如函数 32)(22xxxf并不复杂，但积分后其表达式却很复杂，积分后其原函数F(x)为： ) 322ln(2169321633241)(22222xxxxxxxxF数值积分

3、及matlab实现(3) 被积函数被积函数f(x)没有具体的解析表达式没有具体的解析表达式, 其函数其函数 关系由表格或图形表示。关系由表格或图形表示。 对于这些情况对于这些情况, 要计算积分的准确值都是十分要计算积分的准确值都是十分困难的。由此可见困难的。由此可见, 通过原函数来计算积分有它的通过原函数来计算积分有它的局限性局限性, 因而研究一种新的积分方法来解决因而研究一种新的积分方法来解决Newton-Leibniz公式所不能或很难解决的积分问题公式所不能或很难解决的积分问题, , 这时需要用数值解法来建立积分的近似计算方法。这时需要用数值解法来建立积分的近似计算方法。 将积分区间细分将

4、积分区间细分, ,在每一个小区间内用简单函在每一个小区间内用简单函数代替复杂函数进行积分，这就是数值积分的思想，数代替复杂函数进行积分，这就是数值积分的思想，用代数插值多项式去代替被积函数发用代数插值多项式去代替被积函数发f(x)f(x)进行积分进行积分是本章讨论数值积分的主要内容。是本章讨论数值积分的主要内容。 数值积分及matlab实现 建立数值积分公式的途径比较多建立数值积分公式的途径比较多, 其中最常用的其中最常用的有两种：有两种：（1）由积分中值定理可知，对于连续函数由积分中值定理可知，对于连续函数f(x)，在，在积分区间积分区间a,b内存在一点内存在一点，使得，使得即所求的曲边梯形

5、的面积恰好等于底为即所求的曲边梯形的面积恰好等于底为(b-a)，高为，高为 的矩形面积。但是点的矩形面积。但是点的具体位置一般是未知的的具体位置一般是未知的, 因而因而 的值也是未知的的值也是未知的, 称称 为为f(x) 在区间在区间a,b上上的平均高度。那么只要对平均高度的平均高度。那么只要对平均高度 提供一种算法，提供一种算法，相应地就获得一种数值求积方法相应地就获得一种数值求积方法bafabdxxfba,)()()()(f)(f)(f)(f数值积分及matlab实现三个求积分公式三个求积分公式 梯形公式梯形公式y=f(x)yxab)()()(21)(bfafabdxxfbay=f(x)a

6、byx(a+b)/2 中矩形公式中矩形公式)2()()(bafabdxxfba按照这种思想，可构造出一些求积分值的近似公式。按照这种思想，可构造出一些求积分值的近似公式。例如例如 分别取分别取 和和则分别得到中矩形公式和梯则分别得到中矩形公式和梯形公式。形公式。)(f)2()(baff2)()()(bfaffy=f(x)abab数值积分及matlab实现y=f(x)yab Simpson公式公式(a+b)/2)()2(4)()(61)(bfbafafabdxxfbaf( )的近似值而获得的一种数值积分方法。的近似值而获得的一种数值积分方法。 中矩形公式把中矩形公式把a,b 的中点处函数值的中点

7、处函数值 作为作为平均高度平均高度f( )的近似值而获得的一种数值积的近似值而获得的一种数值积分方法。分方法。 )()(21bfaf)2(bafab(a+b)/2 在这三个公式中在这三个公式中, 梯形公式梯形公式把把f(a), f(b)的加权平均值的加权平均值 作为平均高度作为平均高度 数值积分及matlab实现Simpson公式是以函数公式是以函数f(x)在在a, b, (a+b)/2这三点的函这三点的函数值数值f(a), f(b), 的加权平均值的加权平均值 似值而获得的一种数值积分方法。似值而获得的一种数值积分方法。1()4()()62abfaffb)2(baf作为平均高度作为平均高度f

8、( )的近的近（2）先用某个简单函数先用某个简单函数 近似逼近近似逼近f(x), 用用 代替原被积函数代替原被积函数f(x)，即，即 )(x)(xbabadxxdxxf)()(以此构造数值算法。从数值计算的角度考虑以此构造数值算法。从数值计算的角度考虑,函数函数 应对应对f(x)有充分的逼近程度有充分的逼近程度,并且容易计算其积分。并且容易计算其积分。由于多项式能很好地逼近连续函数由于多项式能很好地逼近连续函数,且又容易计算积且又容易计算积分分,因此将因此将 选取为插值多项式选取为插值多项式, 这样这样f(x)的积分就的积分就可以用其插值多项式的积分来近似代替可以用其插值多项式的积分来近似代替

9、 )(x)(x数值积分及matlab实现2.2 2.2 插值求积公式插值求积公式 设已知设已知f(x)f(x)在节点在节点 有函数值有函数值, ,作作n n次拉格朗日插值多项式次拉格朗日插值多项式 ), 1 , 0(nkxk)(kxfnkkkxlxfxP0)()()()()()()(0kknkjjjkjkxxxxxxxxxl式中式中 )()()(10nxxxxxxx这里这里 多项式多项式P(x)P(x)易于求积易于求积, ,所以可取所以可取 作为作为 的近似值，即的近似值，即 badxxP)(badxxf)(数值积分及matlab实现knkkbaknkkbaknkkbabaAxfdxxlxfd

10、xxlxfdxxPdxxf 000)()()()()()()(bakkbakkdxxxxxdxxlA)()()()(其中其中 称为求积系数。给出如下定义称为求积系数。给出如下定义。 定义定义1 1 求积公式求积公式 nkkkbaxfAdxxf0)()(其系数其系数 时，则称求积公式为插值时，则称求积公式为插值求积公式。求积公式。 bakkdxxlA)(4)(4)数值积分及matlab实现设插值求积公式的余项为设插值求积公式的余项为 , ,由插值余项定理得由插值余项定理得 )(fRbanbadxxnfdxxPxffR)()!1()()()()() 1(ba,其中其中 当当f(x)f(x)是次数不

11、高于是次数不高于n n的多项式时，有的多项式时，有 =0,=0,求积公式求积公式(4)(4)能成为准确的等式。由于闭区能成为准确的等式。由于闭区间间a,ba,b上的连续函数可用多项式逼近，所以一个上的连续函数可用多项式逼近，所以一个求积公式能对多大次数的多项式求积公式能对多大次数的多项式f(x)f(x)成为准确等式，成为准确等式，是衡量该公式的精确程度的重要指标，为此给出以是衡量该公式的精确程度的重要指标，为此给出以下定义。下定义。 0)()1(xfn)( fR数值积分及matlab实现定义定义2 （代数精度）（代数精度） 设求积公式（设求积公式（4）对于一）对于一 切次数小于等于切次数小于等

12、于m的多项式的多项式( (mxxxxf, 1)(2mmxaxaxaaxf2210)(是准确的，而对于次数为是准确的，而对于次数为m+1m+1的多项式是不准确的，的多项式是不准确的，则称该求积公式具有则称该求积公式具有m m次代数精度（简称代数精度）次代数精度（简称代数精度） 或或)数值积分及matlab实现定理定理1 n+1个节点的求积公式个节点的求积公式 为插值型求积公式的充要条件是公式为插值型求积公式的充要条件是公式 至少具有至少具有n次代数精度。次代数精度。 nkkkbaxfAdxxf0)()(数值积分及matlab实现例例1 设积分区间设积分区间a, b为为0, 2，取时，取时 时时,

13、 , 分别用梯形和辛卜生公式分别用梯形和辛卜生公式 xexxxxxf, 1)(43220)2()0()(ffdxxf20)2() 1 (4)0(31)(fffdxxf计算其积分结果并与准确值进行比较计算其积分结果并与准确值进行比较解解: :梯形公式和辛卜生的计算结果与准确值比梯形公式和辛卜生的计算结果与准确值比 较如下表所示较如下表所示 数值积分及matlab实现 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.67 6

14、.421 从表中可以看出从表中可以看出, ,当当f(x)是是 时时, ,辛辛卜生公式比梯形公式更精确卜生公式比梯形公式更精确 432,xxx 一般说来，代数精度越高，求积公式越精确。一般说来，代数精度越高，求积公式越精确。梯形公式和中矩形公式具有梯形公式和中矩形公式具有1 1次代数精度，辛卜生公次代数精度，辛卜生公式有式有3 3次代数精度次代数精度。下面以梯形公式为例进行验证。下面以梯形公式为例进行验证 数值积分及matlab实现babfafabdxxf)()(2)(取取f(x)f(x)=1时，时， abababdxba) 11 (2,1两端相等两端相等 取取f(x)=xf(x)=x时时, ,

15、 )(21)(2),(212222abbaababxdxba取取f(x)=xf(x)=x2 2 时时, , baabbabaababdxx)(21)(2),(312222332两端不相等两端不相等 所以梯形公式只有所以梯形公式只有1 1次代数精度。次代数精度。 两端相等两端相等 数值积分及matlab实现构造插值求积公式有如下特点：构造插值求积公式有如下特点：复杂函数复杂函数f(x)的积分转化为计算多项式的积分的积分转化为计算多项式的积分 求积系数求积系数Ak只与积分区间及节点只与积分区间及节点xk有关，而与被有关，而与被积函数积函数f(x)无关，可以不管无关，可以不管f(x)如何，预先算出如

16、何，预先算出Ak的值的值 n+1个节点的插值求积公式至少具有个节点的插值求积公式至少具有n次代数精度次代数精度 求积系数之和求积系数之和 可用此检验计算求积系数的正确性可用此检验计算求积系数的正确性 abAnkk0数值积分及matlab实现3 牛顿牛顿柯特斯柯特斯(Newton-Cotes)求积公式求积公式 在插值求积公式在插值求积公式nkkkbabaxfAdxxPxxf0)()(d)(中中,当所取节点是等距时称为牛顿当所取节点是等距时称为牛顿-柯特斯公式柯特斯公式其中其中 插值多项式插值多项式 求积系数求积系数 )()()(0nkkkxfxlxPbakkdxxlA)(这里这里 是插值基函数。

17、即有是插值基函数。即有 )(xlkdxxxxxdxxlAbankiiikibakk 0)(数值积分及matlab实现将积分区间将积分区间 a,b 划分为划分为n等分等分, 步长步长求积节点为求积节点为 为了计为了计算系数算系数Ak, 由于由于 , 所以所以nabh), 1 ,0(nkkhaxkhikxxik)( nknnkkkkkkhknkxxxxxxxx)!( !) 1()()()(110作变量代换作变量代换 当当 时时,有有 ,于是可得于是可得 thaxkbax,nt, 0dxxxxxdxxlAbankiiikibakk 0)(dthhntktkttthknknnnkn0)() 1)(1(

18、) 1()!( !) 1(dtitknnkabnnkiikn 00) )()!( !) 1()(数值积分及matlab实现dtitknnkCnnkiiknk00)()!(!)1( ( k=0,1,n ) 代入插值求积公式代入插值求积公式( (4)有有 nkkkbaxfCabxxf0)()(d)(称为牛顿称为牛顿- -柯特斯求积公式柯特斯求积公式,C,Ck k称为柯特斯系数称为柯特斯系数引进记号引进记号kkCabA)( ( k=0,1,n ) 则则数值积分及matlab实现容易验证容易验证 10nkkC bakkkkdxxlAAabC)(1 nkbaknkkdxxlabC00)(1111)(10

19、 babankkdxabdxxlab显然显然, , C Ck k是不依赖于积分区间是不依赖于积分区间a,ba,b以及被积函数以及被积函数f(x)f(x)的常数的常数, ,只要给出只要给出n,n,就可以算出柯特斯系数就可以算出柯特斯系数, ,譬譬如当如当n=1n=1时时 1011002121) 1(! 1! 011tdtCdttC数值积分及matlab实现当当n=2=2时时 202061)2)(1(!2!02)1(dtttC201132)2(! 1! 12) 1(dtttC200261)1(!0!22)1(dtttC数值积分及matlab实现4 4 几个低阶求积公式几个低阶求积公式 在牛顿在牛顿

20、- -柯特斯求积公式中柯特斯求积公式中n=1,2,4=1,2,4时，就分别时，就分别得到下面的梯形公式、辛卜生公式和柯特斯公式。得到下面的梯形公式、辛卜生公式和柯特斯公式。(1)(1) 梯形公式梯形公式 当当n=1=1时，牛顿时，牛顿- -柯特斯公式就是梯形公式柯特斯公式就是梯形公式 )()()(21)(bfafabdxxfba定理定理2 （梯形公式的误差）设（梯形公式的误差）设f(x)在在 a,b 上具有连上具有连续的二阶导数，则梯形公式的误差（余项）为续的二阶导数，则梯形公式的误差（余项）为),()(12)()(31bafabfR 数值积分及matlab实现（2 2） 辛卜生公式辛卜生公式

21、 当当n=2=2时，牛顿时，牛顿- -柯特斯公式就是辛卜生公式（或柯特斯公式就是辛卜生公式（或 称抛物线公式）称抛物线公式） )()2(4)()(61)(bfbafafabdxxfba定理定理3 3（辛卜生公式的误差）设在（辛卜生公式的误差）设在a,ba,b上具有连续上具有连续的四阶导数，则辛卜生求积公式的误差为的四阶导数，则辛卜生求积公式的误差为 ),()(2880)()()4(52bafabfR定理证明从略。定理证明从略。 数值积分及matlab实现（3 3） 柯特斯公式。柯特斯公式。 当当n=4=4时，牛顿时，牛顿- -柯特斯公式为柯特斯公式为 )(7)(32)(12)(32)(790)

22、(43210 xfxfxfxfxfabdxxfba定理定理4 4（柯特斯公式的误差）设在（柯特斯公式的误差）设在a,ba,b上具有连上具有连续的续的6 6阶导数，则柯特斯求积公式的误差为阶导数，则柯特斯求积公式的误差为 ),()(49458)()6(74bafabfR定理的证明从略。定理的证明从略。 数值积分及matlab实现例例11 分别用梯形公式、辛卜生公式和柯特斯分别用梯形公式、辛卜生公式和柯特斯 公式计算定积分公式计算定积分 的近似值的近似值 ( (计算结果取计算结果取5 5位有效数字位有效数字) ) 15 . 0dxx(1) (1) 用梯形公式计算用梯形公式计算 4267767. 0

23、 170711. 025. 0)1 () 5 . 0(25 . 01d15 . 0ffxx(2) (2) 用辛卜生公式用辛卜生公式 /).(.d.xx43093403. 0 103866. 0411707. 0121数值积分及matlab实现(3) (3) 用柯特斯公式计算，系数为用柯特斯公式计算，系数为 , 17875. 03275. 012625. 0325 . 07 905 . 01d15 . 0 xx43096407. 0793326.2939223.1029822.2594975. 41801积分的准确值为积分的准确值为 43096441. 032d15 . 02315 . 0 xx

24、x可见，三个求积公式的精度逐渐提高。可见，三个求积公式的精度逐渐提高。 数值积分及matlab实现例例12 12 用辛卜生公式和柯特斯公式计算定积分用辛卜生公式和柯特斯公式计算定积分3123d)572(xxxx的近似值的近似值, ,并估计其误差并估计其误差( (计算结果取计算结果取5 5位小数位小数) ) 解解: : 辛卜生公式辛卜生公式 322036225941613)(24)(6bfbafafabS由于由于 由辛卜生公式余项由辛卜生公式余项 572)(23xxxxf0)()4(xfbafabfR,),(2880)()()4(5知其误差为知其误差为 0)(fR数值积分及matlab实现例例1

25、2 用辛卜生公式和用辛卜生公式和柯特斯柯特斯公式计算定积分公式计算定积分3123d)572(xxxx的近似值的近似值,并估计其误差并估计其误差(计算结果取计算结果取5位小数位小数) 解解:柯特斯公式柯特斯公式 知其误差为知其误差为 0)(fR322097812532912835327451) 3 (7) 5 . 2(32) 2(12) 5 . 1 (32) 1 (79013fffffC数值积分及matlab实现例例12 用辛卜生公式和用辛卜生公式和柯特斯柯特斯公式计算定积分公式计算定积分3123d)572(xxxx的近似值的近似值,并估计其误差并估计其误差(计算结果取计算结果取5位小数位小数)

26、 该定积分的准确值该定积分的准确值 ,这个例子告诉我这个例子告诉我们，对于同一个积分，当们，对于同一个积分，当n2时，公式却是精确的，时，公式却是精确的，这是由于辛卜生公式具有三次代数精度，这是由于辛卜生公式具有三次代数精度，柯特斯公柯特斯公式具有五次代数精度式具有五次代数精度，它们对被积函数为三次多项，它们对被积函数为三次多项式当然是精确成立的。式当然是精确成立的。 3220I数值积分及matlab实现数值积分基本原理数值积分基本原理 求解定积分的数值方法多种多样，如简单的梯形法、辛普生求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson) 法、牛顿柯特斯法、牛顿柯特斯(Ne

27、wton-Cotes)法等都是经常采用的方法。法等都是经常采用的方法。它们的基本思想都是将整个积分区间它们的基本思想都是将整个积分区间a,b分成分成n个子区间个子区间xi,xi+1，i=1,2,n，其中，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。这样求定积分问题就分解为求和问题。数值积分及matlab实现数值积分的实现方法数值积分的实现方法1变步长辛普生法变步长辛普生法基于变步长辛普生法，基于变步长辛普生法，MATLAB给出了给出了quad函数来求函数来求定积分。该函数的调用格式为：定积分。该函数的调用格式为： I,n=quad(fname,a,b,tol,trace)其中

28、其中fname是被积函数名。是被积函数名。a和和b分别是定积分的下限和分别是定积分的下限和上限。上限。tol用来控制积分精度，缺省时取用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非控制是否展现积分过程，若取非0则展现积分过则展现积分过程，取程，取0则不展现，缺省时取则不展现，缺省时取trace=0。返回参数。返回参数I即定即定积分值，积分值，n为被积函数的调用次数。为被积函数的调用次数。数值积分及matlab实现 例例1 求定积分。求定积分。 (1) 建立被积函数文件建立被积函数文件fesin.m。function f=fesin(x)f=exp(-0.5

29、*x).*sin(x+pi/6); (2) 调用数值积分函数调用数值积分函数quad求定积分。求定积分。S,n=quad(fesin,0,3*pi)S = 0.9008n = 77数值积分及matlab实现2牛顿柯特斯法牛顿柯特斯法基于牛顿柯特斯法，基于牛顿柯特斯法，MATLAB给出了给出了quad8函数来求定积分。该函数的函数来求定积分。该函数的调用格式为：调用格式为：I,n=quad8(fname,a,b,tol,trace)其中参数的含义和其中参数的含义和quad函数相似，只是函数相似，只是tol的缺省值取的缺省值取10-6。 该函数可以该函数可以更精确地求出定积分的值，且一般情况下函数

30、调用的步数明显小于更精确地求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。函数，从而保证能以更高的效率求出所需的定积分值。数值积分及matlab实现例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数值积分及matlab实现例例3 分别用分别用quad函数和函数和quad8函数求定积分的近似值，函数求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。并在

31、相同的积分精度下，比较函数的调用次数。调用函数调用函数quad求定积分：求定积分：format long;fx=inline(exp(-x);I,n=quad(fx,1,2.5,1e-10)I = 0.28579444254766n = 65数值积分及matlab实现 调用函数调用函数quad8求定积分：求定积分：format long;fx=inline(exp(-x);I,n=quad8(fx,1,2.5,1e-10)I = 0.28579444254754n = 33数值积分及matlab实现3被积函数由一个表格定义被积函数由一个表格定义在在MATLAB中，对由表格形式定义的函数关系的求

32、定积分问题中，对由表格形式定义的函数关系的求定积分问题用用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数值积分及matlab实现1.3 二重定积分的数值求解使用MATLAB提供的dblquad函数就可以直接求出上述二重定积分的数值解。该函数的调用格式为：I=dblquad(f,a,b,c,d,tol,trace)该函数求f(x,y)在a,bc,d区域上的二重定积分。参数tol，trace的用法与函数quad完全相同。数值积分及matlab实现例例5 计算二重定积分计算二重定积分(1) 建立一个函数文件建立一个函数文件fxy.m：function f=fxy(x,y)global ki;ki=ki+1; %ki用于统计被积函数的