9.数值积分-龙贝格算法+高斯求积公式

4.4龙贝格算法变步长梯形求积法算法简单，但收敛速度较慢，可以利用梯形法算法简单的优点，形成一个新算法，这就是龙贝格求积公式。龙贝格公式又称逐次分半加速法。根据积分区间分成n等份和2n等份时的误差估计式(4.8)可得)(3122nnnTTTI所以积分值的误差大致等于,如果用对进行修正时，与之和比更接近积分真值,所以可以将看成是对误差的一种补偿,因此可得到具有更好效果的式子.nT2)(312nnTT)(312nnTTnT2)(312nnTTnT2nT2)(312nnTTnT2nnnnnTTTTTT3134)(31222(6.9)考察与n等份辛卜生公式之间的关系。将复化梯形公式TnS)()(2)(211bfxfafhTnkkn梯形变步长公式102)(2221nkknnxfhTT代入表达式得TnnknkkkSbfxfxfafhT1010)()(2)(4)(621241(4.10)33nnnSTT故这就是说，用梯形法二分前后两个积分值和作线性组合，结果却得到复化Simpson公式计算得到的积分值。nTnT2nS221()3nnnTTTT再考察Simpson公式。其截断误差与成正比，因此，如果将步长折半，则误差减至，即有4h1611612nnSISI由此可得nnSSI15115162可以验证,上式右端的值其实等于Cn，就是说，用辛卜生公式二等份前后的两个积分值Sn和S2n作线性组合后，可得到柯特斯公式求得的积分值Cn，即有nnnSSC15115162(4.11)用同样的方法，根据柯特斯公式的误差公式，可进一步导出龙贝格公式nnnCCR63163642(4.12)在变步长的过程中运用(4.10)、(4.11)和(4.12)，就能将粗糙的梯形值Tn逐步加工成精度较高的Simpson值Sn、柯特斯值Cn和龙贝格值Rn或者说，将收敛缓慢的梯形值序列Tn加工成收敛迅速的龙贝格值序列Rn，这种加速方法称为龙贝格算法（龙贝格公式）.4.4.3龙贝格求积法算法实现（1）龙贝格求积法计算步骤①用梯形公式计算积分近似值②按变步长梯形公式计算积分近似值将区间逐次分半,令区间长度)()(21bfafabT),2,1,0(2kabhk102)(2221nkknnxfhTT计算)2(kn③按加速公式求加速值322nnnnTTTS1522nnnnSSSC6322nnnnCCCR梯形加速公式：Simpson加速公式：龙贝格求积公式：④精度控制；直到相邻两次积分值nnRR2（其中ε为允许的误差限）则终止计算并取Rn作为积分的近似值，否则将区间再对分，重复②，③，④的计算，直到满足精度要求为止。badxxf)(例4.16用龙贝格算法计算定积分要求相邻两次龙贝格值的偏差不超过解:由题意102d14xxI510214)(,1,0xxfba3)24(21)1()0(211ffT1.351621321)(21212112fTT13118.3)56.2764.3(411.321)()(4121434124ffTT13899.3)()()()(81218785838148ffffTT14094.3)()()()()()()()(16121161516131611169167165163161816ffffffffTT1333.33134121TTS14157.33134242TTS14159.33134484TTS14159.331348168TTS14212.31511516121SSC14159.31511516242SSC14159.31511516484SSC14158.36316364121CCR14159.36316364242CCR由于，于是有00001.012RR14159.3d14102xxI例用Romberg算法计算积分dxxI1021424(),0,11fxabx111[(0)(1)](42)322Tff116,25f21111222TTf1214133STT解(2)由则(1)(3)11633.1253.13333.再计算4211132444TTff3.13118.由此得再计算2424133STT3.14157,1211611515CSS3.14212.84111357288888TTffff3.13899.从而算得4844133STT3.14159,2421611515CSS3.14159,1216416363RCC3.14158.(4)把区间再分半，重复上面的计算步骤，可得163.14094,T83.14159,S23.14159.R210.00001,RR12043.14159.1Idxx由已精确到小数点后五位，故可取上述计算结果参见下表：k区间个数2k2kT12kS22kC32kR012341248163.000003.100003.131183.138992.140943.13333.141573.141593.141593.142123.141593.141593.141583.14159龙贝格Romberg算法：???………………T1=)0(0TT8=)3(0TT4=)2(0TT2=)1(0TS1=)0(1TR1=)0(3TS2=)1(1TC1=)0(2TC2=)1(2TS4=)2(1T4.6高斯（Gauss）型求积公式4.6.1高斯积分问题的提出在前面建立牛顿-柯特斯公式时，为了简化计算，对插值公式中的节点限定为等分的节点，然后再定求积系数，这种方法虽然简便，但求积公式的精度受到限制。我们已经知道，过n+1个节点的插值形求积公式至少具有n次代数精度，我们不仅要问，是否存在具有最高代数精度的求积公式呢？若有，最高代数精度能达到多少呢？让我们先看一个例子：在构造形如的两点公式时,如果限定求积节点,那么所得插值求积公式111100)()()(xfAxfAdxxf（4.13）1,110xx11)1()1()(ffdxxf的代数精度仅为1。但是,如果对式(4.13)中的系数和节点都不加限制，那么就可适当选取和,使所得公式的代数精度。事实上，若要使求积公式(4.13)对函数都准确成立，只要和满足方程组10,AA10,xx10,xx10,AA1m32,,,1)(xxxxf10,xx10,AA03202311300211200110010xAxAxAxAxAxAAA333310010xxAA解之得11)33()33()(ffdxxf代入(4.13)即得(4.14)可以验证，所得公式（4.14）是具有3次代数精度的插值型求积公式。这个例子告诉我们，只要适当选择求积节点，可使插值型求积公式的代数精度达到最高。这就是本节要介绍的高斯求积公式。同理，对于一般的插值求积公式)()(0kbankkxfAdxxf（4.15）只要适当地选取其2n+2个待定参数xk和,就可使它的代数精度达到2n+1次。),,1,0(,nkAk定义4.3若插值求积公式（4.15）具有2n+1次代数精度，则称之为高斯求积公式，并称相应的求积节点为高斯点。),,1,0(,nkxk可以证明，n个节点的高斯求积公式具有最高不超过2n+1次的代数精度，这就是我们所要讨论的具有最高代数精度的插值型求积公式。4.6.2高斯求积公式的构造与应用像构造两点高斯求积公式(4.14)一样,对于插值型求积公式(4.15),分别取用代定系数法来确定参数xk和从而构造n+1个点高斯求积公式。但是，这种做法要解一个包含2n+2个未知数的非线性方程组，其计算工作量是相当大的。一个较简单的方法是：(1)先利用区间a,b上的n+1次正交多项式确定高斯点(2)然后利用高斯点确定求积系数为简单起见,对求积公式(4.15)的求积区间a,b转换成-1,1的形式，作变换12,,,1)(nxxxf),,1,0(,nkAk),,1,0(,nkbaxk),,1,0(nkAk22batabx就可将求积区间a,b变换到-1,1上，这时bababatabdbatabfdxxf)22()22()(11)22(2dtbatabfab11)(2)(dttabdttfba即有)22()(batabft其中插值求积公式节点一经确定，相应的求积系数就确定了，因此关键在于确定节点。定理4.5节点是高斯点的充要条件是：以这些点为零点的多项式),,1,0(,nkxknkkxxx0)()(与任意次数不超过n的P(x)均正交badxxxP0)()((4.16)由定理4.5可知,如能找到满足公式(4.16)的n+1次多项式,则求积公式的高斯点就确定了，进而就可确定相应的高斯求积公式。为此需要引入勒让得（Legendre）多项式及其相关结论)(x定义4.4一个仅以区间-1,1上的高斯点为零点的n+1次多项式称为Legendre多项式。),,1,0(,nkxk定理4.6若是高斯点，则以这些点为根的多项式是最高次幂系数为1的勒让得多项式，即=),,1,0(,nkxk)(x)(~)1(xLn)(x)(~)1(xLn112101)1()!22()!1()(~),()(nnnnknkdxxdnnxLxxx其中从定理可以看出,当n给定,xk就确定了。P122表4-7给出当积分区间是-1,1时，2个点至5个点的高斯求积公式的节点、系数和余项,其中-1,1，需要时可以查用。nxk(n)Ak(n)Rn1022±0.577350313±0.77459675/9=0.555555608/9=0.88888894±0.86113630.3478548±0.33998100.65214525±0.90617990.2369269±0.53846930.478628700.5688889)(157501)6(f3472875)()8(f1237732650)()10(fGauss-Legendre点及系数表)(31f)(1351)4(f例4.17利用三点高斯求积公式计算的近似值。115.1dxx)7745967.0(5555556.0)0(8888889.0)7745967.0(5555556.0)(11fffdxxf由所求公式得05.18888889.07745967.05.15555556.05.111dxx7745967.05.15555556.03997.2高斯求积公式是高精度求积公式，其求积系数，,求积公式也是数值稳定的。),,1,0(,0nkAk解:由表4-7可知，得到三点高斯型求积公式为但它明显的缺点是当n改变时，系数和节点几乎都在改变，因而应用起来十分不便。同时其余项涉及高阶导数，要利用它们来控制精度也十分困难，因此在实际计算中较多采用复合求积的方法。譬如，先把积分区间a,b分成m个等长的小区间,然后在每个小区间上使用同一低阶（如两点的、三点的…）高斯型求积公式算出积分的近似值，将它们相加即得积分的近似值。),,1,0(,1mkxxkkbadxxf)(§4.8数值微分一、中点方法与误差分析.导数值在某点的