数值分析讲义第四章数值积分

《数值分析讲义第四章数值积分》由会员分享，可在线阅读，更多相关《数值分析讲义第四章数值积分（77页珍藏版）》请在装配图网上搜索。

1、第四章 数值积分与数值微分 /* Numerical Integration and differentiation*/近似计算 ba dxxfI )( 1 引言 对f()采用不同的近似计算方法，从而得到各种不同的求积公式。 以上三种方法都是用被积函数值的线性组合来表示积分值。推广，一般地有0( ) ( )nb k ka kf x dx A f x求积节点求积系数，与被积函数无关像这样，将积分用若干节点上被积函数值的线性组合来表示的数值积分公式称为机械求积公式。 0 ( ) ( )nb k ka kR f f x dx A f x 求积误差 机械型求积公式的构造归结为，确定求积节点xk和求积系

2、数Ak，使在某种意义下精确度较高。总之，要解决三个问题：精确度的度量标准；如何构造具体的求积公式；1.具体求积公式构造出来后，误差如何估计？ 以上三种方法都是用被积函数值的线性组合来表示积分值。推广，一般地有0( ) ( )nb k ka kf x dx A f x求积节点求积系数，与被积函数无关像这样，将积分用若干节点上被积函数值的线性组合来表示的数值积分公式称为机械求积公式。 0 ( ) ( )nb k ka kR f f x dx A f x 求积误差 机械型求积公式的构造归结为，确定求积节点xk和求积系数Ak，使在某种意义下精确度较高。总之，要解决三个问题：精确度的度量标准；如何构造具

3、体的求积公式；1.具体求积公式构造出来后，误差如何估计？ 以上三种方法都是用被积函数值的线性组合来表示积分值。推广，一般地有0( ) ( )nb k ka kf x dx A f x求积节点求积系数，与被积函数无关像这样，将积分用若干节点上被积函数值的线性组合来表示的数值积分公式称为机械求积公式。 0 ( ) ( )nb k ka kR f f x dx A f x 求积误差 机械型求积公式的构造归结为，确定求积节点xk和求积系数Ak，使在某种意义下精确度较高。总之，要解决三个问题：精确度的度量标准；如何构造具体的求积公式；1.具体求积公式构造出来后，误差如何估计？ 定义：代数精度若某个求积公

4、式对次数 m 阶的多项式准确成立，而对m+1 阶的多项式不一定准确成立。即对应的误差满足：R Pk =0 对任意 k m 阶的多项式成立，且 R Pm+1 0 对某个 m+1 阶多项式成立，则称此求积公式的代数精度为 m 。代数精度与误差的关系：代数精度越高，求积误差越小。结论：问题1 由上面代数精度条件确定求积公式可分两种情形：若事先给定求积节点xk(k=0,n),例如被积函数以表的形式给出时xk确定，可令m=n,由上式确定n+1个系数Ak即可-待定系数法和插值法。1.若xk和Ak都可选择，令m=2n +1，确定xk和法Ak -Gauss法要使求积公式具有m阶代数精度，则它对1,x,xm均准

5、确成立，即 0 2 20 1 10 121 1n kkn k kkn m m mk kk A b aA x b aA x b am m+1个方程，2n+2个未知数问题2 Case 1-方法1 1 插值型求积 公式思路利用插值多项式 则积分易算。)()( xfxPn 在a, b上取 a x0 x1 xn b，做 f 的 n 次插值多项式 ，即得到 nk kkn xlxfxL 0 )()()( ba ba knk k dxxlxfdxxf )()()( 0 Ak ba kj xx xxk dxA jk j )( )(由 决定，与 无关。节点 f (x)插值型积分公式/*interpolatory

6、quadrature*/ ba nk kxn ba nba nba nk kk dxxxnf dxxRdxxLxf xfAdxxf fR 0)1( 0 )()!1( )( )()()( )()( 误差 Case 1-方法2 1 Newton-Cotes Formulae例：对于a, b上1次插值，有)()()(1 bfafxL ab axba bx )()()( 2221 bfafdxxfAA abbaab 考察其代数精度。f(x)a bf(a) f(b)梯形公式/* trapezoidal rule*/解：逐次检查公式是否精确成立代入 P0 = 1： ba abdx1 112 ab=代入 P

7、1 = x ：=代入 P2 = x2 ：2 22 abba dxx 2 baab 32 33 abba dxx 222 baab 代数精度 = 1 1 Newton-Cotes FormulaeTh1.形如 的求积公式至少有 n 次代数精度 该公式为插值型（即： ）nk kk xfA0 )( ba kk dxxlA )( 当节点等距分布时：ninabhhiaxi ,.,1,0, dxxx xxA nxx ij ji ji 0 )( )( n jiinn ji dtjtinin abdthhji hjt 00 )()!(! )1)()( )(令htax Cotes系数)(n iC注：Cotes

8、系数仅取决于 n 和 i，可查表得到。与 f (x) 及区间a, b均无关。 2 Newton-Cotes 公式 ( ) 00( ) ( )nb nka kf x dx b a C f x kh NewtonCotes formula 1 Newton-Cotes Formulae21,21 )1(1)1(0 CCn = 1: )()(2)( bfafabdxxfba Trapezoidal RuledxbxaxffR ba x )(!2 )( /* 令 x = a+th, h = ba, 用中值定理 */1,)(121 3 abhbafh 代数精度 = 1n = 2: 61,32,61 )2

9、(2)2(1)2(0 CCC )()(4)(6)( 2 bffafabdxxf baba Simpsons Rule代数精度 = 32,),(,)(901 )4(5 abhbafhfR n = 4: Cotes Rule, 代数精度 = 5, )(9458 )6(7 fhfR 0 1 2 3 4( ) 7 ( ) 32 ( ) 12 ( ) 32 ( ) 7 ( )90ba b af x dx f x f x f x f x f x 偶数阶N-C公式具有n+1阶代数精度n = 3: Simpsons 3/8-Rule, 代数精度 = 3, )(803 )5(5 fhfR 对称节点的系数相同Co

10、tes公式是用不同节点的函数值（高度）的加权平均来近似区间的平均高度注：当n8时，Cotes系数有负，造成公式不稳定，因此常用低阶Cotes公式。 ( 1) 0 020 0 0 0 /22 /2 0 ( ) ( ) ( )( 1)!, ( ) even, /2 integer, let /2, we have ( /2 )n n nb bx k ka ak k nnnk k nnn n k fR f x x dx x x dxnx x jh x x th R f h t k dtn n t u nR f h u n k du 证明：只需证明n为偶数时， N-C公式对f(x)=xn+1的余项R(

11、f)=0即可。因 f(n+1)(x)=(n+1)!, 由余项公式得Th2. n为偶数时， N-C公式至少具有n+1阶代数精度。 0 0/20 /2/2 /21/2 /2 /2 /2/2 /2 ( ) ( /2 ) ( /2)( /2) ( /2)( ) 1 ( )( ) odd 0.n nk kn nk j nn nnj n j n n nj n j n let H u u n k u k nu k n u j let j k nH u u j u ju j u j H uH u R f 注：当n 为偶数时，Cotes公式具有n+1阶精度，与n+1阶Cotes公式精度相同，但少计算一个节点上的

12、函数值， 因此一般常用偶数阶Cotes公式。 偶数阶N-C公式具有n+1阶代数精度N-C公式具有n阶代数精度余项R=o(h n+2) Hint：construct a interpolation polynomial of order 5, H(x), satisfying H(a)=f(a), H(b)=f(b), H (k)(a+b)/2) = f (k)(a+b)/2). HW: p.151-152 #1-6 数值稳定性的一般概念 N-C的稳定性 3 复合求积 /* Composite Quadrature */Havent we had enough formulae? Whats u

13、p now?Oh come on, you dont seriously consider h=(ba)/2 acceptable, do you?Why cant you simply refi e the partition ifyou have to be s picky?Dont ou forget the oscillatory nature of high-degree polynomias!Uh-oh高次插值有Runge 现象，故采用分段低次插值 分段低次合成的 Newton-Cotes 复合求积公式。 复合梯形公式：),.,0(, nkhkaxnabh k 在每个 上用梯形公式

14、：, 1 kk xx nkxfxfxxdxxf kkxx kkkk ,.,1,)()(2)( 111 11 )()(2)(2 nk k bfxfafh ba nk kk xfxfhdxxf 1 1 )()(2)( = Tn),(),()(12 )()(12)(12 2 121 3 bafabh nfabhfhfR nk knk k /*中值定理*/ 2 Composite Quadrature 复化 Simpson 公式：),.,0(, nkhkaxnabh k )()(4)(6)( 1211 kkkxx xfxfxfhdxxfkk kx 21kx 1kx 4 4 4 4 4)()(2)(4)

15、(6)( 10 10 121 nk nk kkba bfxfxfafhdxxf = Sn)(2180 )4(4 fhabfR 注：为方便编程，可采用另一记法：令 n = 2n 为偶数， 这时 ，有hkaxhnabh k ,2 )()(2)(4)(3 kodd keven kkn bfxfxfafhS 2 Composite Quadrature 收敛速度与误差估计：定义 若一个复化积分公式的误差满足 且C 0，则称该公式是 p 阶收敛的。 ChfR ph lim0 3 21 0 12 : ( ) ( )12 12 1 1( ) ( ) ( )12 12n nk kk kba h hR f f

16、hfR f f x dx f b f ah 复 化 梯 形 公 式 /*中值定理*/类似的，可得 4 4 (5) (5)6 6 1: ( ) ( )180 2 2: ( ) ( )945 4R f f b f ahR f f b f ah 复 化 Simpson公 式复 化 Cotes公 式 2阶收敛4阶收敛6阶收敛 例1：计算dxx 10 1 4 2解： )1()(2)0(161 718 fxffT k k 8kxk 其中= 3.138988494 )1()(2)(4)0(241 odd even4 fxfxffS kk 8kxk 其中= 3.141592502运算量基本相同 Q: 给定精度

17、 ，如何取 n ?例如：要求 ，如何判断 n = ？ | nTI 2 ( ) ( )12hR f f b f a 上例中若要求 ，则610| nTI 622 106|)0()1(|12| hffhfRn00244949.0h即：取 n = 409 2 Composite Quadrature 事后误差估计式，可用来判断迭代是否停止。始步长h 0.510-2 4 龙贝格积分 /* Romberg Integration */复化梯形公式算法简单，但精度较差，收敛速度（2阶收敛）较慢，如何提高收敛速度？ 注：按上面规律，可以构造线性组合系数为的新的积分公式，但当m4时，前一个系数接近于1，后一个系

18、数接近于0，这样构造出的新公式与前一个公式结果差别不大，反而增加计算量，因此实际上常做到Romberg公式为止。4 1,4 1 4 1mm m 例：计算dxx 10 1 4 2已知对于 = 106 须将区间对分 9 次，得到 T512 = 3.14159202由 来计算 I 效果是否好些？nnnn TTTTI 3134144 22 考察412 nnTI TI 48 3134 TT = 3.141592502 = S4一般有： nnn STT 144 2 nnn CSS 144 222 nnn RCC 144 323Romberg 序列 Romberg 算法： ? ? 0, k=0,1,n定理(

19、 )( )2 2 202 :( )( ) 2( ) ( ), ( )( ) 0jk jx xk x xj kk nb k i k i k i kia ib k kaproofl x Gauss Lagrangel x nl x dx Al x l xl x dx A 设 是 以 点 为 节 点 的 插 值 基 函 数 ，则 是 一 个 次 多 项 式 ， 因 此 求 积 公 式 对 它 准 确 成 立 ， 即又求积系数的另一种计算方法结论：Gauss型积分公式是数值稳定的。（稳定性分析类似于n 7的N-C公式） 区 间 -1,1上 权 函 数 (x)=1的 Gauss型 求 积 公 式 ,称

20、为Gauss-Legendre求 积 公 式 ,其 Gauss点 为 Legendre多 项 式 的零 点 . 几 种 Gauss型 求 积 公 式 (1) Gauss-Legendre求 积 公 式 公 式 的 Gauss点 和 求 积 系 数 可 在 数 学 用 表 中 查 到 .n xk Ak n xk Ak1 0 2 6 0.9324695142 0.6612093865 0.2386191861 0.17132449240.36076157300.46791393462 0.5773502692 13 0.77459666920 0.55555555560.8888888889 7

21、0.9491079123 0.7415311856 0.40584515140 0.12948496620.27970539150.38183005050.41795918374 0.8611363116 0.3399810436 0.34785484510.6521451549 8 0.9602898565 0.7966664774 0.5255324099 0.1834346425 0.10122853630.22238103450.31370664590.36268378345 0.9061798459 0.53846931010 0.23692688510.47862867050.56

22、88888889 用 3点 Gauss公 式 计 算 积 分 解 查 表 得 x1=-0.7745966692,x2=0,x3=0.7745966692, A1=A3=0.5555555556,A2=0.8888888889, 所 以 有 Gauss-Legendre求 积 公 式 的 余 项 为 )1,1(,)()12()!2( )!(2 )2(3 412 nn fnn nfR例 13 11 .cos xdxI 68300355.1coscoscos 332211 xAxAxAI误 差 为 53 47 103492.6)cos(7)!6( 62 fR实 际 上 ,I=2sin1=1.6829

23、4197, 误 差 为 |R|=6.15810-5 . 用 Simpson公 式 ,则 有 I1.69353487, 误 差 为 |R|=1.0610 -2 . 由 于因 此 ,a,b上 权 函 数 (x)=1的 Gauss型 求 积 公 式 为 用 3点 Gauss公 式 计 算 积 分结 果 远 比 Simpson公 式 的 结 果 精 确 . ba tabbaxdttabbafabdxxf )2 )()()22(2)( 11 ba ni ii xabbafAabdxxf 1 )22(2)(例 14 10 2 .1 4 dxxI 解 这 里 Gauss点 和 积 分 系 数 与 上 例 相

24、 同 ,所 以 31 210 2 141068.32/)1(1 4211 4 i ii xAdxxI求 积 误 差 可 表 示 为 ),(,)()12()!2( )!()( )2(3 412 bafnn nabfR nn 区 间 0,)上 权 函 数 (x)=e-x的 Gauss型 求 积 公 式 ,称为 Gauss-Laguerre求 积 公 式 ,其 Gauss点 为 Laguerre多 项 式的 零 点 . (2) Gauss-Laguerre求 积 公 式 公 式 的 Gauss点 和 求 积 系 数 可 在 数 学 用 表 中 查 到 .n xk Ak n xk Ak2 0.5858

25、8643763.4142135623 0.85355339050.1464466094 5 0.26356031971.41340305913.59642577107.085810005812.6408008442 0.52175561050.39866681100.07594244970.00361175870.00002337003 0.41577455672.2942803602602899450829 0.71109300990.27851773350.0103892565 6 0.22284660411.18893210162.9927363260 5.77514356919.837

26、467418315.9828739806 0.45896467930.41700083070.11337338200.01039919750.00026101720.00000089854 0.32254768961.74576110114.53662029699.3950709123 0.60315410430.35741869240.03888790850.0005392947 Gauss-Laguerre求 积 公 式 为 求 积 公 式 的 误 差 为 由 于 所 以 ,对 0, +)上 权 函 数 (x)=1的 积 分 ,也 可 以 构 造 类 似的 Gauss-Laguerre求

27、积 公 式 : 0 1 )()( ni iix xfAdxxfe ),0(,)()!2( )!( )2(2 nfnnfR 00 )()( dxxfeedxxf xx 0 1 )()( ni ixi xfeAdxxf i 区 间 (-,)上 权 函 数 (x)= (3) Gauss-Hermite求 积 公 式 公 式 的 Gauss点 和 求 积 系 数 可 在 数 学 用 表 中 查 到 .n xk Ak n xk Ak2 0.7071067811 0.8862269254 6 0.4360774119 1.3358490704 2.3506049736 0.72462959520.1570

28、6732030.00453000993 1.22474487130 0.29540897511.81635900064 0.5246476232 1.6506801238 0.80491409000.0813128354 7 0.8162878828 1.6735516287 2.65196135630 0.42560725260.05451558280.00097178120.81026461755 0.9585724646 2.02018287040 0.39361932310.01995324210.9453087204 2xe 的 Gauss型 求 积 公 式 ,称 为 Gauss-Hermite求 积 公 式 , 其 Gauss点 为 Hermite多 项 式的 零 点 . Gauss-Hermite求 积 公 式 为 ni iix xfAdxxfe 1 )()(2 ni ixi xfeAdxxf i1 )()( 2或 6 Numerical Differentiation h=0.8，精度像这样，将微分的计算归结为在若干节点上函数值的线性组 合的数值微分方法称为机械求导方法。 二次插值 f(xn)7.Gauss型积分公式8. HW: p.152-153 #7-12