第三章数值积分

1、一、数值求积的基本思想一、数值求积的基本思想)()()(aFbFdxxfba 积分积分 只要找到被积函数只要找到被积函数 f (x)原函数原函数F(x)，便有，便有牛顿牛顿莱布尼兹莱布尼兹(NewtonLeibniz)公式公式 baxxfId)(实际困难实际困难：大量的被积函数（：大量的被积函数（ , sin x2 等）等）, 找不到用初等函找不到用初等函数表示的原函数数表示的原函数；另外；另外, f (x)是（测量或数值计算出的）一张数是（测量或数值计算出的）一张数据表时，据表时，牛顿牛顿莱布尼兹公式莱布尼兹公式也也不能直接运用不能直接运用。xxsin 积分中值定理：在积分中值定理：在a,

2、b内存在一点内存在一点 ，有，有 f( )成立。成立。 )(d)(abxxfba 1 引言引言第三章第三章 数值积分数值积分(Numerical Integration) 就是说就是说, 底为底为b- -a 而高为而高为f( )的的矩形面积矩形面积恰恰等于所求等于所求曲边梯形的面积曲边梯形的面积 .问题问题 在于点在于点的具体位置一般是不知道的，因而的具体位置一般是不知道的，因而 难以准确算出难以准确算出 f( )的值的值我们将我们将f ( )称为区间称为区间a, b上的平均高度这样上的平均高度这样,只要对只要对平均高度平均高度f( )提供一种算法提供一种算法,相应地便获得一种数值求积方法相应

3、地便获得一种数值求积方法 如果用两端点的如果用两端点的“高度高度”f(a)与与f(b)的算术平均作为平均高度的算术平均作为平均高度f ( ) 的近似值，这样导出的求积公式的近似值，这样导出的求积公式 : 便是我们所熟悉的便是我们所熟悉的梯形公式梯形公式(trapezoidal rule). )()(2bfafabT 2)(bafabR2bac 而如果改用区间中点而如果改用区间中点 的的“高度高度”f (c)近似地取代平近似地取代平均高度均高度f ( )，则又可导出所谓，则又可导出所谓中矩形公式中矩形公式(今后简称矩形公式今后简称矩形公式)：(1.1)(1.2) 更一般地，我们可以在区间更一般地

4、，我们可以在区间a，b上适当选取某些节点上适当选取某些节点 xk ，然后然后用用 f (xk )加权平均得到平均高度加权平均得到平均高度 f ()的近似值的近似值，这样构造出的，这样构造出的求积公式具有下列形式求积公式具有下列形式式中式中 xk 称为称为求积节点求积节点；Ak 称为称为求积系数求积系数，亦称为伴随节点，亦称为伴随节点 xk 的的权权权权Ak 仅仅与节点仅仅与节点xk 的选取有关，而不依赖于被积函数的选取有关，而不依赖于被积函数 f(x)的的具体形式具体形式 ban0kkkxfAdxxf)()(使积分公式具有通用性使积分公式具有通用性 这类数值积分方法通常称作这类数值积分方法通常

5、称作机械求积机械求积， 其特点是将积分求其特点是将积分求值问题归结为函数值的计算，这就避开了牛顿值问题归结为函数值的计算，这就避开了牛顿莱布尼兹公式莱布尼兹公式需要寻求原函数的困难需要寻求原函数的困难(1.3)二、二、代数精度的概念代数精度的概念 数值求积方法是近似方法，为要保证精度，我们自然希望求积数值求积方法是近似方法，为要保证精度，我们自然希望求积公式能对公式能对“尽可能多尽可能多”的函数准确地成立，这就提出了所谓代数精的函数准确地成立，这就提出了所谓代数精度的概念度的概念 定义定义 1 如果某个求积公式对于次数如果某个求积公式对于次数m的多项式均能准确地成的多项式均能准确地成立，但对于

6、立，但对于m+1次多项式就不一定准确，则称该求积公式具有次多项式就不一定准确，则称该求积公式具有m次代次代数精度数精度 一般地，欲使求积公式一般地，欲使求积公式 具有具有m次代数次代数精度，只要令它对于精度，只要令它对于f (x) = 1，x，xm 都能准确成立，这就要求都能准确成立，这就要求 bankkkxfAxxf0)(d)( . )(11;)(21;1122mmmkkkkkabmxAabxAabA例例1： 考察其代数精度。考察其代数精度。 f(x)abf(a)f(b)梯形公式梯形公式解：解：逐次检查公式是否精确成立逐次检查公式是否精确成立代入代入 P0 = 1： baabdx111 2

7、ab=代入代入 P1 = x ：=代入代入 P2 = x2 ： 222abbadxx 2baab 3233abbadxx 222baab 代数精度代数精度 = 1)()(2)(bfafabdxxfba 例例2 试构造形如试构造形如 f(x)dx A0f(0)+ A1f(h)+ A2f(2h) 的数值的数值求积公式求积公式,使其代数精度尽可能高使其代数精度尽可能高,并指出其代数精度的阶数并指出其代数精度的阶数.3h0解解: 令公式对令公式对 f(x)=1,x, x2 均准确成立均准确成立,则有则有3h=A0+ A1+ A2h2=0 + A1h+ A22h9h3=0 + A1h2+ A24h229

8、故求积公式的形式为故求积公式的形式为解之得解之得 A0= h, A1=0, A2= h. 94 34 f(x)dx f(0) + f(2h)3h49h43h0由公式的构造知由公式的构造知,公式公式至少至少具有具有2次代数精度次代数精度; 而当而当f(x)=x3时时,公式的左边公式的左边= h4, 右边右边=18h4, 公式的左边公式的左边 右边右边,说明说明此公式对此公式对 f(x)=x3不能准确成立不能准确成立.因此因此,公式只具有公式只具有2次代数次代数精度精度.814三、三、求积公式的收敛性与稳定性求积公式的收敛性与稳定性 定理定理3表明，只要求积系数表明，只要求积系数Ak0 (k0,1

9、,n)，就能保证，就能保证计算的稳定性计算的稳定性 定义定义2 在求积公式在求积公式 中，若中，若 其中其中 ，则称求积公式是收敛的，则称求积公式是收敛的 由于计算由于计算 f (xk)可能有误差可能有误差,实际得到实际得到 定义定义3 对任给对任给 e e 0，若，若 (k=0,1, ,n), 就有就有 , 则称求积公式是稳定的则称求积公式是稳定的. bankkkxfAxxf0)(d)(e e |)(|00nkkknkkkfAxfA)(max11 iinixxhkkfxf)(0 ，只要，只要 .)(,kkkkfxff 即即 bankkkhnxxfxfAd)()(lim00 定理定理3 若求积

10、公式若求积公式(13)中系数中系数Ak0 (k0，1，n)，则此求积公式是稳定的则此求积公式是稳定的近似近似计算计算 badxxfI)(思思路路利用利用插值多项式插值多项式 则积分易算。则积分易算。)()(xfxPn 在在a, b上取上取 a x0 x1 xn b，做，做 f 的的 n 次插值次插值多项式多项式 ，即得到，即得到 nkkknxlxfxL0)()()( babaknkkdxxlxfdxxf)()()(0Ak bakjxxxxkdxAjkj)()(由由 决定，决定，与与 无关。无关。节点节点 f (x)插值型积分公式插值型积分公式Interpolatory quadratureba

11、nkkxnbanbanbankkkdxxxnfdxxRdxxLxfxfAdxxffR0)1(0)()!1()()()()()()(误差误差bandxxP)(2 插值型的求积公式插值型的求积公式与与Newton-Cotes 公式公式关键是关键是f(x)如果求积公式是插值型的如果求积公式是插值型的, 按余项式按余项式, 对于次数对于次数 n的多项式的多项式 f (x)，其余项其余项R f 等于等于0，因而这时求积公式至少具有，因而这时求积公式至少具有n次代数精度次代数精度定理定理1：形如形如 的求积公式至少有的求积公式至少有 n 次代数精度次代数精度 该该公式为公式为插值型插值型（即：（即： ）

12、nkkkxfA0)( bakkdxxlA)(为便于计算，一般取为便于计算，一般取等距离节点等距离节点得到近似公式得到近似公式:一、一、Newton-Cotes 公式公式2、 把把a, b二等分，作二等分，作2次插值，有次插值，有)()(4)( )(66bffafdxxfbaabba此公式称为此公式称为辛普森（辛普森（Simpson）公式）公式。badxxL)(21、 对于对于a, b上上1次插值，有次插值，有)()()(1bfafxLabaxbabx )()()(2221bfafdxxfAAabbaab 此即此即梯形公式梯形公式。 节点节点等距分布等距分布：ninabhhiaxi,., 1,

13、0, dxxxxxAnxxijjiji 0)()( nijinnijdtjtininabdthhjihjt00)()!( !) 1)()()(令令htax Cotes系数系数)(niC注：注：Cotes 系数仅取决于系数仅取决于 n 和和 i，可查表得到。与可查表得到。与 f (x) 及区及区间间a, b均无关。均无关。 3、 把把a, b n 等分，用插值等分，用插值Ln(x)近似近似 f(x)积分，有积分，有当当n=4时时, 牛顿牛顿-柯特斯公式特别称作柯特斯公式特别称作柯特斯公式柯特斯公式,其形式为其形式为 )(7)(32)(12)(32)(79043210 xfxfxfxfxfabC

14、21,21)1(1)1(0 CCn = 1:)()(2)(bfafabdxxfba Trapezoidal RuledxbxaxffRbax)(!2)( /* 令令 x = a+th, h = b a, 用中用中值定理值定理 */1, , )(1213abhbafh 代数精度代数精度 = 1n = 2:61,32,61)2(2)2(1)2(0 CCC)()(4)(6)(2bffafabdxxfbaba Simpsons Rule代数精度代数精度 = 32,),(,)(901)4(5abhbafhfR n = 3: Simpsons 3/8-Rule, 代数精度代数精度 = 3,)(803)5(

15、5 fhfR n = 4: Cotes Rule, 代数精度代数精度 = 5,)(9458)6(7 fhfR n 为为偶数阶偶数阶的的Newton-Cotes 公式至少有公式至少有 n+1 次代数精度。次代数精度。二、几种低阶求积公式的余项二、几种低阶求积公式的余项三、三、偶阶求积公式的代数精度偶阶求积公式的代数精度 作为插值型的求积公式，作为插值型的求积公式，n 阶的牛顿阶的牛顿-柯特斯公式至柯特斯公式至少具有少具有n 次的插值精度（定理次的插值精度（定理1）。实际的代数精度还可）。实际的代数精度还可进一步提高，一般地，可以证明下述定理进一步提高，一般地，可以证明下述定理: 定理定理 2 当

16、阶当阶 n 为偶数时，牛顿为偶数时，牛顿-柯特斯公式柯特斯公式至少有至少有 n+1 次代数精度次代数精度 . nkknknxfCabI0)()()(注：注：由公式知，当由公式知，当n8时，柯特斯系数出现负值，这时时，柯特斯系数出现负值，这时，初始数据误差将会引起计算结果误差增大，即计算不，初始数据误差将会引起计算结果误差增大，即计算不稳定。因此，实际计算不用稳定。因此，实际计算不用n8的牛顿的牛顿-柯特斯公式柯特斯公式 .估计估计截断误差截断误差为为解解 用用梯形公式梯形公式计算计算:=2.1835估计估计截断误差截断误差为为=0.6796用用Simpson公式公式计算：计算：=2. 0263

17、例例3 试分别使用梯形公式和试分别使用梯形公式和Simpson公式计算积分公式计算积分 的近的近似值，并估计截断误差似值，并估计截断误差.=198.4306890. 0)(max2880)12()4(2152 xfRx=0.068903 3 复化求积公式复化求积公式高次插值有高次插值有Runge 现象现象，故采用分段低次插值，故采用分段低次插值 分段低次合成的分段低次合成的 Newton-Cotes 复合复合求积公式。求积公式。一、复化梯形公式一、复化梯形公式:),., 0(,nkhkaxnabhk 在每个在每个 上用梯形公式：上用梯形公式：,1kkxx 11)()(2)(2nkkbfxfaf

18、h bankkkxfxfhdxxf11)()(2)(=Tn),(),()(12)()(12)(1221213bafabhnfabhfhfRnkknkk /*中值定理中值定理*/nkxfxfxxdxxfkkxxkkkk,., 1,)()(2)(111 二、复化辛普森公式二、复化辛普森公式:),., 0(,nkhkaxnabhk )()(4)(6)(1211 kkkxxxfxfxfhdxxfkkkx21 kx1 kx44444 )()(2)(4)(6)(1010121 nknkkkbabfxfxfafhdxxf= Sn)(2180)4(4 fhabfR 注：注：为方便编程，可采用另一记法：令为方便

19、编程，可采用另一记法：令 n = 2n 为偶数，为偶数， 这时这时 ，有，有hkaxhnabhk ,2 )()(2)(4)(3 koddkevenkknbfxfxfafhS三、收敛速度与误差估计：三、收敛速度与误差估计：定义定义 若一个积分公式的误差满足若一个积分公式的误差满足 且且C 0，则则称该公式是称该公式是 p 阶收敛阶收敛的。的。 ChfRphlim0)(,)(,)(642hOChOShOTnnn例例4：计算计算dxx 10142 解：解： )1()(2)0(161718fxffTkk8kxk 其中其中= 3.138988494 )1()(2)(4)0(241oddeven4fxfx

20、ffSkk8kxk 其中其中= 3.141592502运算量基运算量基本相同本相同问题问题: 给定精度给定精度 e e，如何取，如何取 n ?例如：要求例如：要求 ，如何判断，如何判断 n = ？e e |nTI)()(122 fabhfR ? nkkhfh12)(12 )()(12)(1222afbfhdxxfhba 上述上述例例4中若要求中若要求 , 则则610| nTI622106| )0() 1 (|12| | hffhfRn00244949.0 h即：取即：取 n = 409通常采取将区间通常采取将区间不断对分不断对分的方法，即取的方法，即取 n = 2k上述上述例例4中中2k 40

21、9 k = 9 时，时，T512 = 3.14159202例例4中：中：S4 = 3.141592502注意到区间再次对分时注意到区间再次对分时412)()(12122fRhafbffRnn 412 nnTITI)(3122nnnTTTI 可用来判断迭代可用来判断迭代是否停止。是否停止。(1)(2)(3)事后误差估计事后误差估计一、梯形法的递推化一、梯形法的递推化逐次分半法逐次分半法 上一节介绍的复化求积方法对提高精度是行之有效的，但上一节介绍的复化求积方法对提高精度是行之有效的，但在使用求积公式之前必须给出合适的步长，在使用求积公式之前必须给出合适的步长，步长步长取得取得太大精度太大精度难以

22、保证难以保证，步长太小步长太小则会导致则会导致计算量计算量的的增加增加，而事先给出一个，而事先给出一个恰当的步长又往往是困难的恰当的步长又往往是困难的 实际计算中常常实际计算中常常采用变步长的计算方案采用变步长的计算方案，即在步长，即在步长逐次分逐次分半半(即步长二分即步长二分)的过程中，反复利用复化求积公式进行计算，的过程中，反复利用复化求积公式进行计算，直至所求得的积分值满足精度要求为止直至所求得的积分值满足精度要求为止 设将求积区间设将求积区间a，b分成分成n等分，则一共有等分，则一共有n+1个分点，按个分点，按梯形公式计算积分值梯形公式计算积分值Tn，需要提供，需要提供n+1个函数值如

23、果将求积个函数值如果将求积区间再二分一次，则分点增至区间再二分一次，则分点增至2n+1个，我们来个，我们来考察考察二分二分前后两前后两个积分值个积分值之间的之间的联系联系4 4 龙贝格求积公式龙贝格求积公式逐次分半逐次分半计算计算方案方案的实现的实现: 注意到每个子区间注意到每个子区间xk，xk+1经过二分只增加了一个分经过二分只增加了一个分点点 xk+1/2( xk+xk+1)/2，用复化梯形公式求得该子区间上的，用复化梯形公式求得该子区间上的积分值为积分值为 101021102110122)12(221)(221)(2)()(4nknnkknnkknkkknhkafhTxfhTxfhxfx

24、fhT)()(2)(4121 kkkxfxfxfh这里这里 代表二分前的步长代表二分前的步长. .将每个子区间上的积分值将每个子区间上的积分值相加得相加得nabh 二、龙贝格算法二、龙贝格算法).,()(212);,()(12)(222bafhabTIbafhabTIfRnnn 有有：根据复化梯形公式的余项表达式根据复化梯形公式的余项表达式. )(31.41)()(222nnnnnTTTITITIff 整理后可得：整理后可得：，则有，则有假定假定 可见，可见，利用利用两种步长两种步长计算的结果能估计截断误差计算的结果能估计截断误差.若将该截断若将该截断误差加到计算结果中误差加到计算结果中，nn

25、nnnTTTTTT3134)(31222 就得出就得出“改进的梯形求积公式改进的梯形求积公式”：事后误差事后误差估计估计例：例：计算计算dxx 10142 已知对于已知对于e e = 10 6 须将区间对分须将区间对分 9 次，得到次，得到 T512 = 3.14159202由由 来计算来计算 I 效果是否好些？效果是否好些？nnnnTTTTI313414422 483134TT = 3.141592502 = S4改进梯形求积公式改进梯形求积公式的右边实际是的右边实际是nnknkkknkknkknkknnnkknnnSbfxfxfafhxfhbfxfafhxfhTTxfhTTT 101121

26、102111102110212)()(2)(4)(6)(2)()(2)(231)(231)(221431)4(31这就是说用这就是说用梯形法二分前后的两个积分值梯形法二分前后的两个积分值Tn与与T2n的的线性组合线性组合的结果的结果得到得到复化辛普森法求积公式复化辛普森法求积公式nnnnnTTTTS141144313422 类似的情况，用辛普森法二分前后的两个积分值类似的情况，用辛普森法二分前后的两个积分值Sn与与S2n的线性组合的结果可得到的线性组合的结果可得到复化柯特斯求积公式复化柯特斯求积公式nnnnnSSSSC151151614114422222 重复同样的手续，用柯特斯法二分前后的两

27、个积分值重复同样的手续，用柯特斯法二分前后的两个积分值Cn与与C2n的线性组合的结果可得到的线性组合的结果可得到龙贝格龙贝格(Romberg)求积公式求积公式nnnnnCCCCR631636414114423233 我们在变步长的过程中运用加速公式，就能将粗糙的梯我们在变步长的过程中运用加速公式，就能将粗糙的梯形值形值Tn逐步加工成精度较高的辛普森值逐步加工成精度较高的辛普森值Sn 、柯特斯值、柯特斯值Cn和龙和龙贝格值贝格值Rn .一般有：一般有：nnnSTT 1442nnnCSS 144222nnnRCC 144323 Romberg 算法：算法： e e ? e e ? e e ? T1

28、 =)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(1TRomberg 序列序列kk2kT212 kS22 kC32 kR0 20=1 T11 21=2 T2 S12 22=4 T4 S2 C13 23=8 T8 S4 C2 R14 24=16 T16 S8 C4 R25 25=32 T32 S16 C8 R4 区间等分数区间等分数 梯形序列梯形序列 辛普森序列辛普森序列 柯特斯序列柯特斯序列 龙贝格序列龙贝格序列 龙贝格求积算法可用下表来表示：龙贝格求积算

29、法可用下表来表示： 例例5 用龙贝格方法计算椭圆用龙贝格方法计算椭圆 x2/4 + y2 l 的周长，使结果的周长，使结果具有五位有效数字具有五位有效数字 分析分析 为便于计算，先将椭圆方程采用参数形式表示为便于计算，先将椭圆方程采用参数形式表示, ,再根再根据弧长公式将椭圆周长用积分形式表示由于计算结果要求具据弧长公式将椭圆周长用积分形式表示由于计算结果要求具有五位有效数字，因此需要估计所求积分值有几位整数，从而有五位有效数字，因此需要估计所求积分值有几位整数，从而确定所求积分值的绝对误差限最后再应用龙贝格方法计算积确定所求积分值的绝对误差限最后再应用龙贝格方法计算积分分 解解 令令 x 2

30、cosq q，y sinq q ， 则椭圆的周长为则椭圆的周长为Iyxl4d sin314d42022022 q qq qq qq qq q.10125. 01081)(1021)(4422d sin3124451202 fRIfRIlI的的截截断断误误差差为为故故计计算算，的的截截断断误误差差为为则则需需结结果果有有五五位位有有效效数数字字，有有一一位位整整数数，要要求求，因因此此由由于于 q qq q 下表给出了用龙贝格方法计算积分下表给出了用龙贝格方法计算积分I= 1+1+3sin2q q dx 的过程的过程. /20kk2kT212 kS22 kC32 kR4322 kkRR0 1 2

31、.356 1941 2 2.419 921 2.441 1632 4 2.422 103 2.422 830 2.421 608 3 8 2.422 112 2.422 115 2.422 067 2.422 074 4 16 2.422 112 2.422 112 2.422 112 2.422 113 0.000 0395 32 2.422 112 2.422 112 2.422 112 2.422 112 0.000 001 0.125 10- -4 故积分故积分I 2.422112, 椭圆周长的近似值为椭圆周长的近似值为l = 4I 9.6884。三、理查森三、理查森(Richards

32、on)外推加速法外推加速法 上面讨论说明由梯形公式出发上面讨论说明由梯形公式出发, 将区间将区间a, b逐次二分逐次二分可提高求积公式的精度可提高求积公式的精度, 上述加速过程还可继续下去上述加速过程还可继续下去. 下面我们讨论其下面我们讨论其理论依据理论依据. ,)(24221 llhhhIhT .)(122nabhbafhabTIn ， 22hTTn若记若记Tn = T(h), 当区间当区间a, b分为分为2n等分时等分时, 有有 , 则则可见可见I = T(h)的误差为的误差为O(h2). llhhhIhT2422121642 3)(24)(1hThThT 若记若记 ，则，则 将梯形公式

33、按余项展开将梯形公式按余项展开. 由误差公式有由误差公式有 62411)(hhIhT 6416262411hhIhT 显然显然T1(h)与与 I 近似的阶为近似的阶为O(h4) . 就是就是辛普森公式辛普森公式序列序列Sn, S2n, . .， 2),(11hThT这样构造的这样构造的 )(1412144)(11hThThTmmmmmm 则又可进一步从余项中则又可进一步从余项中消去消去 h4 项，这样构造出的项，这样构造出的 ，其实就是，其实就是柯特斯公式柯特斯公式序序列，它与列，它与 I 的逼近阶为的逼近阶为O(h6) . )(2hT)(151 21516)(112hThThT 若令若令 ,

34、 一般地，若记一般地，若记T0(h) = T(h)，经过，经过m (m =1,2,)次加速次加速后，则有后，则有如此继续下去，每加速一次，误差的量级便提高如此继续下去，每加速一次，误差的量级便提高2阶阶. )21(141144)(1)1(1)()(0)()(0，次次加加速速值值，可可得得的的序序列列表表示示以以次次后后求求得得的的梯梯形形值值，且且表表示示二二分分设设以以 kTTTmTTkTkmmkmmmkmkkmk. ., 321.数数表表来来计计算算构构造造一一个个三三角角形形数数表表根根据据公公式式可可以以逐逐行行龙龙贝贝格格序序列列公公式式辛辛普普森森、柯柯特特斯斯、即即可可得得到到加

35、加速速、若若取取算算法法上上式式也也称称为为龙龙贝贝格格求求积积Tm Romberg 算法算法 可以证明，如果可以证明，如果 f (x) 充分光滑，那么充分光滑，那么T 数表每一列的数表每一列的元素及对角线元素均收敛到所求的积分值元素及对角线元素均收敛到所求的积分值 I ，即，即ITmITkmmkmk )()(lim)(lim，固定固定 理查德森理查德森外推法外推法利用利用低低阶公式产生阶公式产生高高精度的结果。精度的结果。设对于某一设对于某一 h 0，有公式，有公式 T0(h) 近似计算某一未知值近似计算某一未知值 I。由。由Taylor展开得到：展开得到： T0(h) I = 1 h +

36、2 h2 + 3 h3 + i 与与 h 无关无关现将现将 h 对分，得：对分，得：( () )( () )( () ).)(3232222120 hhhhIT Q：如何将公式精度由如何将公式精度由 O(h) 提高到提高到 O(h2) ？.432112)()(23322020 hhIhTTh 即：即：.12)()(2)(32210201 hhIhTThTh .)(42312 hhIhT 12)()(221212 hTTh.)(2211 mmmhhIhT 12)()(2121 mmhmmhTTHW: p.159 6,7,8.3) 在构造在构造Newton-Cotes公式公式时，限定用积分区间的时

37、，限定用积分区间的等分点等分点作为求积节点作为求积节点，这样做虽然使问题的处理过程得以简化，但，这样做虽然使问题的处理过程得以简化，但同时也同时也限制了精度限制了精度。 求积公式含有求积公式含有2n+2个待定参数个待定参数xk、Ak(k0,1,n)若用若用待定系数法确定它们待定系数法确定它们, 则最好需要则最好需要2n+2个独立的条件联立方个独立的条件联立方程组求解程组求解, 从而易知求积公式的从而易知求积公式的最大代数精度最大代数精度可达到可达到2n+1次次. 在节点数目固定为在节点数目固定为n 的条件下，能否通过的条件下，能否通过适当选取求积适当选取求积节点节点xk的位置以及相应的求积系数

38、的位置以及相应的求积系数Ak，使求积公式，使求积公式具有尽可能高具有尽可能高(最高最高)的代数精度？的代数精度？ bankkkxfAxxf0)(d)(这类求积公式称为这类求积公式称为高斯高斯(Gauss)求积公式求积公式。4 高斯求积公式高斯求积公式 将节点将节点 x0 xn 以及系数以及系数 A0 An 都作为待定系数。都作为待定系数。令令 f (x) = 1, x, x2, , x2n+1 代入可求解，得到的公式代入可求解，得到的公式具有具有2n+1 次代数精度。这样的节点称为次代数精度。这样的节点称为Gauss 点点，公式称为公式称为Gauss 型求积公式型求积公式。 baxxxfId)

39、()( 为使问题更具一般性为使问题更具一般性,我们研究带权积分我们研究带权积分 (x)为权函数为权函数, Ak(k0,1,n)为不依赖于为不依赖于f (x)的求积系数的求积系数, xk (k0,1,n)为求积节点为求积节点. bamnkmkknmxxxxA)2.5(.12,1,0d)(0 要使要使(5.1)具有具有2n+1次代数精度，则需要满足次代数精度，则需要满足 bankkkxfAdxxfx0)()()( 构造具有构造具有2n+1次代数精度的求积公式次代数精度的求积公式(5.1) 从例中可看到求解非线性方程组从例中可看到求解非线性方程组(5.2)较复杂，通常较复杂，通常n2就很难求解故一般

40、不通过解方程就很难求解故一般不通过解方程(5.2)求求 xk 及及 Ak (k0,1, , n)例：例：求求 的的 2 点点 Gauss 公式。公式。dxxfx)(10 解：解：设设 ，应有，应有 3 次代数精度。次代数精度。 101100)()()(xfAxfAdxxfx代入代入 f (x) = 1, x, x2, x3 31130092211200721100521032xAxAxAxAxAxAAA2776. 03891. 02899. 08212. 01010 AAxx不是线性方程组，不是线性方程组，不易求解。不易求解。 而从研究而从研究高斯点的基本特性高斯点的基本特性来着手解决来着手解

41、决Gauss 求积公式求积公式的构造问题的构造问题由插值型公式构由插值型公式构造知造知,关键求关键求xk，一、高斯点的基本特性一、高斯点的基本特性0)(d)()()(1 banxxxxxP 证明：证明： “” x0 xn 为为 Gauss 点点, 则公式则公式 至少有至少有 2n+1 次代数精度。次代数精度。 bankkkxfAdxxfx0)()()( 对任意次数对任意次数不大于不大于n 的多项式的多项式 Pm(x)， Pm(x) w(x)的次数的次数不大于不大于2n+1，则代入公式应则代入公式应精确成立精确成立： nkkkmkbamxwxPAdxxwxPx0)()()()()( 0= 0 “

42、” 要证明要证明 x0 xn 为为 Gauss 点，即要证公式对任意次点，即要证公式对任意次数数不大于不大于2n+1 的多项式的多项式 Pm(x) 精确成立，即证明：精确成立，即证明： nkkmkbamxPAdxxPx0)()()( 设设)()()()(xrxqxwxPm bababamdxxrxdxxqxwxdxxPx)()()()()()()( 0 nkkkxrA0)( nkkmkxPA0)( x0 xn 为为 Gauss 点点 与任意次数与任意次数不大于不大于n 的多项式的多项式 P(x) （带权）正交（带权）正交。 nkkxxxw0)()(定理定理求求 Gauss 点点 求求w(x)的

43、的零点零点 Gauss 公式的余项：公式的余项： bankkkxfAdxxffR0)()(/* 设设P为为f 的过的过x0 xn的插值多项式的插值多项式 */ bankkkxPAdxxf0)()(/*只要只要P 的阶数不大于的阶数不大于2n+1，则下一步，则下一步等式成立等式成立*/dxxPxfdxxPdxxfbababa)()()()( 插值多项式插值多项式的余项的余项Q：什么样的什么样的插值多项式插值多项式在在 x0 xn 上有上有 2n+1 阶？阶？A：Hermite 多项式！多项式！ 满足满足)()(),()(kkkkxfxHxfxH badxxHxffR)()(),(,)()!22(

44、)()()!22()(2)12(2)12(badxxwnfdxxwnfbanbaxn 二、高斯求积公式的余项二、高斯求积公式的余项三、高斯求积公式的稳定性与收敛性三、高斯求积公式的稳定性与收敛性 定理定理6 高斯求积公式高斯求积公式(5.1)的求积系数的求积系数 Ak (k0,1,n)全是正的全是正的 由本定理及定理由本定理及定理2，则得，则得 推论推论 高斯求积公式高斯求积公式(5.1)是稳定的是稳定的. 定理定理7 设设 f (x)C a，b，则高斯求积公式，则高斯求积公式(5.1)是收敛是收敛 的，即的，即 nkbakknxxxfxfA0.d)()()(lim 正交多项式族正交多项式族

45、0, 1, , n, 有性质：任意次数不大有性质：任意次数不大于于n 的多项式的多项式 P(x) 必与必与 n+1 正交。正交。若取若取 w(x) 为其中的为其中的 n+1，则，则 n+1的根的根就是就是 Gauss 点。点。再解上例：再解上例： 101100)()()(xfAxfAdxxfxStep 1：构造正交多项式构造正交多项式 2设设cbxxxaxxx 2210)(,)(, 1)( 53 a0)(10 dxaxx0),(10 1021102100)(53(0),(0)(0),(dxcbxxxxdxcbxxx 215910 cb即：即：215910)(22 xxx 四、常用的高斯型求积公

46、式四、常用的高斯型求积公式Step 2：求求 2 = 0 的的 2 个根，即为个根，即为 Gauss 点点 x0 ，x1221/20)9/10(9/1021;0 xStep 3：代入代入 f (x) = 1, x 以求解以求解 A0 ，A1解解线性线性方程组，方程组，简单。简单。结果与前一方法相同：结果与前一方法相同：2776. 0,3891. 0,2899. 0,8212. 01010 AAxx 利用此公式计算利用此公式计算 的值的值 10dxexx2555. 1 10dxexx2899. 08212. 0102776. 03891. 010eeeAeAxx 注：注：构造正交多项式也可以利用

47、最小二乘数据拟合中介绍构造正交多项式也可以利用最小二乘数据拟合中介绍过的递推式进行。过的递推式进行。 特殊正交多项式族：特殊正交多项式族： Legendre 多项式族：多项式族：1)( x 定义在定义在 1, 1上，上，kkkkkxdxdkxP)1(!21)(2 满足：满足： lklkPPklk1220),(xPP 10, 1由由 有递推有递推11)12()1( kkkkPxPkPk以以 Pn+1 的根为节点的求积公式称为的根为节点的求积公式称为Gauss-Legendre 公式公式。 Chebyshev 多项式族：多项式族：211)(xx 定义在定义在 1, 1上，上，) arccos( c

48、os)(xkxTk Tn+1 的根的根为为 2212cosnkxkk = 0, , n以此为节点构造公式以此为节点构造公式 1102)()(11nkkkxfAdxxfx称为称为 Gauss-Chebyshev 公式公式。注意到积分端点注意到积分端点 1 可能是积分可能是积分的的奇点奇点，用普通，用普通Newton-Cotes公公式在端点会出问题。而式在端点会出问题。而Gauss公公式可能避免此问题的发生。式可能避免此问题的发生。其它公式见教材其它公式见教材p.144-148注：注：一般一般a,b上的积分可化为上的积分可化为-1,1上特殊高斯公式进行计算。上特殊高斯公式进行计算。5 5 数值微分

49、数值微分 数值微分的数值微分的概念概念 数值微分的数值微分的计算方法计算方法 原始概念近似原始概念近似: :中点法及外推法中点法及外推法 函数近似函数近似: :插值型的求导公式插值型的求导公式 函数相互关系转化函数相互关系转化: :利用数值积分求导利用数值积分求导 数值微分的数值微分的误差分析误差分析 泰勒展式估计泰勒展式估计 事后误差估计事后误差估计 基本关系转化基本关系转化 数值微分数值微分就是就是用函数值的线性组合近似函数在某点用函数值的线性组合近似函数在某点的导数值的导数值一、中点法和外推法一、中点法和外推法 按导数定义按导数定义 , 是差商是差商 当当 时的极限时的极限取取差商差商作

50、为作为导数导数的近似值的近似值，建立简单的数值微分方法，建立简单的数值微分方法:)(0 xf hxfhxf)()(00 0h( () )( () )( () )hxfhxfxf000 (6.1)向后差商近似导数向后差商近似导数(6.2)(6.3)中心差商近似导数中心差商近似导数( () )( () )( () )hhxfxfxf 000( () )( () )( () )hhxfhxfxf2000 ( ( ) )( () )( () )hhxfhxfhD200 容易看出，就精度而言，以（容易看出，就精度而言，以（6.3）式更为可取，称）式更为可取，称(6.4)为为 的的中点公式中点公式, 其中

51、其中h为一增量，称为为一增量，称为步长步长 这种数值这种数值微分方法称为微分方法称为中点方法中点方法, 它是前两种方法的算术平均它是前两种方法的算术平均)(0 xf 分别将分别将在在 x=a 处做处做Taylor展开有展开有)(haf )(! 5)(! 4)(! 3)(! 2)()()() 5(5) 4(432afhafhafhafhafhafhaf代入代入D(h)得得 )(! 5)(! 3)()()5(42afhafhafhD,6)()(2MhhDaf 其中其中)(maxxfMhax 现在来考虑中点公式现在来考虑中点公式 的截断误差，的截断误差，hhafhafhD2)()()( (6.5)所

52、以所以截断误差截断误差 )(! 5)(! 3)()() 5(42afhafhafhD(6.6)从截断误差的角度来看，步长从截断误差的角度来看，步长h越小，计算结果越准确。且越小，计算结果越准确。且 但从计算角度看，但从计算角度看，h 越小，越小， f (a+h)与与 f (a- -h) 越接近，直越接近，直接相减会造成有效数字的严重损失接相减会造成有效数字的严重损失(参看第参看第1章第章第4节节)。因此。因此, 从舍入误差的角度来看，步长从舍入误差的角度来看，步长 h 不宜太小。不宜太小。 所以所以, 在在实际计算时实际计算时，通常，通常采用采用二分步长二分步长及及误差事后估误差事后估计法计法

53、, 在变步长的过程中实现步长的自动选择，在保证截断在变步长的过程中实现步长的自动选择，在保证截断误差满足的精度要求的前提下选取取尽可能大的步长。误差满足的精度要求的前提下选取取尽可能大的步长。kD(h)01239103.017652.791352.736442.722812.718282.71828 解解 这里采用的计算公式是这里采用的计算公式是 计算结果见表计算结果见表6.1，表中，表中 k 代表二分的次代表二分的次数，步长数，步长 。二分。二分 9 次得结果次得结果 D= 2.71828，它的每一数字都是有效数字，它的每一数字都是有效数字(所所求导数的准确值为求导数的准确值为e=2.718

54、2818)。( ( ) )heehDhh211 kh28 . 0M 例例6.1 用用变步长变步长的中点方法求的中点方法求 在在x=1处的导数值处的导数值,设取设取h= 0.8起算。起算。xe表表6.1 计算结果计算结果 210)()(hxfhD 4)()2()()(00 xfhDxfhD21041)()2(hxfhD )()2(34)()(0hDhDhDxf 事后误差估计事后误差估计 我们看到，中点公式具有如下形式我们看到，中点公式具有如下形式 (6.7)式中的系数均与步长无关。式中的系数均与步长无关。 6342210)()(hhhxfhD 若将步长二分，则有若将步长二分，则有 (6.8) 6

55、34221064116141)()2(hhhxfhD 41h 若令若令 (6.10)则进一步消去则进一步消去D1(h)误差主项误差主项 , 有有 )(151)2(1516)(112hDhDhD 6102)()(hxfhD2141h则可消去误差主项则可消去误差主项 ，得，得 )(31)2(34)(1hDhDhD 624101)()(hhxfhD取取(6.7)与与(6.8)加加权平均权平均(6.9)重复同样的手续重复同样的手续，再导出下列加速公式，再导出下列加速公式 这种加速过程还可继续下去。这种加速过程还可继续下去。这种加速方法通常称作这种加速方法通常称作Richardson外推加速法外推加速法

56、。 )(631)2(6364)(223hDhDhD 例例6.2 运用加速公式和加工例运用加速公式和加工例6.1的结果。的结果。表表6.2 Richardson外推加速法计算结果外推加速法计算结果)(1hD)(2hD)(3hDh D(h)0.80.40.20.13.017652.791352.736442.722812.7159172.7181372.7182672.7182852.7182762.71828 解解 计算结果计算结果见表见表6.2。这里，。这里，加速的效果同样加速的效果同样是相当显著的。是相当显著的。hhaGafafe ee ee e 2)()()(21最小最小,步长步长h不宜太

57、大不宜太大,也不宜太小也不宜太小. 其其最优步长最优步长应为应为 它表明它表明h 越小越小, 舍入误差舍入误差)(af 越大越大, 故它是病态的故它是病态的. 用中点用中点公式计算公式计算)(af 的的误差上界误差上界为为,6)(2hMhhEe e 要使误差要使误差E(h)3opt/3Mhe e 则计算则计算 当当 f (a+h) 及及 f (a-h) 分别有舍入误差分别有舍入误差 e e1 及及 e e2 时，若令时，若令 e e)(af 的舍入误差上界为的舍入误差上界为 21,maxe ee e 注注: 中点心公式及其加速方法适合用表达式表示的函数。对中点心公式及其加速方法适合用表达式表示

58、的函数。对于列表函数，则宜使用插值方法等导出数值求导公式。于列表函数，则宜使用插值方法等导出数值求导公式。 二、插值型的求导公式二、插值型的求导公式 x x0 x1 x2 xn y y0 y1 y2 yn 对于列表函数对于列表函数 y = f (x)： 插值多项式插值多项式 y = Pn(x)作为它的近似作为它的近似, 我们取我们取统称统称插值型的求导公式插值型的求导公式)()(xPxfn )(xPn )(xf 的近似值，的近似值，作为作为建立的数值公式建立的数值公式 依据插值余项定理，依据插值余项定理，求导公式求导公式(6.11)的的余项余项为为)(dd)!1()()()!1()()()()

59、1(11)1( nnnnnfxnxxnfxPxf式中式中. )()(01 niinxxx (6.11) 我们限定我们限定:求某个节点求某个节点 xk 上的导数值，上的导数值，上面的第二项变为零上面的第二项变为零，这时有余项公式，这时有余项公式)()!1()()()(1)1(knnknkxnfxPxf (6.12) 1 1两点公式两点公式 已给两节点已给两节点 x0, x1 上的函数值上的函数值 f (x0), f (x1)，做线性插值，做线性插值)()()(101001011xfxxxxxfxxxxxP 记记 x1 x0 = h，对上式两端求导，有，对上式两端求导，有)()(1)(101xfx

60、fhxP )()(1)( )()(1)(01110101xfxfhxPxfxfhxP 于是有下列求导公式：于是有下列求导公式：)(2)()(1)()(2)()(1)(011010 fhxfxfhxffhxfxfhxf 而利用余项公式而利用余项公式(6. 4)知，带知，带余项的两点公式余项的两点公式是：是： 下面我们仅仅考察下面我们仅仅考察节点处的导数值节点处的导数值为简化讨论为简化讨论, 假定假定所给的节点是等距的所给的节点是等距的 设已给出三节点设已给出三节点x0, xl=x0+h, x2=x0+2h上的函数值上的函数值,做二次插值做二次插值)()()()()()()()()()(21202