数值积分课件(《计算方法》)

1、数值计算方法2022-5-11第7章 数值积分 l1 插值型求积公式插值型求积公式 l2 复化求积公式复化求积公式l3 龙贝格龙贝格(Romberg)求积方求积方法法数值计算方法2022-5-12 1 插值型求积公式插值型求积公式 l在一元函数的积分学中,我们已经熟知,若函数f(x)在区间a, b上连续且其原函数为F(x) ,则可用牛顿莱布尼兹公式 ( )( )( )baf xF bF a(71)来求定积分。数值计算方法2022-5-13 公式(71)虽然在理论上或在解决实际问题中都起了很大的作用,但它并不能完全解决定积分的计算问题。因为定积分的计算常常会碰到以下三种情况:l (1)被积函数f

2、(x)的原函数F(x)不易找到。许多很简单的函数,例如 2sin1, , lnxxexx其原函数都不能用初等函数表示成有限形式。数值计算方法2022-5-14l (2)被积函数f(x)没有具体的解析表达式。其函数关系由表格或图形表示,无法求出原函数。 41badxx 其被积函数 的原函数就比较复杂,从数值计算角度来看,计算量太大。411xl(3)尽管f(x)的原函数能表示成有限形式但其表达式相当复杂。例如定积分 数值计算方法2022-5-15 图 7.1 l 如图7.1,若用左矩形近似地代替曲边梯形,则得到左矩形公式 (72)( )() ( )baf x dxb a f a数值计算方法2022

3、-5-16l同样可得到右矩形公式：( )() ( )baf x dxba f b(73) 数值计算方法2022-5-17图 7.2l 如图7.2,若用梯形的面积近似地代替曲边梯形的面积,则得到计算定积分的梯形公式 ( ) ( )( )2babaf x dxf af b(74) 数值计算方法2022-5-18l 如图7.3,若用抛物线代替曲线f(x),则可得到抛物线公式(或辛普生公式)( ) ( )4 ()( )2babaabf x dxf aff bb(75) 图7.3 数值计算方法2022-5-19此外此外,众所周知的众所周知的梯形公式梯形公式:l I(f)(b-a)f(a)+f(b)/2和

4、和 Simpson公式公式:l I(f)(b-a)f(a)+4f(a+b)/2)+f(b)/6则分别可以看作用则分别可以看作用 a, b, c=(a+b)/2, 三点三点高度的加权平均值高度的加权平均值 f(a)+f(b)/2 和和 f(a)+4f(c)+f(b)/6作为作为平均高度平均高度f()的近似值的近似值.数值计算方法2022-5-110 更一般地更一般地,取区间取区间a,b内内n+1个点个点 xi,(i=0,1, 2,n) 处的高度处的高度f(xi) (i=0,1,n)通过通过加权平均加权平均的方法的方法近似地得出平均高度近似地得出平均高度f(),这类求积方法称为这类求积方法称为机械

5、求机械求积积:)()()(0ibaniixfabdxxf数值计算方法2022-5-111 或写成或写成:数值积分公式数值积分公式求积系数求积系数 求积节点求积节点 )()(0kbankkxfAdxxf(1)数值计算方法2022-5-112记记)2()()(0knkknxfAfI)3(, )()()()()(0bankkknxfAdxxffIfIfR称称(2)为数值求积公式为数值求积公式,(3)为求积公式余项为求积公式余项(误差误差). 构造或确定一个求积公式，要讨论解决的问构造或确定一个求积公式，要讨论解决的问题有题有(i) 确定求积系数确定求积系数Ak和求积节点和求积节点xk ；(ii)求积

6、公式的误差估计和收敛性求积公式的误差估计和收敛性 为了构造形如式为了构造形如式(2)的求积公式的求积公式,需要提供一需要提供一种种判定求积方法精度高低准则判定求积方法精度高低准则数值计算方法2022-5-113求积公式的代数精度求积公式的代数精度定义定义1 称求积公式称求积公式(2)具有具有m次次代数精度代数精度,如果它满如果它满足如下两个条件足如下两个条件: (i)对所有次数对所有次数 m次的多项式次的多项式 ,有有 (ii)存在存在m+1次多项式次多项式 ,使得使得)(xPm0)()()(mnmmPIPIPR)(1xPm0)()()(111mnmmPIPIPR定义定义1中的条件中的条件(i

7、),(ii)等价于等价于:0)()()0(, 0)()()()(1mknkkxRiimkxIxIxRi数值计算方法2022-5-114插值型求积公式插值型求积公式 在积分区间在积分区间a,b 上上取取n+1个节点个节点xi，i=0,1,2,n,作作f(x)的的n次代数插值多项式次代数插值多项式（拉格朗（拉格朗日插值公式）日插值公式）:则有则有 为插值余项为插值余项于是有于是有njjjnxfxlxL0)()()()()()(xRxLxfnn)()!1()()(1)1(xwnfxRnn数值计算方法2022-5-115取取称称(4)式为插值型求积公式式为插值型求积公式,其中其中求积系数求积系数Ak由

8、由(5) 式确定式确定. bajnjbajbanbanbadxxRxfdxxldxxRdxxLdxxf)()()()()()(0(4)(5) babaknkkdxxlxfdxxf)()()(0Ak0()()nbikaikiikxxAdxxx由由 节点节点 决定，决定，与与 f(x) 无关。无关。数值计算方法2022-5-116数值计算方法2022-5-117推论推论1 求积系数满足求积系数满足:0(1)0 ( )() ( )( )()()(1)!nbbkknaaknnbxkakR ff x dxA f xf xLx dxfxxdxn误误 差差定理定理1 形如形如 的求积公式至少有的求积公式至少

9、有 n 次代数精度次代数精度 该该公式为公式为插值型插值型（即：（即： ） nkkkxfA0)( bakkdxxlA)(abAnjj0数值计算方法2022-5-118l 现用第六章介绍的插值多项式Pn(x)来代替被积函数f(x),即有 1.1 1.1 牛顿牛顿柯特斯公式柯特斯公式(NewtonCotes)(NewtonCotes) ( )( )bbaaf x dxx dx( )( )bbnaaf x dxP x dx取节点为等距,即 a=x0 x1xn=b l 建立数值积分公式最基本的思想是选取一个既简单又有足够精度的函数(x), 用(x)代替被积函数f(x),于是有 数值计算方法2022-5

10、-119l 利用拉格朗日插值多项式10,0,1,2,10,1,2,kkibahxx knnxxihin( )( )( )nnf xpxRx(76) 其中 000(1)1( )( )()( )( )( )( , )(1)!nnnkniiiiikikk innnxxP xl x yyxxfR xxa bn (77) 数值计算方法2022-5-120l 这里yi=f(xi),对式(76)两边积分得 (1)1000( )( )( )1( )( )(1)!( )bbbnnaaannbbnkinaaikikk iniinif x dxpx dxR x dxxxdx yfx dxxxna yRf 数值计算方

11、法2022-5-121为牛顿牛顿柯特斯柯特斯(Newton-Cotes)(Newton-Cotes)求积公式求积公式,Rn(f)为牛顿柯特斯求积公式的余项。0(1)1 7-81( )( )( ) 7-9(1)!nbkiakikk ibnnnaxxadxxxRffx dxn （）（）我们称 0( ) 7-10nbiiaif x dxa y（）数值计算方法2022-5-122( )000()nnbnnkiiakkikk ik ixxbaskadxdsba cxxnik令 x=x0+sh , 0sn dx=hds=(b-a)/nds(711) 0,1,2,in()00001(1)()!()!nnni

12、kkininnkkiskCdsniksk dsinin数值计算方法2022-5-123(1)12(1)00( )( )( )(1)!( )(),( , )(1)!nbnannnnjfR fwx dxnhfsj dsa bnNewton-Cotes公式的误差为公式的误差为:与与x有关有关注意注意:由由(7-11)式确定的式确定的Cotes系数只与系数只与i和和n有关有关,与与f(x)和积分区间和积分区间a,b无关无关,且且满足满足:( )01nniiC(7-9)数值计算方法2022-5-124l 称Ci(n)为柯特斯求积系数。很显然,当n=1时,可算得1(1)001(1)101(1)212Csd

13、sCsds 此时式(710)为 ( ) ( )( )2babaf x dxf af b(712) 这是梯形公式梯形公式。数值计算方法2022-5-125 l 当n=2时,可得2(2)002(2)102(2)2011(1)(2)4614(2)2611(1)46CssdsCs sdsCs sds 于是 ( ) ( )4 ()( )62babaabf x dxf aff b(713) 这是抛物线抛物线(Simpson)公式公式。数值计算方法2022-5-126l 当n=3时,3(3)003(3)103(3)203(3)3011(1)(2)(3)18813(2)(3)6813(1)(3)6811(1)

14、(2)188CsssdsCs ssdsCs ssdsCs ssds 代入(710)式得到求积公式 0123( ) ( ) 3 ( ) 3 ( )( ) 7-148bab af x dxf xf xf xf x（）数值计算方法2022-5-127类似地可分别求出n=4,5,时的柯特斯系数,从而建立相应的求积公式。具体结果见表71。 l 从表中可以看出,当n7时,柯特斯系数为正;从n8开始,柯特斯系数有正有负。因此,当n8时,误差有可能传播扩大,牛顿 柯特斯求积公式不宜采用。l 柯特斯系数Ci (n) 仅与n和i有关,与被积函数f(x)无关,且满足 ( )01nniiC(715) 事实上,式(71

15、0)对f(x)=1是准确成立的。数值计算方法2022-5-128 表 71 数值计算方法2022-5-129l 定理定理 当阶数当阶数n为偶数时为偶数时, Newton-Cotes公式公式(8)至少具有至少具有n+1次代数精度次代数精度.证明证明 只需验证当只需验证当n为偶数时为偶数时,Newton-Cotes公式公式对对f(x)=xn+1的余项为零的余项为零.由于由于f(x)=xn+1,所以所以f(n+1)(x)=(n+1)! .由式由式(7-9)得得 nnjndtjthfR002)()(引进变换引进变换t=u+n/2,因为因为n为偶数为偶数,故故n/2为整数为整数,于是有于是有 2202)

16、2()(nnnjndujnuhfR据此可断定据此可断定R(f)=0,因为上述被积函数是个奇函数因为上述被积函数是个奇函数.数值计算方法2022-5-130Newton-Cotes公式的数值稳定性公式的数值稳定性 现在讨论现在讨论舍入误差舍入误差对计算结果产生的影响对计算结果产生的影响.设设用公式用公式 近似计算积分近似计算积分 时时,其中计算函数值其中计算函数值f(xj)有误差有误差j (j=0,1,2,n).设设计算计算Cj(n)没有误差没有误差,中间计算过程中的舍入误差也中间计算过程中的舍入误差也不考虑不考虑,则在式则在式(10 )的计算中的计算中,由由j引起的误差为引起的误差为njjnj

17、nxfCabfI0)()()()(badxxffI)()(10)数值计算方法2022-5-131njjnjnjjjnjnjjnjnCabxfCabxfCabe0)(0)(0)()()()()()(如果如果Cj(n)都是正数都是正数,并设并设|max0jnj)(|)(|0)(abCabenjnjn故故en是有界的是有界的,即由即由j引起的误差受到控制引起的误差受到控制,不超过不超过的的(b-a)倍倍,保证了保证了数值计算的稳定性数值计算的稳定性.而当而当n7时时,Cj(n)将出现将出现负数负数,njnjC0)(|保证数值稳定性保证数值稳定性. 因此高阶公式不宜采用因此高阶公式不宜采用,有实用有实

18、用价值的仅仅是几种价值的仅仅是几种低阶的求积公式低阶的求积公式.将随将随n增大增大,因而因而不能不能则有则有数值计算方法2022-5-132解解 利用梯形公式10.5xdx10.510.5(0.51)20.4267767xdx利用抛物线公式10.510.5(0.540.751)60.43093403xdx例1 试分别用梯形公式和抛物线公式计算积分:原积分的准确值 31120.50.520.430964413xdxx数值计算方法2022-5-133l 现对牛顿柯特斯求积公式所产生的误差作一个分析。由式(79),牛顿柯特斯求积公式的余项为 (1)11( )( )( )(1)!bnnnaRffx d

19、xn 1.2 1.2 误差估计误差估计l 易知,牛顿柯特斯求积公式(710)对任何不高于n次的多项式多项式是准确成立的。这是因为 f(n+1)()0 故 Rn(f)0数值计算方法2022-5-134l 牛顿牛顿柯特斯求积公式的代数精确度至少为柯特斯求积公式的代数精确度至少为n n。通常在基点个数相等的情况下,代数精确度愈高,求积公式就愈精确。l 一般说来,若某个求积公式对于次数不高于m的多项式都准确成立(即Rn(f)0),而对于某一次数为m+1的多项式并不准确成立,则称这一求积公式的代数精确度代数精确度为m。l 定理定理1 (梯形公式的误差)设f(x)在区间a, b上具有连续的二阶导数,则梯形

20、求积公式的误差为31()( )( ),( , ) 7-1612b aR ffa b （）数值计算方法2022-5-135由于 1(x)=(x-a)(x-b)11( )( )()()2baRffxa xb dx证 由式(79)知,梯形公式的余项为 l 1(x)在区间(a, b)内不变号,f()是x的函数且在a,b上连续,故根据积分第二中值定理知,存在某一(a, b)使131( )( )()()21()( )12baRffxaxb dxbaf 数值计算方法2022-5-136l 定理定理2 2 (抛物线公式的误差)设f(x)在a, b上有连续的四阶导数,则抛物线公式的误差为(717) 5(4)2(

21、)( )( ),( , )2880baRffa b 数值计算方法2022-5-137(1)(1)0111,22CCn = 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 = 4:

22、 Cotes Rule, 代数精度代数精度 = 5,)(9458)6(7 fhfR 数值计算方法2022-5-138复合求积公式复合求积公式高次插值有高次插值有Runge 现象现象, ,高阶高阶Newton-Cotes公式会出现公式会出现数值不稳定数值不稳定, ,低阶低阶Newton-Cotes公式有时又不能满足精度公式有时又不能满足精度要求要求. .解决这个矛盾的办法是将积分区间解决这个矛盾的办法是将积分区间 a, ,b 分成若干小区间分成若干小区间, ,在每个小区间上在每个小区间上用低阶求积公式计算用低阶求积公式计算, ,然后将它们加起来然后将它们加起来, ,这就是这就是复合求积方法复合求

23、积方法. .数值计算方法2022-5-1392 复合求积公式 l 2.1 复合梯形公式 对于定积分(71),将积分区间a, b分成n个相等的子区间xi,x i+1,这里步长1, 0,1,1iibahxxinn在每一个子区间 xi,x i+1 上使用梯形公式,则 数值计算方法2022-5-140 复化梯形公式积分法复化梯形公式积分法数值计算方法2022-5-141相加后得 1311( ) ( )()( )212( ,)iixiiixiiihhf x dxf xf xfx x(718) 110311100( )( ) ( )()( )212iinbxaxinniiiiif x dxf x dxhh

24、f xf xf(719) 若f(x)在a,b上连续,由连续函数的介值定理,存在某一(a,b)使得10( )()niiffn数值计算方法2022-5-142因而 331032( )()1212()()12niihnhffbafn 于是得到复合梯形公式 1013( )12( ) ()()2()2()( )( )12nbniainhf x dxf xf xf xbaRffn (721)其余项为 数值计算方法2022-5-143例2 若用复合梯形公式计算积分 10 xe dx3()12()()( )12nbaRffn 则 当0 x1时,有 因为 又 解 由余项(721)式( )12( )( )( ),

25、 1, ( )( )12xnf xf xfxeb afxeeRfn 问积分区间要等分多少才能保证有五位有效数字?数值计算方法2022-5-144l 由于原积分的准确值具有一位整数,因此要使近似积分值有五位有效数字,只需取n满足 4224110122610enne 两边取对数得 整理后得到 lg2lg4lg4lglglg1.82662bneebn0811.67n 取 n=68.数值计算方法2022-5-145l 类似复合梯形公式的做法,把区间a,b分成n个相等的子区间x2i,x2i+2(i=0,1,n-1),设每个子区间上的中点为x2i+1(i=0,1,n-1),且 1,0,1,212iibah

26、xx inn2225(4)22122222( ) ()4 ()()( ),390(,)iixiiiixiiihhf x dxf xf xf xfxx(722) 2.2 复合抛物线公式l 在每一个子区间 x2i ,x2i+2 上利用抛物线公式得数值计算方法2022-5-146 复化复化Simpson公式积分法公式积分法数值计算方法2022-5-147l相加后得 22210511(4)2212200( )( ) ()4 ()()( )390iinbxaxinniiiiiif x dxf x dxhhf xf xf xf(723) 数值计算方法2022-5-148l 若f(4)(x)在a, b上连续

27、,则5551(4)(4)(4)40()( )( )( ), 90902880 ( , )niiihnhbafffna b 从而得到复合抛物线公式 1202212200( ) ()()4()2()3nnbniiaiihf x dxf xf xf xf x(724) l 其余项为 5( )(4)24()( )( )2880nbaRffn (725)数值计算方法2022-5-149 图 7.4 复合抛物线公式框图数值计算方法2022-5-150sin( )xfxx的数据表,分别用复合梯形公式和复合抛物线公式计算 10sinxIdxxl 例3 已知函数 x f (x)0 11/8 0.99739782

28、/8 0.98961583/8 0.97672674/8 0.95885105/8 0.93615566/8 0.90885167/8 0.87719251 0.8414709 数值计算方法2022-5-151l 解 用复合梯形公式,这里10.1258h 10sin0.125 (0)2 (0.125)(0.25)2(0.375)(0.5)(0.625)(0.75)0.94569086xdxfffxffff数值计算方法2022-5-152l 用复合抛物线公式可得10sin0.125 (0)4 (0.125)(0.375)3(0.625)(0.875)2 (0.25)(0.5)(0.75)(1)0

29、.946083305xdxfffxffffff 比较上面两个结果比较上面两个结果T8和和S4,它们都需要提供它们都需要提供9个个点上的函数值工作量基本相同点上的函数值工作量基本相同,然而精度却差别很大然而精度却差别很大. 同积分的同积分的准确值准确值I(f)=0.9460831比较比较,复化梯形法复化梯形法的结果的结果T8=0.9456909只有只有两位有效数字两位有效数字, 而复化而复化Simpson法的结果法的结果S4=0.9460832却有却有六位有效数字六位有效数字.数值计算方法2022-5-153 复化梯形公式：复化梯形公式：,(0,., )ibahxaihinn在每个在每个 上用梯

30、形公式：上用梯形公式：,1iixx111( ) ( )(),0,.,12iixiiiixxxf x dxf xf xin11( )2()( )2nkihf af xf b110( ) ( )()2nbiiaihf x dxf xf x= Tn1321002() ()()1212()( ),( , )12niniiifhhR ffbanhba fa b /*中值定理中值定理*/数值计算方法2022-5-154,(0,1,.,2 )2jbahxaj hjnn222121( ) ()4 ()()3njjjjhI ff xf xf x22jx21jx2 jx222121( ) ()4 ()()3nnj

31、jjjhSff xf xf x 复化复化 Simpson 公式：公式：121211( ) ( )4()2()( )3nnbjjajjhf x dxf af xf xf b数值计算方法2022-5-1552.3 变步长公式l 前面介绍的复合梯形公式和复合抛物线公式的步长都是预先确定的。它的主要缺点是事先很难估计出n的大小(或步长h的大小),使结果达到预先给定的精度。在实际计算中,我们常常借助于计算机来完成积分步长h的自动选择,即采用变步长求积公式。具体地讲,就是将步长逐次折半,反复利用复合求积公式,直到满足精度要求为止。 数值计算方法2022-5-156l 下面介绍变步长复合抛物线公式(变步长复

32、合梯形公式留给读者作为练习)。 121221220 ()4 ()()3mmmkkkkhSf xf xf x(726) 其中 2mmbah再把每个子区间分成两半,用 21122mmmbahh 逐次将区间a,b分成2,4,2m等分,并按复合抛物线公式逐次计算积分得到S1,S2,Sm,而数值计算方法2022-5-157 作步长,按复合抛物线公式计算出积分的近似值S2m。对于相邻两次的积分近似值Sm、S2m,考察 222mmmmmSSdSSS当S2m1 当S2m1 (727) 设预先给定的精度为,若 d 则以S2m作为所要求的积分近似值,否则继续将区间分半,利用复合抛物线公式求积分,直到满足预给的精度

33、为止。数值计算方法2022-5-158 图 7.5 变步长复合抛物线公式数值计算方法2022-5-159 图7.5 变步长复合抛物线公式数值计算方法2022-5-1603 龙贝格(Romberg)积分方法 l我们已经知道,当被积函数f(x)在区间a,b上连续时,要使得复合梯形公式或复合抛物线公式比较精确地代替定积分( )baf x dxl可将分点(即基点)加密,也就是将区间a,b细分,然后利用复合梯形公式或复合抛物线公式求积。数值计算方法2022-5-161l若用Tm表示把a,b作m等分并按复合梯形公式求积的结果,将每一小段再对分,令新的小段的长h=h/2,则T2m与Tm之间有如下关系: 21

34、12(21) mmmkTThf akh(528) 其中 数值计算方法2022-5-162l 另外,若用Sm表示把a,b分成m(偶数)个小段按复合抛物线公式计算的结果,那么只要把Sm中的m改为2m,h改为h就有20122221201222212024222(4224)3(244442)3(222)3MmmmmmmmmhSffffffhffffffhfffff 从Tm的定义可得到关系式 22441mmmTTS(529) 数值计算方法2022-5-163l我们再举一个计算上半单位圆面积的例子(它的准确面积为/2)。现用内接正多边形的逼近方法来计算。 l如图5.6,图(a) 、 图(b)是用同样的内接

35、正多边形计算上半单位圆的面积。图(a)是用梯形方法计算其面积,图(b)是用三角形方法计算其面积。 数值计算方法2022-5-164 图 5.6 数值计算方法2022-5-165l 设正多边形边数为n=2k,则由图(b)利用三角形公式算得面积为 ( )201235111133524(1)350112sincossin2sin242112()()23! 25! 222()()23!25!222()()23!25!2kkkkkkkkkkkkkkkknnTnnnT同理 数值计算方法2022-5-166l 如果组合一下,就会得到更精确的结果,即 (1)( )( )00154441225!(2 )kkkk

36、TTT同理 (2)(1)(1)00151 4441225!(2)kkkkTTT数值计算方法2022-5-167l 再以类似方法组合得 2(1)( )( )11164411()2(2 )kkkkTTTO 这样继续下去,其值越来越接近上半单位圆面积/2。这种方法可以用到计算定积分 ( )baf x dx数值计算方法2022-5-168l 为了推广公式(529)和上述计算上半单位圆面积的组合方法,我们引进龙贝格求积算法。l 龙贝格求积算法本来是利用所谓外推法构造出的一种计算积分的方法。为了避免从外推引入而带来理论上的麻烦,我们将直接从构造一个T数表开始。l 首先将a,b依次作20,21,22,等分,

37、记,0,1,2,2iibahi数值计算方法2022-5-169l 按复合梯形公式(520)算得的值相应地记为T(k)0(k=0,1,2,);把按式(529)算得的S2m依次记为T(k)1(k=0,1,2,崐),而这每一个S2m又理解为由T2m与Tm的线性组合得到的改进值,即(1)( )( )0014,0,1,2,41kkkTTTk 我们可按照类似的方法继续进行改进,也即由S2m与Sm的线性组合得到改进值,依次记为T(k)2(k=0,1,2,),即 2(1)( )( )10224,0,1,2,41kkkTTTk数值计算方法2022-5-170l 这样就可构造出一个数表 (5-30)数值计算方法2022-5-171l 其中除第0列(即最左一列)的T(k)0是按复合梯形公式计算外,其余各列都按下述规则(对m) (1)( )( )111,2