Chap数值积分实用PPT课件

1、3.1 数值积分公式及其代数精度1. 数值积分公式的一般形式函数值为f (xk) (k = 0 , 1, . . . , n) , 取上述这些函数值的带权的和0( )()nnkkkIfA f x(1.1)作为 的近似值, ( )baf x dx0( )(),nbkkakf x dxA f x即令(1.2)设 f (x) 在节点 处的01naxxxb第1页/共76页称上式为数值积分公式。 xk (k = 0 , 1 , . . . , n) 称为求积节点。权 Ak 又称求积系数, Ak 仅与 xk 的选取有关。 0( )( )()nbkkakR ff xA f x称 为数值求积公式的余项（或截断

2、误差）。0( )(),nbkkakf x dxA f x(1.2)这类数值积分公式也称为机械求积公式，其特点是利用一些函数值的组合算出积分的近似值。 这就避开了牛顿-莱伯尼兹公式需要寻求原函数的困难。第2页/共76页 2.数值积分公式的代数精度 利用余项 R(f) 可以描述数值积分公式的精度, 而刻画其精度的另一概念是代数精度.定义1.1 若数值积分公式对于一切次数 m 的代数多项式, 都准确成立, 称其至少具有 m 次代数精度; 若数值积分公式对于一切次数m 的代数多项式都准确成立, 而对于某个 m +1 次代数多项式不准确成立, 则称此求积公式具有m 次代数精度.定理1.1数值积分公式具有

3、 m次代数精度的充分必要条件是当 f (x) =1, x , x2 ,. . . , xm 时数值积分公式准确成立, 而当 f (x) = xm+1 时其不准确成立.按以上定义易知第3页/共76页中的待定系数 A0 , x1 , A1 , 使其代数精度尽量高, 并指出所确定的求积公式的代数精度. 101102( )(0)()(1)3f x dxA ff xA f解：令 f (x) = 1 , x , x2 使之准确成立, 则有011121121,321,3221.33AAxAxA1011,21.6xAA101211( )(0)(1)6326( )f x dxfff例1.1. 确定下列数值积分公

4、式第4页/共76页取 f (x) = x3 时，上式左边 = 右边 = 1/4 。取 f (x) = x4 时，上式左边 = 1/5，右边 = 5/24，左边 右边。所以确定的求积公式具有3 次代数精度。101211( )(0)(1)6326( )f x dxfff下面确定数值积分公式的代数精度。若用类似的方法确定一般的数值积分公式如何？0( )d()nbkkakf xxA f x第5页/共76页22d1d2dd112222120121102200nabxxxAiabxxxAabxxxAabxAnnbannknkkiibainkikkbankkkbankk若有解，则得到的数值积分公式(1.2)

5、至少有2n+1次代数精度。0( )d()nbkkakf xxA f x但这是一个非线性方程组。很难求解！第6页/共76页3 插值型数值积分公式设 f (x) 在插值节点 a x0 x1 . . . 0，m为正整数，步长h =(b - a)/2m。即将积分区间分割成2m等份。Step 2. 计算 12012)()(2mmkkkxfxfhT这里 .2 , 1 , 0,mkkhkaxStep 3. 计算 1202122)(2211mmmkkxfhTT这里 . 12, 1 , 0,)21(21mkkhkax将每一个小子区间二等分，即步长折半。 Step 4. 如果 mmTT221，则停止，输出值 ,1

6、2mT否则，置 m = m+1，h : = h/2 ，转到Step 3。 第32页/共76页例2.3 用变步长的梯形公式计算积分解: 对于 , 定义 f (0) = 1, 首先在区间 0 , 1 上用梯形公式（即步长 h = 1）,求得sin( )xf xx11 (0)(1)0.920 7355,2Tff将 0 , 1 对分, 它的中点函数值 , 则有10.95885102( )f211110.939 7933.222( )TTf6221001TT如果 不成立，则 h := h/2 = 1/2 ，计算10sin xIdxx（精确到10-6）第33页/共76页如此继续下去, 计算结果如下表120

7、2122112)(221kkxfhTT32T6221012TT如果 不成立，则 h := h/2 = 1/4 ，继续计算。)()(22123212xfxfhT)43()41(2212ffhT.9445135. 0第34页/共76页kkT 2kT 2k0123456789100.920 735 50.939 793 30.944 513 50.945 690 90.945 985 00.946 059 60.946 076 90.946 081 50.946 082 70.946 083 00.946 083 1从上表可看出, 将积分区间对分了10次, 求得 I 的近似值为0.9460831 (

8、积分精确值为0.9460831. . . ), 可见收敛速度比较缓慢。为了加速变步长梯形公式的收敛速度，我们采用外推策略。第35页/共76页3. Richardson外推算法 若用一个步长为 h 的函数 I1( h ) 去逼近问题 I , 设其截断误差可表示为 为了提高逼近的精度，选取 q 为满足 的正数, 将上式(1)中的 h 换为 qh , 则有 110Pq其中 )2 , 1( ,)1(kak是与h无关的常数,并且 kkPPPP1210，12(1)(1)(1)112( )kpppkII ha ha ha h(2.1)由 (1) 可知I1(h)逼近I的误差为 )(1PhO。12(1)(1)(

9、1)112()()()()kPPPkII qhaqhaqhaqh(2.2)第36页/共76页211211111)(1)1(121()(1)1kkPPPPPkPPPPPqqqqIaahqqI qhq I hhq(2.2) 式减上式 , 得式(1)两端同乘以 1Pq得 )()1()1(2)1(112111kPkPPPPhahahaqhIIq12(1)(1)(1)112( )kpppkII ha ha ha h (2.1)12(1)(1)(1)112()()()()kPPPkII qhaqhaqhaqh(2.2)第37页/共76页211211111)(1)1(121()(1)1kkPPPPPkPPP

10、PPqqqqIaahqqI qhq I hhq则 I2(h) 逼近I 误差降为2().PO h11112()( )( )1PPI qhq I hIhq,1,111112)1()2()1(2)2(2PPPkkPPPqqqaaqqqaak令其中 是与h无关的常数,则有), 3 , 2(,)2(iaikPkPPhahahahII)2()2(3)2(2232)(，如此继续。第38页/共76页 一般地, 选取 q 为满足 的正数, 由此得到序列10, (1,2,)mPqm1()( )( )1(1,2,)mmPmmmPIqhq IhIhmq则 Im+1 ( h ) 逼近 I 的误差由下面的定理给出。定理

11、2.1 设 I1 ( h ) 逼近 I 的截断误差由下式给出12(1)(1)(1)112( )kpppkII ha ha ha h则 Im+1 ( h )逼近 I 的截断误差为12(1)(1)112( )mmppmmmmmIIhahah其中 是与 h 无关的常数。(1)(1,2,)miaimm这种利用序列Im+1(h)逐步加速去逼近 I 的方法称为Richardson外推算法Richardson外推公式第39页/共76页4. Romberg 算法 Romberg 算法是利用变步长的梯形求积序列 外推加速来逼近积分真值的算法.2kT考虑积分( ).baIf x dx由复化梯形公式有11( )2(

12、)( ) ,2nnkkhITf af xf b现在将 Tn 记为 T1( h ), 则T2n 为 T1(h/2) 且111( )( )2()( ) .2nkkhT hf af xf b10212)(221nkknnxfhTT第40页/共76页 设 f (x) 在区间 a , b 上任意次可微, 根据 Euler-Maclaurin公式有242112( ).kkIT hhhh其中 是与 h 无关的常数。(1,2)kk因为Pm = 2m , 带入上式整理后得易知 Tm+1( h ) 逼近 I 的误差为 O ( h2(m+1) ) ，这种算法称为 Romberg算法。11( )22( )112( )

13、 ( )( )PmPmmmmhTThTh，1(1,2,)4( )2( )41( )mmmmmhTThhmT。1()( )( )1(1,2,)mmPmmmPIqhq IhIhmq则有 选取 利用 Richardson 外推公式, 1,2q 第41页/共76页知T2 ( h ) = Sn，当m =1时，由上式得21141( ),323( )( )hT hTT h则 T2 ( h ) 逼近 I 的误差为 O ( h4 )。12,2()nhTT由 T1 ( h ) = Tn 和241.33nnnSTT故有这是因为 nnTT31342nnkknTxfhT31)(2213410211011021)()(2

14、31)(64nkkknkkxfxfhxfh10121)()(4)(6nkkkkxfxfxfhnS从二分前后两个复化梯形值生成复化 Simpson值 Sn , 将误差 O( h2 )变为O( h4 ), 从而提高了逼近精度第42页/共76页 当 m = 2 时 , 322161( ),15215( )( )hT hTT h则 T3 ( h ) 逼近 I 的误差为 O ( h6 )。22,2( )nhTS2161.1515nnnCSS由T2 ( h ) = Sn ,可以证明 T3 ( h ) = Cn , 故有能从二分前后两个复化Simpson值生成复化Cotes值Cn ,将误差O( h4 )变为

15、O( h6 ), 从而提高了逼近精度。第43页/共76页433641( ),63263( )( )hT hTT h则 T4 ( h ) 逼近 I 的误差为 O ( h8 )。 26416.363nnnRCC上式称为 Romberg 公式, 利用此公式能从二分前后的两个复化 Cotes值生成 Romberg 值 Rn , 且 Rn 逼近 I 的误差为 O ( h8 ) . 32,2( )nhTC由 T3 ( h ) = Cn 知则有令 Rn = T4 ( h )， 当 m = 3 时 ,第44页/共76页 这样可以从变步长的梯形序列 出发, 可逐次 求 得 Simpson 序列 , Cotes

16、序列 , Romberg 序列 ，利用 Romberg 序列 还可以继续外推，但由于继续外推后构成的新的求积序列与原来的序列差别不大, 故通常只外推到 Romberg 序列为止。2kT2kS2kC2kR2kR第45页/共76页 T 1 T 2 T 4 T 8 S 1 S 2 S 4 C 1 C 2 R 10123k2kT12kS22kC32kR其中 表示计算顺序,k代表二分次数. 如果 f ( x ) 在区间 a , b 上充分光滑, 可以证明上表中各列都收敛到积分 所以当同一列中相邻两个数所以当同一列中相邻两个数之差的绝对值小于预给精度之差的绝对值小于预给精度时终止计算时终止计算.( ),a

17、bf x dxRomberg 算法的计算过程241.33nnnSTT2161.1515nnnCSS26416.363nnnRCC第46页/共76页例2.4用Romberg算法计算下列积分 (精确到10-6)10sin xIdxx解: 由变步长梯形公式求得二分 3 次的复化梯形值 T2 , T4 , T8 , 它们精度都很低. 利用 Romberg 算法对其进行加工, 结果列于下表.从此表可以看出, 运用上述二分 3 次的复化梯形值, 采用 Romberg 算法加速三次获得了变步长梯形公式二分 10 次才能获得的结果, 因此加速效果是相当明显的.0123k2kT12kS22kC32kR0.920

18、 735 50.939 793 30.944 513 50.945 690 90.946 145 90.946 086 90.946 083 40.946 083 00.946 083 10.946 083 1241.33nnnSTT2161.1515nnnCSS26416.363nnnRCC241.33nnnSTT2161.1515nnnCSS第47页/共76页3.3 Gauss型求积公式本节讨论在插值节点个数一定的条件下,代数精度最高的插值型求积公式Gauss型求积公式。为了方便，本节考虑加权型积分公式：baxxfxd)()(第48页/共76页则0( ) ( )( )( )( )( ) (

19、).bbnaanbkkakx f x dxx L x dxxlx f xdx设 f (x) 在插值节点 处的函01naxxxb数值为 f (xk), 作 n 次 Lagrange 插值多项式0( )( ) ()nnkkkL xlx f x1 Gauss型求积公式本节考虑带权的积分 其中 为权函数. 若 ，即为通常的积分。 baxxfxd)()(0)(x1)(x第49页/共76页其中 (3.2) ( ) ( ).bkkaAx lx dx上式称为带权的插值型积分公式, Ak 是其求积系数.从而有)3.1(, )(d)()(0nkkkbaxfAxxfx这里取消对积分节点xk是等距的限制。第50页/共

20、76页定理3.1 求积公式(3.1)的代数精度最高不超过2n + 1次.证明：考虑2n+2次多项式：它只在节点xi（i = 0,1,n）处为零，在其它点处均大于零，所以2021)()()(niinxxxxf.0d)()(d)()(21banbaxxxxxfx)3.1(, )(d)()(0nkkkbaxfAxxfx第51页/共76页对于2n+2次多项式 不准确成立，可知(3.1)的其代数精度至多为2n+1。)()(21xxfn而 ， 0)()(0210nkknknkkkxAxfA故插值型的数值积分公式 (3.1) )3.1(, )(d)()(0nkkkbaxfAxxfx第52页/共76页定义 3

21、.1 若插值型求积公式(3.1) 13.(, )(d)()(0nkkkbaxfAxxfx度,则称该求积公式是 Gauss 积分公式, 相应的求积节点 xk ( k = 0, 1, , n ) 称为Gauss点。其中bakkxxlxAd)()( ，具有 2n +1 次代数精1.1 Gauss积分公式和Gauss点为了确定Gauss点和求积系数Ak，需要利用正交多项式理论。第53页/共76页定义 3.2 给定区间 a, b 和权函数 (x) 及多项式序列0,(,)( )( )( )0,bjkjkakjkggx gx gx dxAjk0( )(0,1,2,),kjkjjgxa xk其中首项系数0,k

22、a 若 gk(x) 满足则称 为在区间a ,b上带权 (x) 的正交多项式序列, gk(x) 称为k 次正交多项式. 0)(kkxg定理3.2 n次正交多项式 gn(x)与任意次数不超过n-1的多项式 P(x)在区间a, b上带权(x)正交 。 第54页/共76页1. 1 Gauss积分公式定理 3.3 带权Gauss 求积公式(3.1)中的求积节点 xk (k = 0, 1, . . ., n) 是Gauss 点的充分必要条件是 与任意次数不超过 n 的多项式 P(x) 均在区间 a, b带权 正交, 即 。从而在a , b上带权 的 n + 1次正交多项式 Pn+1(x) 的零点即为Gau

23、ss 点。10( )()nnkkxxx( )x( )x) 13.(, )(d)()(0nkkkbaxfAxxfx0d)()()(1banxxxPx第55页/共76页111111d)()()()(d)()(xxxxxxxxlxAknknkk111111d)()()()(xxAxxxAxknnknn)3 . 3(d)()()()(1111xxPxxxPxknkn设求积节点xk (k = 0, 1, . . ., n)是正交多项式 的零点，其首项系数为 ，则 )(1xPn1nA.)()()(01111nkknnnnxxAxAxP于是求积系数1.1 Gauss 积分公式第56页/共76页定理3.4 G

24、auss型求积公式(3.1)的求积系数Ak) 13.(, )(d)()(0nkkkbaxfAxxfx满足下面的性质：0( ).nbkakAx dx., 1 , 00d)()(2njxxlxAbajj，取f (x) = 1.d)()()(202bajnkkjkjxxlxxlAA下面讨论Gauss型求积公式的截断误差（余项）。第57页/共76页1.2 Gauss积分公式的余项, 收敛性与稳定性定理 3.5 设函数 f C2n+2a , b, 则 Gauss 求积公式的余项为 0(22)21( )( ) ( )()( )( )( ),( , ).(22)!nbkkaknbnaR fx f x dxA

25、 f xfxx dxa bnGauss积分公式的余项第58页/共76页Gauss 求积公式的数值稳定性对于某一给定的算法，原始数据的误差（如舍入误差）为，且假设在运算过程中的其它误差都是由引起的，如果误差在一定条件下能够得到控制，即数值计算结果的误差至多是原始数据误差的同阶无穷小量，则称该算法是数值稳定的；否则称该算法是数值不稳定的。 第59页/共76页Gauss 求积公式的数值稳定性 设 f (xk) 的近似值为() (0,1, ),kf xkn记0max |.kk n 由 Gauss 求积公式和Ak 0, 则有误差估计00000()()|( ).| |nnnnkkkkkkkkkkkknbk

26、akA f xA f xAAAx dxC()(0,1,),),kkkf xf xkn令第60页/共76页其中( )baCx dx是一个大于零的常数. 由此可知 Gauss 求积公式是数值稳定的.00000()()|( ).| |nnnnkkkkkkkkkkkknbkakA f xA f xAAAx dxC第61页/共76页定理 3.6 对任意的 f Ca , b, 则 Gauss 求积公式均收敛, 即有0lim()( ) ( ).nbkkankA f xx f x dx对于 Gauss 求积公式的收敛性, 有如下的定理第62页/共76页2 几种常用的Gauss型求积公式Gauss-Legend

27、re求积公式 Legendre多项式 是定义在区间 -1 , 1上, 权函数为 的正交多项式。其标准形式为:02( )1,1( )(1)(1,2,) .2!nnnnnP xdP xxnn dx1)(x第63页/共76页可以证明, Legendre 多项式有下列递推关系：0111( )1,( ),(1)( )(21)( )( ),1,2,nnnP xP xxnPxnxP xnPxn22334245351( )(31)21( )(53 )21( )(35303)81( )(637015 )8P xxP xxxP xxxP xxxx 由上式可推出：当 n 为偶数时, Legendre 多项式 Pn(

28、x) 为偶函数；当 n 为奇数时, Legendre 多项式 Pn(x) 为奇函数。第64页/共76页所以，当n =0 时，一次勒让德(Legendre)多项式 P1( x ) = x的零点（Gauss点）x0 = 0，取其为求积节点，由。)0(2)(d)(0011fxfAxxf1111d)()()(xxPxxxPAknknk得A0 = 2。从而得到一点Gauss-Legendre 求积公式第65页/共76页当n = 1时，二次勒让德（Legendre）多项式 ) 13(21)(22xxP它有两个零点(Gauss点) , 取它们为求积节点，由（6.5.9）确定出A0 = A1 = 1。从而得到二点Gauss-Legendre 求积公式 33。)33()33()(11ffdxxf)3 . 3(d)()()(1111xxPxxxPAknknk第66页/共76页常见Gauss-Lege