数值积分及matlab实现

数值积分与微分2009.4.22数值积分和数值微分1引言我们知道,若函数f(x)在区间[a,b]上连续且其原函数为F(x),则可用Newton-Leibnitz公式baaFbFdxxf)()()(求得定积分求定积分的值,Newton-Leibnitz公式无论在理论上还是在解决实际问题上都起了很大作用，但它并不能完全解决定积分的计算问题，因为积分学涉及的实际问题极为广泛，而且极其复杂，在实际计算中经常遇到以下三种情况：(1)被积函数f(x)并不一定能够找到用初等函数的有限形式表示的原函数F(x)，例如：Newton-Leibnitz公式就无能为力了dxedxxxx10102sin和(2)还有被积函数f(x)的原函数能用初等函数表示，但表达式太复杂，例如函数32)(22xxxf并不复杂，但积分后其表达式却很复杂，积分后其原函数F(x)为：)322ln(2169321633241)(22222xxxxxxxxF(3)被积函数f(x)没有具体的解析表达式,其函数关系由表格或图形表示。对于这些情况,要计算积分的准确值都是十分困难的。由此可见,通过原函数来计算积分有它的局限性,因而研究一种新的积分方法来解决Newton-Leibniz公式所不能或很难解决的积分问题,这时需要用数值解法来建立积分的近似计算方法。将积分区间细分,在每一个小区间内用简单函数代替复杂函数进行积分，这就是数值积分的思想，用代数插值多项式去代替被积函数发f(x)进行积分是本章讨论数值积分的主要内容。建立数值积分公式的途径比较多,其中最常用的有两种：（1）由积分中值定理可知，对于连续函数f(x)，在积分区间[a,b]内存在一点ξ，使得即所求的曲边梯形的面积恰好等于底为(b-a)，高为的矩形面积。但是点ξ的具体位置一般是未知的,因而的值也是未知的,称为f(x)在区间[a,b]上的平均高度。那么只要对平均高度提供一种算法，相应地就获得一种数值求积方法bafabdxxfba,)()()()(f)(f)(f)(f三个求积分公式①梯形公式y=f(x)yxab)()()(21)(bfafabdxxfbay=f(x)abyx(a+b)/2②中矩形公式)2()()(bafabdxxfba按照这种思想，可构造出一些求积分值的近似公式。例如分别取和则分别得到中矩形公式和梯形公式。)(f)2()(baff2)()()(bfaffy=f(x)ababy=f(x)yab③Simpson公式(a+b)/2)()2(4)()(61)(bfbafafabdxxfbaf()的近似值而获得的一种数值积分方法。中矩形公式把[a,b]的中点处函数值作为平均高度f()的近似值而获得的一种数值积分方法。)()(21bfaf)2(bafab(a+b)/2在这三个公式中,梯形公式把f(a),f(b)的加权平均值作为平均高度Simpson公式是以函数f(x)在a,b,(a+b)/2这三点的函数值f(a),f(b),的加权平均值似值而获得的一种数值积分方法。1(()4()())62abfaffb)2(baf作为平均高度f()的近（2）先用某个简单函数近似逼近f(x),用代替原被积函数f(x)，即)(x)(xbabadxxdxxf)()(以此构造数值算法。从数值计算的角度考虑,函数应对f(x)有充分的逼近程度,并且容易计算其积分。由于多项式能很好地逼近连续函数,且又容易计算积分,因此将选取为插值多项式,这样f(x)的积分就可以用其插值多项式的积分来近似代替)(x)(x2.2插值求积公式设已知f(x)在节点有函数值,作n次拉格朗日插值多项式),,1,0(nkxk)(kxfnkkkxlxfxP0)()()()()()()(0kknkjjjkjkxxxxxxxxxl式中)())(()(10nxxxxxxx这里多项式P(x)易于求积,所以可取作为的近似值，即badxxP)(badxxf)(knkkbaknkkbaknkkbabaAxfdxxlxfdxxlxfdxxPdxxf000)()()()()()()(bakkbakkdxxxxxdxxlA)()()()(其中称为求积系数。给出如下定义。定义1求积公式nkkkbaxfAdxxf0)()(其系数时，则称求积公式为插值求积公式。bakkdxxlA)((4)设插值求积公式的余项为,由插值余项定理得)(fRbanbadxxnfdxxPxffR)()!1()()()()()1(ba,其中当f(x)是次数不高于n的多项式时，有=0,求积公式(4)能成为准确的等式。由于闭区间[a,b]上的连续函数可用多项式逼近，所以一个求积公式能对多大次数的多项式f(x)成为准确等式，是衡量该公式的精确程度的重要指标，为此给出以下定义。0)()1(xfn)(fR定义2（代数精度）设求积公式（4）对于一切次数小于等于m的多项式(mxxxxf,,,,1)(2mmxaxaxaaxf2210)(是准确的，而对于次数为m+1的多项式是不准确的，则称该求积公式具有m次代数精度（简称代数精度）或)定理1n+1个节点的求积公式为插值型求积公式的充要条件是公式至少具有n次代数精度。nkkkbaxfAdxxf0)()(例1设积分区间[a,b]为[0,2]，取时时,分别用梯形和辛卜生公式xexxxxxf,,,,,1)(43220)2()0()(ffdxxf20)2()1(4)0(31)(fffdxxf计算其积分结果并与准确值进行比较解:梯形公式和辛卜生的计算结果与准确值比较如下表所示f(x)1xx2x3x4ex准确值222.6746.406.389梯形公式计算值2248168.389辛卜生公式计算值222.6746.676.421从表中可以看出,当f(x)是时,辛卜生公式比梯形公式更精确432,,xxx一般说来，代数精度越高，求积公式越精确。梯形公式和中矩形公式具有1次代数精度，辛卜生公式有3次代数精度。下面以梯形公式为例进行验证babfafabdxxf)()(2)(取f(x)=1时，abababdxba)11(2,1两端相等取f(x)=x时,)(21)(2),(212222abbaababxdxba取f(x)=x2时,baabbabaababdxx))((21)(2),(312222332两端不相等所以梯形公式只有1次代数精度。两端相等构造插值求积公式有如下特点：(1)复杂函数f(x)的积分转化为计算多项式的积分(2)求积系数Ak只与积分区间及节点xk有关，而与被积函数f(x)无关，可以不管f(x)如何，预先算出Ak的值(3)n+1个节点的插值求积公式至少具有n次代数精度(4)求积系数之和可用此检验计算求积系数的正确性abAnkk03牛顿—柯特斯(Newton-Cotes)求积公式在插值求积公式nkkkbabaxfAdxxPxxf0)()(d)(中,当所取节点是等距时称为牛顿-柯特斯公式其中插值多项式求积系数)()()(0nkkkxfxlxPbakkdxxlA)(这里是插值基函数。即有)(xlkdxxxxxdxxlAbankiiikibakk0)(将积分区间[a,b]划分为n等分,步长求积节点为为了计算系数Ak,由于,所以nabh),,1,0(nkkhaxkhikxxik)(nknnkkkkkkhknkxxxxxxxx)!(!)1()())(()(110作变量代换当时,有,于是可得thaxkbax,nt,0dxxxxxdxxlAbankiiikibakk0)(dthhntktkttthknknnnkn0)()1)(1()1()!(!)1(dtitknnkabnnkiikn00))(()!(!)1()(dtitknnkCnnkiiknk00))(()!(!)1((k=0,1…,n)代入插值求积公式(4)有nkkkbaxfCabxxf0)()(d)(称为牛顿-柯特斯求积公式,Ck称为柯特斯系数引进记号kkCabA)((k=0,1…,n)则容易验证10nkkC∵bakkkkdxxlAAabC)(1∴nkbaknkkdxxlabC00)(1111)(10babankkdxabdxxlab显然,Ck是不依赖于积分区间[a,b]以及被积函数f(x)的常数,只要给出n,就可以算出柯特斯系数,譬如当n=1时1011002121)1(!1!011tdtCdttC当n=2时202061)2)(1(!2!02)1(dtttC201132)2(!1!12)1(dtttC200261)1(!0!22)1(dtttC4几个低阶求积公式在牛顿-柯特斯求积公式中n=1,2,4时，就分别得到下面的梯形公式、辛卜生公式和柯特斯公式。(1)梯形公式当n=1时，牛顿-柯特斯公式就是梯形公式)()()(21)(bfafabdxxfba定理2（梯形公式的误差）设f(x)在[a,b]上具有连续的二阶导数，则梯形公式的误差（余项）为),()(12)()(31bafabfR（2）辛卜生公式当n=2时，牛顿-柯特斯公式就是辛卜生公式（或称抛物线公式）)()2(4)()(61)(bfbafafabdxxfba定理3（辛卜生公式的误差）设在[a,b]上具有连续的四阶导数，则辛卜生求积公式的误差为),()(2880)()()4(52bafabfR定理证明从略。（3）柯特斯公式。当n=4时，牛顿-柯特斯公式为)(7)(32)(12)(32)(790)(43210xfxfxfxfxfabdxxfba定理4（柯特斯公式的误差）设在[a,b]上具有连续的6阶导数，则柯特斯求积公式的误差为),()(49458)()6(74bafabfR定理的证明从略。例11分别用梯形公式、辛卜生公式和柯特斯公式计算定积分的近似值(计算结果取5位有效数字)15.0dxx(1)用梯形公式计算4267767.0]170711.0[25.0)]1()5.0([25.01d15.0ffxx(2)用辛卜生公式]/).(.[.d.xx43093403.0]103866.0411707.0[121(3)用柯特斯公式计算，系数为,,,,]17875.03275.012625.0325.07[905.01d15.0xx43096407.0]793326.2939223.1029822.2594975.4[1801积分的准确值为43096441.032d15.02315.0xxx可见，三个求积公式的精度逐渐提高。例12用辛卜生公式和柯特斯公式计算定积分3123d)572(xxxx的近似值,并估计其误差(计算结果取5位小数)解:辛卜生公式322036225941613)(24)(6bfbafafabS由于由辛卜生公式余项572)(23xxxxf0)()4(xfbafabfR,),(2880)()()4(5知其误差为0)(fR例12用辛卜生公式和柯特斯公式计算定积分3123d)572(