数值分析方法第六章

1、Chapter 5 数值微分与数值积分2先看一个实例:已知20世纪美国人口的统计数据为(单位:百万)年份 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990人口 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 试计算美国20世纪的年增长率则人口的增长率为时刻的人口为若记),(txt)()(txdtdxtr如何求dtdx 1 数值微分3 在实际问题中，往往会遇到某函数在实际问题中，往往会遇到某函数f(x)f(x)是是用表格用表格表示的表示的, ,用通常的导数定义无法求导用通常的导数

2、定义无法求导, ,因此要寻求其他因此要寻求其他方法近似求导。常用的数值微分方法有方法近似求导。常用的数值微分方法有: :一一. . 运用差商求数值微分运用差商求数值微分2 2 运用插值函数求数值微分运用插值函数求数值微分三三. . 运用样条插值函数求数值微分运用样条插值函数求数值微分四四. . 运用数值积分求数值微分运用数值积分求数值微分4运用差商求数值微分运用差商求数值微分分析分析：微分和积分是一对互逆的数学运算。：微分和积分是一对互逆的数学运算。复习复习：导数：导数 是差商是差商 当当 时的极限，在精度要求不高的情形下，可时的极限，在精度要求不高的情形下，可以简单地取差商作为导数的近似值以

3、简单地取差商作为导数的近似值( )fa()( )f ahf ah0h 基本原理：用差商代替微分基本原理：用差商代替微分5可可用用差差商商来来逼逼近近导导数数充充分分小小时时当当,h处在点根据导数定义x,hhxfhxfhhxfxfhxfhxfxfhhh2)()(lim)()(lim)()(lim)( 0006( )()( )f af ahfah()()( )2f ahf ahfah向后差商近似计算：向后差商近似计算：()( )( )f ahf afah向前差商近似计算：向前差商近似计算：中心差商近似计算：中心差商近似计算： 差商公式差商公式7hxfhxfxf)()()( 000由Taylor展开

4、hxxfhxhfxfhxf002000),( ! 2)( )()(误差)()( ! 2)()()( )(000hOfhhxfhxfxfxR向前差商8hhxfxfxf)()()( 000由Taylor展开hxxfhxhfxfhxf002000),( ! 2)( )()(误差)()( ! 2)()()( )(000hOfhhhxfxfxfxR向后差商9hhxfhxfxf2)()()( 000由Taylor展开23000010102300002020()()()()( ),2!3!()()()()(),2!3!hhf xhf xhfxfxfxxhhhf xhf xhfxfxfxhx误差)()( 6)

5、( )( 12 2)()()( )(22212000hOfhffhhhxfhxfxfxR中心差商10分析分析 1、步长过大则截断误差显著，但如果步长太小又会导、步长过大则截断误差显著，但如果步长太小又会导致舍入误差的增长。致舍入误差的增长。 2、在实际计算时，在保证截断误差满足精度要求的前、在实际计算时，在保证截断误差满足精度要求的前提下选取尽可能大的步长提下选取尽可能大的步长 3、事先给出一个合适的步长往往是很困难的。通常在、事先给出一个合适的步长往往是很困难的。通常在在在变步长变步长的过程中实现步长的的过程中实现步长的自动选择自动选择。11主要方法：插值多项式，讨论函数主要方法：插值多项式

6、，讨论函数f(x)导数的近似值求法。导数的近似值求法。 插值型求导公式插值型求导公式)()(xPxfn 定义：若函数定义：若函数f(x)在节点在节点 处的函数处的函数值已知，作值已知，作f(x)的的n次插值多项式次插值多项式 ，并用，并用 近似代近似代替替f(x)，即，即) 1, 2 , 1(nixi)(xPn由于由于 是多项式，容易求其导数，故对应于是多项式，容易求其导数，故对应于f(x)的每一个插值多项式的每一个插值多项式 ，建立一个数值微分公式，建立一个数值微分公式)(xPn)(xPn( )( )nfxP x这样建立起来的数值微分公式，称为这样建立起来的数值微分公式，称为插值型数值微分公

7、式。)(xPn12(1)(1)11( )( )( )( )( )( )(1)!(1)!nnnnnxfdff xP xxnndx其中其中 之间，上式第一项是之间，上式第一项是 的近的近似值。似值。nxxx,10在( )ifx即使即使 与与f(x)的函数值处处相差甚微，但两个的函数值处处相差甚微，但两个函数的导数在某些点上的值仍可能有很大的差异函数的导数在某些点上的值仍可能有很大的差异所以要对误差进行分析所以要对误差进行分析)(xPn13125102050100025020015010050300annealing time t, 1000 sSeebeck Coefficient a a, m

8、mV/Ka a 010203040dtda adtda a误差放大 !)exp()(nKtBAt a a14(1)1( )( )( )( )(1)!nininiffxP xxn如果只计算结点处的导数值如果只计算结点处的导数值( )( )inifxP x数值微分公式的余项为数值微分公式的余项为1510001101( )( )( )( )2, , ( )( )( )( )2f xf xhf xfhx xf xf xhf xfh（n=1）011010110( )()()xxxxp xf xf xxxxx100101( )( )( )( )( )( )f xf xf xhf xf xf xh16200

9、122102220121() 3 ()4 ( )()( )231( )()()( )261() () 4 ( ) 3 ()( )23hf xf xf xf xfhhf xf xf xfhhf xf xf xf xfh02,x x类似地，可以得到类似地，可以得到n=2的带余项的三点公式的带余项的三点公式172(4)00121222(4)101222(4)20121221() ()2 ()() ()()61() ()2 ()()( )121() ()2 ()() ()()6hfxf xf xf xhffhhfxf xf xf xfhhfxf xf xf xhffh 02,x x类似地，可以得到类似

10、地，可以得到n=2的带余项的二阶三点公式的带余项的二阶三点公式18利用三次样条插值函数构造数值微分公式利用三次样条插值函数构造数值微分公式对于给定函数表对于给定函数表 和适当地边界条件，有三次样条插值函数和适当地边界条件，有三次样条插值函数S(x)，并用并用S(x)近似代替近似代替f(x)，即，即由于由于S(x)是一个分段三次多项式，在各子区间是一个分段三次多项式，在各子区间 上容易求出其导数，故建立一个数上容易求出其导数，故建立一个数值微分公式值微分公式.,)()(baxxSxf,1iixx), 2 , 1(ni192211111( ) ( )()() 2,()6iiiiiiiiiiifxS

11、xM xxMxxhhf xxMM( )iiSxM注意注意1111( ) ( )()()(,1,2,)iiiiiiiifxSxMxxMxxhxxxi数值微分公式数值微分公式20例：例： 利用函数利用函数 在节点在节点 上的函数值和边界条件上的函数值和边界条件S(-1)=0.0740, S(1)=-0.0740构造三次样条插值函数构造三次样条插值函数S(x)，并用它来计算，并用它来计算f(x)和和f(x)在下列点在下列点处的近似值。处的近似值。22511)(xxf1 0.1 (0,1,20)ixi i )100, 1 , 0(02. 01kkxi2122MATLAB中的数值微分中的数值微分数值微分

12、的实现数值微分的实现在在MATLAB中，没有直接提供求数值导数的函数，只有中，没有直接提供求数值导数的函数，只有计算向前差分的函数计算向前差分的函数diff，其调用格式为：，其调用格式为：DX=diff(X)：计算向量：计算向量X的向前差分，的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,n-1。DX=diff(X,n)：计算：计算X的的n阶向前差分。阶向前差分。diff(X,2)=diff(diff(X)。23例例 设设x由由0,2间均匀分布的间均匀分布的10个点组成，求个点组成，求sinx的的13阶差分。阶差分。命令如下命令如下：X=linspace(0,2*pi,10);Y=

13、sin(X);DY=diff(Y); %计算计算Y的一阶差分的一阶差分D2Y=diff(Y,2); %计算计算Y的二阶差分的二阶差分，也可用命令也可用命令diff(DY)计算计算D3Y=diff(Y,3); %计算计算Y的三阶差分的三阶差分，也可用也可用diff(D2Y)或或diff(DY,2)24例例 用不同的方法求函数用不同的方法求函数f(x)的数值导数，并在同一个坐标系中的数值导数，并在同一个坐标系中做出做出f(x)的图像。的图像。程序如下：程序如下：f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2);g=inline(3*x.2+4*x-1

14、)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5);x=-3:0.01:3;p=polyfit(x,f(x),5); %用用5次多项式次多项式p拟合拟合f(x)dp=polyder(p); %对拟合多项式对拟合多项式p求导数求导数dpdpx=polyval(dp,x); %求求dp在假设点的函数值在假设点的函数值dx=diff(f(x,3.01)/0.01; %直接对直接对f(x)求数值导数求数值导数gx=g(x); %求函数求函数f的导函数的导函数g在假设点的导数在假设点的导数plot(x,dpx,x,dx,.,x,gx,-); %作图作图252 数值

15、积分数值积分在高等数学中，牛顿在高等数学中，牛顿-莱布尼兹公式莱布尼兹公式F(x)是被积函数是被积函数f(x)的原函数的原函数完美地解决了定积分问题。完美地解决了定积分问题。)a(F)b(F)x(Fdx)x(fbaba 26 在工程技术和科学研究中，常常遇在工程技术和科学研究中，常常遇到如下情况：到如下情况：1. f(x)的结果复杂，求原函数困难。的结果复杂，求原函数困难。2. f(x)的原函数不存在，或不能用的原函数不存在，或不能用初等函数表示。初等函数表示。3. f(x)无函数式，只有函数表。无函数式，只有函数表。27R1R2F 颗粒尺寸，颗粒尺寸，R 分布密度，分布密度，f 例：计算 R

16、1 R R2 的颗粒数量(分数) 21)(RRdRRfF原始方法：曲线下面面积所占方格数目数值积分雏形 数格子28当当 f (x) 是闭区间上的连续函数时，如何是闭区间上的连续函数时，如何计算下列定积分？计算下列定积分？badxxffI)()(积分中值定理积分中值定理)()()()(fabdxxffIba定积分的计算定积分的计算29上述公式的几何意义上述公式的几何意义底为底为 而高为而高为 的矩形面积恰好等于的矩形面积恰好等于所求曲边梯形的面积所求曲边梯形的面积 。ba( )f( )I f问题问题：如何确定：如何确定 ？I(f)Oxy30)()()()(afabdxxffIba)()()()(

17、bfabdxxffIba)2()()()(bafabdxxffIba( )( )( )( )()2baf af bI ff x dxba左矩形左矩形右矩形右矩形矩形矩形中矩形中矩形定积分常用的近似公式定积分常用的近似公式31 线性近似：梯形公式ab梯形近似xyy = f(x)( )( )( )2babaf x dxf af b简单，但是ab梯形近似xyy = f(x)32 抛物线近似：辛普森(Simpson)公式xyy = f(x)ab2ba 抛物线近似 )(24)(6)(bfbafafabdxxfba通常，抛物线近似优于梯形近似但是，有特例 ab2ba xyy = f(x)33利用梯形公式，

18、利用梯形公式，Simposn公式公式公式计算公式计算 ，并与精确值进行比较。，并与精确值进行比较。dxxI15 . 04267767. 0) 15 . 0(25 . 015 . 0dxxI解：梯形公式解：梯形公式Simposn公式公式10.50.5( 0.5 4 0.75 1) 0.430934036Ixdx43096441. 03215 . 0315 . 0 xdxxI精确值精确值34插值型求积公式插值型求积公式 利用插值多项式来构造求积公式：在积分利用插值多项式来构造求积公式：在积分区间区间a, b上取一组点上取一组点 用用f(x)的的n次插值多项式次插值多项式 来近似代替被积函数来近似代

19、替被积函数f(x) bxxxan 10 nkkkn)x(l )x(f)x(L00( )()( )nbbkkaakf x dxf xlx dx35若记若记 bankjjjkjbakkdxxxxxdx)x(lA0则得数值求积公式则得数值求积公式 nkkkba)x(fAdx)x(f0称为称为插值型求积公式插值型求积公式 36数值积分公式的余项数值积分公式的余项1(1)1101 ( )()( ) ( )( )( )(1)!( , )( )()()()nbnkkaknbbnnaannRff x dxA f xff xL x dxx dxna bxxxxxxx37数值积分公式的代数精度数值积分公式的代数精

20、度nkkkbaxfAdxxf0)()(如果对于任意不高于如果对于任意不高于m次的多项式都准确地成立，次的多项式都准确地成立，而对于而对于 不能准确地成立，则称该求积公式不能准确地成立，则称该求积公式的的代数精度为代数精度为m.1mx数值求积方法是近似方法，为保证精度，希望所提供数值求积方法是近似方法，为保证精度，希望所提供求积公式对于求积公式对于“尽可能多尽可能多”的函数是准确的。的函数是准确的。38几个常用的求积公式的代数精度1.梯形公式的代数精度222( )1( ) ( ( )( )( )2( )11( )()22 ( ( )( )()( )22bbaababbbaaabaf xf x d

21、xdxbabaT ff af bbaf x dxf xxf x dxxdxxbababaT ff af babf x dx当时当时392233322( )11( )()33 ( ( )( )()( )22bbbaaabaf xxf x dxx dxxbababaT ff af babf x dx当时所以梯形公式具有一次的代数精度402.辛普森公式的代数精度2222( )1( )( )4()( )62( )( )1( )()2( )4()( )621(4)()622( )bbaababbaabafxfx dxdxbababaS ff aff bbafx dxS ffxxfx dxxdxbabab

22、aS ff aff bbaababbafx dxS f当时所 以成 立当时所 以成 立4122332222233( )1( )()3 ( ( )4 ()( )62(4()621(222)()63( ) bbaabaf xxf x dxx dxbababaS ff aff bbaababbaaabbbaf x dxS f当时即精确成立423344333332233322344( )1( )()4 ( ( )4 ()( )62(4()621(33)6231()()624( ) bbaabaf xxf x dxx dxbababaS ff aff bbaababbaaaa babbbbaaa bab

23、bbaf x dxS f当时即精确成立可以证明对可以证明对 等式不成立等式不成立辛普森公式具有三次代数精度辛普森公式具有三次代数精度4( )f xx43n +1个结点的插值型数值积分公式至少有个结点的插值型数值积分公式至少有 n 阶代数精度。阶代数精度。(1)( )( )(1)!nbafR fx dxn因此对于次数不大于因此对于次数不大于n的多项式的多项式其余项其余项 0R f 44 求插值型求积公式求插值型求积公式 )21()21()(1011fAfAdxxf 并确定其代数精度。并确定其代数精度。分析：分析：实际上该题目是求实际上该题目是求A0，A1，并确定其，并确定其代数精度代数精度。21

24、,2110 xx因因为为是是插插值值型型的的，且且， 11)21()21()(ffdxxf从而求积公式为从而求积公式为有有2个结点个结点,因而代数精度大于等于因而代数精度大于等于1，从，从m=2开始验证代数精度开始验证代数精度。，从从而而1 m解解 1100)(dxxlA 1111)(dxxlA 11232dxx21)21()21( ff,)(2xxf 对对 11, 1)21(dxx1)21(11 dxx45在等分情况下，称为牛顿在等分情况下，称为牛顿-科茨公式科茨公式将将a, bn等分等分 为步长，节点为等分点为步长，节点为等分点nabh (0,1,2,., )kxakhknx=a+thx

25、xj =(t - j)h4600()()nnkjj ktj hAhdtkj h bankjjjkjbakkdxxxxxdx)x(lA0 x xj =(t - j)hdx =hdt00( )00( 1)()!()!( 1)()()()!()!n knnjj kn knnnkjj khtj dtk nkbatj dtba Cn k nk科茨系数科茨系数47 科特斯系数表48 且满足 注：注：Cotes 系数系数仅取决于仅取决于 n 和和 i，可通过查表得到。与，可通过查表得到。与被积函数被积函数 f (x) 及积分区间及积分区间 a, b 均无关。均无关。q 牛顿牛顿- -科特斯公式科特斯公式ni

26、inibaxfCabxxf0)()()(d)(2) )()(ninniCC(1) 10)(niniC49梯形公式梯形公式 当当n=1, x0=a, x1=b时，有时，有 2201011010abdxabaxdxxxxxAabdxbabxdxxxxxAbabababa2)b(f)a(fabdx)x(fba 50辛普森公式辛普森公式 当当n=2, x0=a, x1=(a+b)/2, x2=b时，有时，有6646120210221012012010210abdx)xx)(xx()xx)(xx(Aabdx)xx)(xx()xx)(xx(Aabdx)xx)(xx()xx)(xx(Abababa)b(f)

27、ba(f)a(fabdx)x(fba24651柯特斯公式柯特斯公式其中其中)4 , 3 , 2 , 1 , 0(4kabkaxk01234()( )7 ( ) 32 ( ) 12 ( ) 32 ( ) 7 ( )90bab af x dxf xf xf xf xf x52 定理定理：若：若 在在a,b上连续，梯形公式的截断上连续，梯形公式的截断误差为误差为)( xf)(12)(31fabfR若若 在在a,b上连续，辛普生公式的截断误差为上连续，辛普生公式的截断误差为)()4(xf)()2(901)4(52fabfR若若 在在a,b上连续，柯特斯公式的截断误差为上连续，柯特斯公式的截断误差为)(

28、) 6 (xf7(6)48()( )9454baRff 53实际计算中，当积分区间较大时，直接使用这些积分公式精度难以保证。为了提高计算精度，采用复合求积的方法。将区间a,b适当分割成若干个子区间，对每个子区间使用求积公式，构成所谓的复化求积公式，这是提高积分精度的一个常用的方法。54xyy = f(x)ab将 a, b 分成若干小区间xk-1, xk, (k = 1, 2, , n)x0 x1x2xnxn-1h前提：节点函数值 f(xk) 已知11( )( )2()( )2bnkkahf x dxf af xf b复合梯形公式复合梯形公式55将积分区间将积分区间 a,b 划分为划分为n等分等

29、分, ,记子区间记子区间 的中点为的中点为 ，在每个小区间上应用辛普森，在每个小区间上应用辛普森公式，则有公式，则有 1,kkxxhxxkk2121)()(4)(6)()(11010211kkkbankxxnkxfxfxfhdxxfdxxfIkk)()(2)(4)(6101121bfxfxfafhnknkkk)()(2)(4)(6101121bfxfxfafhSnknkkkn记记 复合辛普森公式复合辛普森公式5611110042113014( )7 ( )32() 12()9032() 14()7 ( )nnbnakkkknnkkkkhf x dxCf af xf xf xf xf b复合柯特

30、斯公式复合柯特斯公式hxxhxxhxxkkkkkk43;21;41432141其中其中57误差分析误差分析A、复合梯形公式的误差、复合梯形公式的误差)(12)( )( 12)(22fhabafbfhIdxxfnba或B、复合、复合Simpson公式的误差公式的误差4(3)(3)4(4 )1( )() ( )( )1802()()1802bnahfx dxSfbfabahf 或58C、复合、复合Cotes公式的误差公式的误差6(5 )(5 )6( 6 )2()() ( )()94542()()()9454bnahfx dxCfbfabahf 或59 用复化梯形公式计算定积分用复化梯形公式计算定积

31、分 才能使误差不超过才能使误差不超过 10dxeIx51021解解: :取取 , ,则则 , ,又区间长度又区间长度b-a=1b-a=1，对，对复化梯形公式有余项复化梯形公式有余项 xexf)(xexf )(52210211121)(12)( enfhabxRT即即 ,n212.85,n212.85，取，取n=213n=213，即将区间，即将区间0,10,1分为分为213213等份时，用复化梯形公式计算误差等份时，用复化梯形公式计算误差不超过不超过 。 52106en51021问区间问区间0,10,1应分多少等份应分多少等份60p复合求积方法对提高精度是一种行之有效的方法，但是复合求积方法对提

32、高精度是一种行之有效的方法，但是在使用求积公式之前必须先给出步长。在使用求积公式之前必须先给出步长。p步长取得太大精度难以保证，步长太小则会导致计算量步长取得太大精度难以保证，步长太小则会导致计算量的增加，而事先给出一个合适的步长往往是困难的。的增加，而事先给出一个合适的步长往往是困难的。p为了解决这个问题，较好的方法就是采用为了解决这个问题，较好的方法就是采用变步长变步长的计算的计算方法方法p即在步长逐次减半的过程中，反复利用复合求积公式进即在步长逐次减半的过程中，反复利用复合求积公式进行计算，直到所求得的积分值满足精度要求为止。行计算，直到所求得的积分值满足精度要求为止。610123456

33、7891000.511.501234567891000.511.501234567891000.511.562 设将积分区间设将积分区间 a,b n等分，即分成等分，即分成n个子区间，个子区间，一共有一共有n+1个节点，即个节点，即x=a+kh, , k=0,1,，n，步，步长长 。对于某个子区间。对于某个子区间 , ,利用梯形公利用梯形公式计算积分近似值有式计算积分近似值有 nabh1,kkxx)()(21kkxfxfh11011()()2( )2()( )2nnkkknkkhTf xf xhf af xf b对整个区间对整个区间a,ba,b有有63将子区间将子区间 再二等份再二等份, ,取

34、其中点取其中点作新节点作新节点, ,此时区间数增加了一倍为此时区间数增加了一倍为2n,2n,对某个子区对某个子区间间 , ,利用复化梯形公式计算其积分近似值利用复化梯形公式计算其积分近似值 。1,kkxx)(21121kkkxxx1,kkxx)()(2)(4121kkkxfxfxfh对整个区间对整个区间a,ba,b有有 1012)()(2)(421nkkkknxfxfxfhT10101)(2)()(421nkknkkkxfhxfxfh比较比较 和和 有有nTnT2102)(2221nkknnxfhTT64 当把积分区间分成当把积分区间分成n等份，用复化梯形等份，用复化梯形公式计算积分公式计算积

35、分I的近似值的近似值 时，截断误差为时，截断误差为 nT)(122nnnfnababTIR 若把区间再分半为若把区间再分半为2n2n等份，计算出定积分等份，计算出定积分的近似值的近似值 ，则截断误差为，则截断误差为 nT2)(2122222nnnfnababTIR 当当 在区间在区间a,ba,b上变化不大时上变化不大时, ,有有 )(xf )()(2nnff 214nnITIT所以所以 65 可见可见, ,当步长二分后误差将减至当步长二分后误差将减至 , ,将将上式移项整理，可得验后误差估计式上式移项整理，可得验后误差估计式 41)(3122nnnTTTI上式说明，只要二等份前后两个积分值上式

36、说明，只要二等份前后两个积分值和和 相当接近，就可以保证计算结果相当接近，就可以保证计算结果的误差很小，使的误差很小，使 接近于积分值接近于积分值I I。 nTnT2nT2nT266变步长的梯形求积算法实现变步长的梯形求积算法实现（1 1）变步长的梯形求积法的计算步骤）变步长的梯形求积法的计算步骤 变步长梯形求积法。它是以梯形求积公式为变步长梯形求积法。它是以梯形求积公式为基础，逐步减少步长，基础，逐步减少步长，按如下递推公式求二分后的按如下递推公式求二分后的梯形值梯形值102)(2221nkknnxfhTT其中其中Tn和和T2n分别代表二等分前后的积分值分别代表二等分前后的积分值 如果如果

37、, (为给定的误差限为给定的误差限 ) 则则T2n作为积分的近似值作为积分的近似值, 否则继续进行二等分否则继续进行二等分, 即即转转 再计算，直到满足所要求的精度为止，最终取再计算，直到满足所要求的精度为止，最终取二分后的积分值二分后的积分值T2n 作为所求的结果作为所求的结果 nnTT2nnTThh2,267例例 用变步长梯形求积法计算定积分用变步长梯形求积法计算定积分解解: : 先对整个区间先对整个区间 0,10,1 用梯形公式用梯形公式, ,对于对于 10dsinxxxI8410709. 0) 1 (, 1)0(,sin)(ffxxxf所以有所以有 9207355. 0) 1 ()0(

38、211ffT然后将区间二等份然后将区间二等份, ,由于由于 , ,故有故有 12( )0.958851f9397933.0)(21212112fTT进一步二分求积区间进一步二分求积区间, ,并计算新分点上的函数值并计算新分点上的函数值 9088516.0)(,9896158.0)(4341ff68有有 9445135.0)()(4121434124ffTT这样不断二分下去，计算结果如这样不断二分下去，计算结果如P P117117所示。所示。积分的准确值为积分的准确值为0.9460830.946083。 69龙贝格求积公式龙贝格求积公式 变步长梯形求积法算法简单，但精度较差，收敛变步长梯形求积法

39、算法简单，但精度较差，收敛速度较慢，但可以利用梯形法算法简单的优点，形速度较慢，但可以利用梯形法算法简单的优点，形成一个新算法，这就是龙贝格求积公式。龙贝格公成一个新算法，这就是龙贝格求积公式。龙贝格公式又称逐次分半加速法。式又称逐次分半加速法。 根据积分区间分成根据积分区间分成n等份和等份和2n等份时的误差估计等份时的误差估计式可得式可得 )(3122nnnTTTI 积分值积分值 的误差大致等于的误差大致等于 , ,如果用如果用 对对 进行修正时，进行修正时， 与与 之和比之和比 更接近积分真值更接近积分真值nT2)(312nnTT)(312nnTTnT2)(312nnTTnT2nT2702

40、22141()333nnnnnnTTTTTT考察考察 与与n等份辛普森公式等份辛普森公式 之间的关系。之间的关系。将复化梯形公式将复化梯形公式 nS)()(2)(211bfxfafhTnkkn梯形变步长公式梯形变步长公式 102)(2221nkknnxfhTT代入代入 表达式得表达式得 nT121100( )4()2()( )6nnnknkkkhTf af xf xf bSnnnTTS31342故故 nT71nnnnnTTTTTT3134)(31222 (6.9) TnS这就是说，用梯形法二分前后两个积分值这就是说，用梯形法二分前后两个积分值 和和 作线性组合，结果却得到复化辛普森公式计算得到

41、作线性组合，结果却得到复化辛普森公式计算得到的积分值的积分值 。 nTnT2nS72再考察辛普森法。其截断误差与再考察辛普森法。其截断误差与 成正比，因此，成正比，因此，如果将步长折半，则误差减至如果将步长折半，则误差减至 ，即有，即有 4h1/161612nnSISI由此可得由此可得 nnSSI15115162可以验证可以验证, ,上式右端的值其实等于上式右端的值其实等于C Cn n，就是说，用辛，就是说，用辛普森公式二等份前后的两个积分值普森公式二等份前后的两个积分值S Sn n和和S S2n2n 作线性组作线性组合后，可得到柯特斯公式求得的积分值合后，可得到柯特斯公式求得的积分值C Cn

42、 n，即有，即有 nnnSSC1511516273用同样的方法，根据柯特斯公式的误差公式，可进用同样的方法，根据柯特斯公式的误差公式，可进一步导出龙贝格公式一步导出龙贝格公式 nnnCCR63163642p在变步长的过程中运用上述公式，就能将粗糙的梯在变步长的过程中运用上述公式，就能将粗糙的梯形值形值Tn逐步加工成精度较高的辛普森值逐步加工成精度较高的辛普森值Sn、柯特斯值、柯特斯值Cn和龙贝格值和龙贝格值Rnp或者说，将收敛缓慢的梯形值序列或者说，将收敛缓慢的梯形值序列Tn加工成收敛迅加工成收敛迅速的龙贝格值序列速的龙贝格值序列Rn，这种加速方法称为龙贝格算法，这种加速方法称为龙贝格算法74

43、龙贝格求积法计算步骤龙贝格求积法计算步骤 用梯形公式计算积分近似值用梯形公式计算积分近似值 按变步长梯形公式计算积分近似值按变步长梯形公式计算积分近似值 将区间逐次分半将区间逐次分半, ,令区间长度令区间长度 )()(21bfafabT), 2 , 1 , 0(2kabhk102)(2221nkknnxfhTT计算计算)2(kn 按加速公式求加速值按加速公式求加速值 322nnnnTTTS1522nnnnSSSC6322nnnnCCCR梯形加速公式：梯形加速公式： 辛普森加速公式：辛普森加速公式： 龙贝格求积公式：龙贝格求积公式： 75 精度控制；直到相邻两次积分值精度控制；直到相邻两次积分值

44、 nnRR2（其中（其中为允许的误差限）则终止计算并取为允许的误差限）则终止计算并取R Rn n作为积分作为积分 的近似值，否则将区间再对分，重的近似值，否则将区间再对分，重复复 ， 的计算，直到满足精度要求为止。的计算，直到满足精度要求为止。 badxxf)(76例例 用龙贝格算法计算定积分用龙贝格算法计算定积分 要求相邻两次龙贝格值的偏差不超过要求相邻两次龙贝格值的偏差不超过解解: :由题意由题意 102d14xxI510214)(, 1,0 xxfba3)24(21) 1 ()0(211ffT1 . 351621321)(21212112fTT13118. 3)56. 2764. 3(4

45、11 . 321)()(4121434124ffTT13899. 3)()()()(81218785838148ffffTT14094. 3)()()()()()()()(16121161516131611169167165163161816ffffffffTT771333.33134121TTS14157. 33134242TTS14159.33134484TTS14159.331348168TTS14212. 31511516121SSC14159. 31511516242SSC14159. 31511516484SSC7814158. 36316364121CCR14159. 36316

46、364242CCR由于由于 ，于是有，于是有 00001. 012RR14159. 3d14102xxI79本章小结本章小结 本章介绍了积分和微分的数值计算方法，其基本本章介绍了积分和微分的数值计算方法，其基本原理主要是逼近论，即设法构造某个简单函数原理主要是逼近论，即设法构造某个简单函数P(x)P(x)近近似表示似表示f(x)f(x)，然后对，然后对P(x)P(x)求积或求导得到求积或求导得到f(x)f(x)的积分的积分或导数的近似值。基于插值原理，推导了数值积分和或导数的近似值。基于插值原理，推导了数值积分和数值微分的基本公式。数值微分的基本公式。80MATLAB数值积分数值积分基于变步长

47、辛普生法，基于变步长辛普生法，MATLAB给出了给出了quad函数来求定积分。函数来求定积分。该函数的调用格式为：该函数的调用格式为： I,n=quad(fname,a,b,tol,trace)fname是被积函数名。是被积函数名。a和和b分别是定积分的下限和上限。分别是定积分的下限和上限。tol用来控制积分精度，缺省时取用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非控制是否展现积分过程，若取非0则展现积分过程，取则展现积分过程，取0则则不展现，缺省时取不展现，缺省时取trace=0。返回参数返回参数I即定积分值，即定积分值，n为被积函数的调用次数。为被积函数的调用次数。81MATLAB数值积分实例数值积分实例 (1) 建立被积函数文件建立被积函数文件fesin.m。function f=fesin(x)f=exp(-0.5*x).*sin(x+pi/6); (2) 调用数值积分函数调用数值积分函数quad求定积分。求定积分。S,n=quad(fesin,0,3*pi)S = 0.9008n = 7782quad8函数函数基于牛顿柯特斯法，