1、 第四章第四章 数值积分与数值微分（二）数值积分与数值微分（二）第四节第四节 龙贝格求积公式龙贝格求积公式第五节第五节 高斯求积公式高斯求积公式第六节第六节 数值微分数值微分一、梯形法的递推化一、梯形法的递推化 上一节介绍的复化求积方法对提高精度是行之有效的，但上一节介绍的复化求积方法对提高精度是行之有效的，但在使用求积公式之前必须给出合适的步长，步长取得太大精度在使用求积公式之前必须给出合适的步长，步长取得太大精度难以保证，步长太小则会导致计算量的增加，而事先给出一个难以保证，步长太小则会导致计算量的增加，而事先给出一个恰当的步长又往往是困难的恰当的步长又往往是困难的 实际计算中常常采用变步
2、长的计算方案，即在步长逐次分实际计算中常常采用变步长的计算方案，即在步长逐次分半半(即步长二分即步长二分)的过程中，反复利用复化求积公式进行计算，的过程中，反复利用复化求积公式进行计算，直至所求得的积分值满足精度要求为止直至所求得的积分值满足精度要求为止 设将求积区间设将求积区间a，b分成分成n等分，则一共有等分，则一共有n+1个分点，按个分点，按梯形公式计算积分值梯形公式计算积分值Tn，需要提供，需要提供n+1个函数值如果将求积个函数值如果将求积区间再二分一次，则分点增至区间再二分一次，则分点增至2n+1个，我们将二分前后两个积个，我们将二分前后两个积分值联系起来加以考察分值联系起来加以考察
3、4 龙贝格求积公式龙贝格求积公式 注意到每个子区间注意到每个子区间xk，xk+1经过二分只增加了一个分点经过二分只增加了一个分点xk+1/2(xk+xk+1)/2，用复化梯形公式求得该子区间上的积分，用复化梯形公式求得该子区间上的积分值为值为 101021102110122)12(221)(221)(2)()(4nknnkknnkknkkknhkafhTxfhTxfhxfxfhT二、龙贝格算法二、龙贝格算法).,()(212);,()(12)(222bafhabTIbafhabTIfRnnn 有有：根据复化梯形公式的余项表达式根据复化梯形公式的余项表达式.)(31.41)()(222nnnnn
4、TTTITITIff 整理后可得：整理后可得：，则有，则有假定假定 可见，利用两种步长计算的结果能估计截断误差可见，利用两种步长计算的结果能估计截断误差.若将该截断若将该截断误差加到计算结果中，就可得出误差加到计算结果中，就可得出“改进梯形求积公式改进梯形求积公式”：nnnnnTTTTTT3134)(31222 改进梯形求积公式的右边实际是改进梯形求积公式的右边实际是nnknkkknkknkknkknnnkknnnSbfxfxfafhxfhbfxfafhxfhTTxfhTTT 101121102111102110212)()(2)(4)(6)(2)()(2)(231)(231)(221431)
5、4(31这就是说用梯形法二分前后的两个积分值这就是说用梯形法二分前后的两个积分值Tn与与T2n的线性组的线性组合的结果得到复化辛普森法求积公式合的结果得到复化辛普森法求积公式)1(141144313422nnnnnTTTTS 类似的情况，用辛普森法二分前后的两个积分值类似的情况，用辛普森法二分前后的两个积分值Sn与与S2n的的线性组合的结果可得到复化柯特斯求积公式线性组合的结果可得到复化柯特斯求积公式)2(151151614114422222nnnnnSSSSC 重复同样的手续，用柯特斯法二分前后的两个积分值重复同样的手续，用柯特斯法二分前后的两个积分值Cn与与C2n的线性组合的结果可得到龙贝
6、格的线性组合的结果可得到龙贝格(Romberg)求积公式求积公式)3(631636414114423233nnnnnCCCCR 我们在变步长的过程中运用加速公式我们在变步长的过程中运用加速公式(1)、(2)、(3)，就能，就能将粗糙的梯形值将粗糙的梯形值Tn逐步加工成精度较高的辛普森值逐步加工成精度较高的辛普森值Sn、柯特斯、柯特斯值值Cn和龙贝格值和龙贝格值Rn.龙贝格求积算法可用下表来表示：龙贝格求积算法可用下表来表示：三、理查森三、理查森(Richardson)外推加速法外推加速法 上面讨论说明由梯形公式出发，将区间上面讨论说明由梯形公式出发，将区间a，b逐次二分可逐次二分可提高求积公式
7、的精度，上述加速过程还可继续下去，其理论依提高求积公式的精度，上述加速过程还可继续下去，其理论依据是梯形公式的余项展开，即据是梯形公式的余项展开，即 llllhhhIhThhhIhT242212422121642,)(.)(122nabhbafhabTIn ，22hTTn若记若记Tn=T(h)，当区间，当区间a，b分为分为2n等分时，有等分时，有 ，则，则可见可见I=T(h)的误差为的误差为O(h2)阶阶.若记若记 ，则，则 3)(24)(1hThThT 64162)(6241162411hhIhThhIhT 显然显然T1(h)与积分值与积分值 I 近似的阶为近似的阶为O(h4).这样构造的这
8、样构造的，2)(11hThT就是辛普森公式序列就是辛普森公式序列Sn，S2n，.若令若令 ，)(15121516)(112hThThT 则又可进一步从余项中消去则又可进一步从余项中消去 h4 项，这样构造出的项，这样构造出的 ，其实，其实就是柯特斯公式序列，它与积分值就是柯特斯公式序列，它与积分值 I 的逼近阶为的逼近阶为O(h6).如此继如此继续下去，每加速一次，误差的量级便提高续下去，每加速一次，误差的量级便提高2阶，一般地，若记阶，一般地，若记T0(h)=T(h)，经过，经过m(m=1，2，)次加速后，则有次加速后，则有 )(2hT)(1412144)(11hThThTmmmmmm )2
9、1(141144)(1)1(1)()(0)()(0，次次加加速速值值，可可得得的的表表示示序序列列以以次次后后求求得得的的梯梯形形值值，且且表表示示二二分分设设以以 kTTTmTTkTkmmkmmmkmkkmk.)3()2()1(321数数表表表表出出下下列列三三角角形形数数根根据据公公式式可可以以逐逐行行构构造造、加加速速公公式式，即即可可得得到到、算算法法，若若取取上上式式也也称称为为龙龙贝贝格格求求积积Tm 可以证明，如果可以证明，如果 f(x)充分光滑，那么充分光滑，那么T 数表每一列的元数表每一列的元素及对角线元素均收敛到所求的积分值素及对角线元素均收敛到所求的积分值 I，即，即IT
10、mITkmmkmk )()(lim)(lim，固定固定 机械求积公式机械求积公式 含有含有2n+2个待定参数个待定参数xk、Ak(k0，1，n)当当 xk 为等距节点时得到的插值求积公式的代数精为等距节点时得到的插值求积公式的代数精度至少为度至少为n次，如果适当选取次，如果适当选取 xk(k0，1，n)，有可能使求积公式，有可能使求积公式具有具有 2n+1 次代数精度，这类求积公式称为次代数精度，这类求积公式称为高斯高斯(Gauss)求积公式求积公式.为使问题更具一般性，我们研究带权积分为使问题更具一般性，我们研究带权积分 ，这，这里里r r(x)为权函数，类似机械求积公式，它的求积公式为为权
11、函数，类似机械求积公式，它的求积公式为 Ak(k0，1，n)为不依赖于为不依赖于f(x)的求积系数，的求积系数，xk(k0，1，n)为为求积节点，可适当选取求积节点，可适当选取 xk 及及 Ak(k0，1，n)使积分使积分(1)具有具有2n+1次次代数精度代数精度 bankkkxfAxxf0)(d)(baxxxfId)()(r r)1.5()(d)()(0 nkkkbaxfAxxxfr r5 5 高斯求积公式高斯求积公式一、高斯点一、高斯点 定义定义4 如果求积公式如果求积公式(5.1)具有具有2n+1次代数精度，则称其节点次代数精度，则称其节点 xk(k0，1，n)为高斯点，相应公式为高斯点
12、，相应公式(5.1)称为高斯求积公式称为高斯求积公式.根据定义要使根据定义要使(51)具有具有2n+1次代数精度，只要取次代数精度，只要取f(x)xm，对对m0，1，2n+1，(5.1)精确成立，则得精确成立，则得 当给定权函数当给定权函数r r(x)，求出右端积分，则可由，求出右端积分，则可由(5.2)解得解得 Ak 及及 xk(k0，1，n)bamnkmkknmxxxxA)2.5(.12,1,0d)(0r r 从教材例从教材例5看到求解非线性方程组看到求解非线性方程组(5.2)较复杂，通常较复杂，通常n2就很就很难求解故一般不通过解方程难求解故一般不通过解方程(5.2)求求 xk 及及 A
13、k(k0，1，n)，而从分析高斯点的特性来构造高斯求积公式而从分析高斯点的特性来构造高斯求积公式 定理定理5 插值型求积公式插值型求积公式(5.1)的节点的节点 ax0 xlxnb是是高斯点的充分必要条件是以这些节点为零点的多项式高斯点的充分必要条件是以这些节点为零点的多项式 与任何次数不超过与任何次数不超过n的多项式的多项式P(x)带权带权r r(x)正交，即正交，即)()()(101nnxxxxxxx )5.5(0)(d)()()(1 banxxxxxPr r 定理表明在定理表明在a，b上带权上带权r r(x)的的n+1次正交多项式的零点就是次正交多项式的零点就是求积公式求积公式(5.1)
14、的高斯点，有了求积节点的高斯点，有了求积节点 xk(k0，l，n)，再，再利用利用(5.2)对对m0，l，n 成立，则得到一组关于求积系数成立，则得到一组关于求积系数A0，A1，An 的线性方程解此方程则得的线性方程解此方程则得Ak(k0，1，n).也也可直接由可直接由x0，x1，xn 的插值多项式求出求积系数的插值多项式求出求积系数Ak(k=0,1,n).二、高斯求积公式的余项二、高斯求积公式的余项 利用利用 f(x)在节点在节点xk(k=0，1，n)的的埃尔米特插值埃尔米特插值 H2n+1(x)，即，即于是于是 ，两端乘，两端乘r r(x)，并由，并由a到到b积分，积分，则得则得其中右端第
15、一项积分对其中右端第一项积分对2n+1次多项式精确成立，故次多项式精确成立，故 由于由于 0，故由积分中值定理得，故由积分中值定理得(5.1)的余项为的余项为.,1,0,)()(,)()(1212nkxfxHxfxHkknkkn )()!22()()()(21)22(12xnfxHxfnnn bannnkkknxxxnfxfAIfR.d)()()!22()()(21)22(0r r )()(21xxnr r )8.5(.d)()()!22()(21)22(bannnxxxnffRr r bannbafRxxxHxxxfId)()(d)()(12r rr r三、高斯求积公式的稳定性与收敛性三、高
16、斯求积公式的稳定性与收敛性 定理定理6 高斯求积公式高斯求积公式(5.1)的求积系数的求积系数 Ak(k0，1，n)全是正的全是正的 由本定理及定理由本定理及定理2，则得，则得 推论推论 高斯求积公式高斯求积公式(5.1)是稳定的是稳定的.定理定理7 设设 f(x)C a，b，则高斯求积公式，则高斯求积公式(5.1)是收敛是收敛 的，即的，即 nkbakknxxxfxfA0.d)()()(limr r四、常用的高斯型求积公式四、常用的高斯型求积公式1、高斯、高斯勒让德求积公式勒让德求积公式 勒让德多项式勒让德多项式 是区间是区间-1，1上上权函数权函数 r r(x)=1的正交多项式，若以勒让德
17、多项式的零点为高的正交多项式，若以勒让德多项式的零点为高斯点，则求积公式斯点，则求积公式称为称为高斯高斯勒让德求积公式勒让德求积公式.nnnnnxxnxP)1(dd!21)(2 )9.5()(d)(011 nkkkxfAxxf高斯高斯勒让德求积公式勒让德求积公式(5.9)的余项由的余项由(5.8)得得.)11()(901,)(13511)10.5(1,1)()!22)(32()!1(21)(1,1d)()!22()()4(1)4(1)22(343211121)22(少算一个函数值少算一个函数值还小，且比辛普森公式还小，且比辛普森公式，区间为区间为它比辛普森公式余项它比辛普森公式余项时，有时，有
18、当当三章知识三章知识的勒让德多项式，有第的勒让德多项式，有第是最高项系数为是最高项系数为这里这里 ffRffRnfnnnfRxPxxPnffRnnnnnnn2、高斯、高斯切比雪夫求积公式切比雪夫求积公式 切比雪夫正交多项式切比雪夫正交多项式 是区间是区间-1，1上权函上权函数数 的正交多项式，若选取的正交多项式，若选取 n+1 次多项式的零点次多项式的零点)arccoscos()(xnxTn 211)(xx r r时时，只只要要做做变变换换，而而是是一一般般的的区区间间，当当积积分分区区间间不不是是,11ba 22batabx )11.5(d222d)(1,1,11 tbatabfabxxfb
19、aba，这这时时化化为为可可将将nknkxk,2,1,02212cos ，为高斯点，则求积公式为高斯点，则求积公式称为称为高斯高斯切比雪夫求积公式切比雪夫求积公式.相应的求积系数相应的求积系数)12.5()(d1)(0112 nkkkxfAxxxf.,2,1,0,1d)(11112nknxxlxAkk 其中，其中，lk(x)是关于所选节点的拉格朗日插值基函数是关于所选节点的拉格朗日插值基函数.使用时将使用时将 n+1 个个节点公式改为节点公式改为 n 个节点，于是高斯个节点，于是高斯切比雪夫求积公式可写成切比雪夫求积公式可写成)13.5(212cos,)(d1)(1112 nkkknkxxfn
20、xxxf 公式公式(5.13)的余项由的余项由(5.8)可算得可算得)14.5(.)1,1(,)()!2(22)2(2 nnfnfR 与积分相反，数值微分非常困难。积分描述了一个函数的与积分相反，数值微分非常困难。积分描述了一个函数的整体或宏观性质，而微分则描述一个函数在一点处的斜率，这整体或宏观性质，而微分则描述一个函数在一点处的斜率，这是函数的微观性质。因此积分对函数的形状在小范围内的改变是函数的微观性质。因此积分对函数的形状在小范围内的改变不敏感。而微分却很敏感。一个函数小的变化，容易产生相邻不敏感。而微分却很敏感。一个函数小的变化，容易产生相邻点的斜率的大的改变点的斜率的大的改变。由于
21、微分这个固有的困难，所以应尽可能避免数值微分，由于微分这个固有的困难，所以应尽可能避免数值微分，特别是对实验获得的数据进行微分。在这种情况下，最好用最特别是对实验获得的数据进行微分。在这种情况下，最好用最小二乘曲线拟合这种数据，然后对所得到的多项式进行微分。小二乘曲线拟合这种数据，然后对所得到的多项式进行微分。或用另一种方法，对该数据进行三次样条拟合，然后寻找该样或用另一种方法，对该数据进行三次样条拟合，然后寻找该样条函数的微分。条函数的微分。6 6 数值微分数值微分一、中点方法与误差分析一、中点方法与误差分析 数值微分就是用函数值的线性组合近似函数在某点的导数数值微分就是用函数值的线性组合近
22、似函数在某点的导数 值按导数定义可以简单地用差商近似导数，这样立即得到几值按导数定义可以简单地用差商近似导数，这样立即得到几种数值微分公式种数值微分公式hhafhafafhhafafafhafhafaf2)()()()()()()()()(其中其中h为一增量称为步长后一种数值微分方法称为中点为一增量称为步长后一种数值微分方法称为中点方法、它是前两种方法的算术平均但它的误差阶却由方法、它是前两种方法的算术平均但它的误差阶却由O(h)提提高到高到O(h2)上面所给出的三个公式是很实用的尤其是中点公上面所给出的三个公式是很实用的尤其是中点公式更为常用式更为常用为要利用中点公式为要利用中点公式hhaf
23、hafhG2)()()()(af 计算导数计算导数的近似值，首先须选取合适的步长为此需要进行误差分析的近似值，首先须选取合适的步长为此需要进行误差分析分别将分别将在在 x=a 处做泰勒展开有处做泰勒展开有)(haf )(!5)(!4)(!3)(!2)()()()5(5)4(432afhafhafhafhafhafhaf代入代入G(h)得得 )(!5)(!3)()()5(42afhafhafhG 由此得知，从截断误差的角度看，步长越小，计算结果越由此得知，从截断误差的角度看，步长越小，计算结果越准确准确.且且MhhGaf6)()(2 其中其中.)(maxxfMhax 再考察舍入误差按中点公式计算
24、，当再考察舍入误差按中点公式计算，当h很小时，因很小时，因 f(a+h)与与 f(a-h)很接近，直接相减会造成有效数字的严重损失很接近，直接相减会造成有效数字的严重损失(参看第参看第1章第章第4节节)因此，从舍入误差的角度来看，步长不宜太小因此，从舍入误差的角度来看，步长不宜太小则计算则计算 当当 f(a+h)及及 f(a-h)分别有舍入误差分别有舍入误差 1及及 2时，若令时，若令 21,max )(af 的舍入误差上界为的舍入误差上界为hhaGafaf 2)()()(21步长步长h不宜太大，也不宜太小其最优步长应为不宜太大，也不宜太小其最优步长应为 它表明它表明h越小，舍入误差越小，舍入
25、误差)(af 越大，故它是病态的用中点越大，故它是病态的用中点公式计算公式计算)(af 的误差上界为的误差上界为,6)(2hMhhE 要使误差要使误差E(h)最小，最小，3opt/3Mh 二、插值型的求导公式二、插值型的求导公式 x x0 x1 x2 xn y y0 y1 y2 yn 对于列表函数对于列表函数 y=f(x)：运用插值原理，可以建立插值多项式运用插值原理，可以建立插值多项式 y=Pn(x)作为它的近作为它的近似由于多项式的求导比较容易，我们取似由于多项式的求导比较容易，我们取统称插值型的求导公式统称插值型的求导公式)3.6()()(xPxfn )(xPn)(xf 的近似值，这样建
26、立的数值公式的近似值，这样建立的数值公式的值作为的值作为 依据插值余项定理，求导公式依据插值余项定理，求导公式(63)的余项为的余项为)(dd)!1()()()!1()()()()1(11)1(nnnnnfxnxxnfxPxf式中式中.)()(01 niinxxx 如果我们限定求某个节点如果我们限定求某个节点 xk 上的导数值，那么上面的第二上的导数值，那么上面的第二项因式变为零，这时有余项公式项因式变为零，这时有余项公式)4.6()()!1()()()(1)1(knnknkxnfxPxf 下面我们仅仅考察节点处的导数值为简化讨论，假定所给下面我们仅仅考察节点处的导数值为简化讨论，假定所给的节
27、点是等距的的节点是等距的1两点公式两点公式 设已给出两个节点设已给出两个节点 x0，x1 上的函数值上的函数值 f(x0)，f(x1)，做线性，做线性插值得公式插值得公式)()()(101001011xfxxxxxfxxxxxP 对上式两端求导，记对上式两端求导，记x1 x0=h，有，有)()(1)(101xfxfhxP )()(1)()()(1)(01110101xfxfhxPxfxfhxP 于是有下列求导公式：于是有下列求导公式：)(2)()(1)()(2)()(1)(011010 fhxfxfhxffhxfxfhxf 而利用余项公式而利用余项公式(6.4)知，带余项的两点公式是：知，带余
28、项的两点公式是：2三点公式三点公式 设已给出三个节点设已给出三个节点x0，xl=x0+h，x2=x0+2h上的函数值，上的函数值，做二次插值做二次插值)()()()()()()()()()(2120210121012002010212xfxxxxxxxxxfxxxxxxxxxfxxxxxxxxxP 令令 x=x0+th，则，则 )5.6()()12()()44()()32(21)()()1(21)()2()()2)(1(21)(2100221002xftxftxfththxPxfttxfttxfttthxP 这里撇号这里撇号()表示对变量表示对变量x求导数上式分别取求导数上式分别取t=0，1，
29、2，得，得到三种三点公式：到三种三点公式：.)(3)(4)(21)(;)()(21)(;)()(4)(321)(21022201221002xfxfxfhxPxfxfhxPxfxfxfhxP 而带余项的三点求导公式如下：而带余项的三点求导公式如下：.)(3)(3)(4)(21)()6.6();(6)()(21)(;)(3)()(4)(321)(22102220122100 fhxfxfxfhxffhxfxfhxffhxfxfxfhxf 公式公式(6.6)是我们所熟悉的中点公式在三点公式中，它由于少是我们所熟悉的中点公式在三点公式中，它由于少用了一个函数值用了一个函数值 f(x1)而引人注目而引
30、人注目 用插值多项式用插值多项式Pn(x)作为作为 f(x)的近似函数，还可以建立高阶的近似函数，还可以建立高阶数值微分公式：数值微分公式：,2,1,)()()()(kxPxfkmk例如，将式例如，将式(6.5)再对再对t求导一次，有求导一次，有 )()(2)(1)(210202xfxfxfhthxP 于是有于是有 )()(2)(1)(111212hxfxfhxfhxP 而带余项的二阶三点公式如下：而带余项的二阶三点公式如下：)7.6()(12)()(2)(1)()4(211121 fhhxfxfhxfhxf 三、利用数值积分求导三、利用数值积分求导 微分是积分的逆运算，因此可利用数值积分的方
31、法来计算数微分是积分的逆运算，因此可利用数值积分的方法来计算数值微分设值微分设 f(x)是一个充分光滑的函数，设是一个充分光滑的函数，设habhnkkhaxxfxk ,1,0,)()(则有：则有：)8.6()1,1(d)()()(1111 kkxxkknkxxxfxf 对上式右边积分采用不同的求积公式就可得到不同的数值微分对上式右边积分采用不同的求积公式就可得到不同的数值微分公式例如，对公式例如，对 11d)(kkxxxx 用中矩形公式用中矩形公式，2)(bafabR则得则得),(,)()2(241)(2d)(11311 kkkkkxxxxhxhxxkk 从而得到中点微分公式从而得到中点微分公式)(612)()()(211kkkkfhhxfxfxf ),(,)(90)()(4)(3d)(11)4(51111 kkkkkkkxxxxhxxxhxxkk 若对若对(6.8)右端积分用辛普森求积公式，则有右端积分用辛普森求积公式，则有)()(kkxfx )1,1()()(341111 nkxfxfhmmmkkkkk上式略去余项，并记上式略去余项，并记普森数值微分公式普森数值微分公式的近似值为的近似值为mk，则得到辛，则得到辛这是关于这是关于m0，m1，mn这这 n+1 个未知量的个未知量的 n-1 个方程组，个方程组，可用解方程的方法求解可用解方程的方法求解.
