数值第4章数值积分与数值微分

1、数理学院数理学院SCHOOL OF MATHEMATICS AND PHYSICS4 .1 引言引言4 .2 牛顿牛顿-柯特斯公式柯特斯公式4 .3 复化求积公式复化求积公式4 .4 龙贝格求积公式龙贝格求积公式4 .5 高斯求积公式高斯求积公式4 .6 数值微分数值微分第四章第四章 数值积分与数值微分数值积分与数值微分一、数值求积的基本思想一、数值求积的基本思想)()()(aFbFdxxfba 积分积分 只要找到被积函数只要找到被积函数 f (x)原函数原函数F(x)，便有，便有牛顿牛顿莱布尼兹莱布尼兹(NewtonLeibniz)公式公式 baxxfId)(实际困难实际困难：大量的被积函数

2、（：大量的被积函数（ , sin x2 等）等）, 找不到用初等函找不到用初等函数表示的原函数数表示的原函数；另外；另外, 当当f (x)是（测量或数值计算出的）一张是（测量或数值计算出的）一张数据表时，数据表时，牛顿牛顿莱布尼兹公式莱布尼兹公式也也不能直接运用不能直接运用。xxsin 积分中值定理：在积分中值定理：在a, b内存在一点内存在一点 ，有，有 f( )成立。成立。 )(d)(abxxfba 4.1 引言引言 就是说就是说, 底为底为b- -a 而高为而高为f( )的的矩形面积矩形面积恰恰等于所求等于所求曲边梯形的面积曲边梯形的面积 .问题问题 在于点在于点的具体位置一般是不知道的

3、，因而的具体位置一般是不知道的，因而 难以准确算出难以准确算出 f( )的值的值我们将我们将f ( )称为区间称为区间a, b上的平均高度这样上的平均高度这样,只要对只要对平均高度平均高度f( )提供一种算法提供一种算法,相应地便获得一种数值求积方法相应地便获得一种数值求积方法 如果用两端点的如果用两端点的“高度高度”f(a)与与f(b)的算术平均作为平均高度的算术平均作为平均高度f ( ) 的近似值，这样导出的求积公式的近似值，这样导出的求积公式 : 便是我们所熟悉的便是我们所熟悉的梯形公式梯形公式 . )()(2bfafabT 2)(bafabR2bac 而如果改用区间中点而如果改用区间中

4、点 的的“高度高度”f (c)近似地取代平近似地取代平均高度均高度f ( )，则又可导出所谓，则又可导出所谓中矩形公式中矩形公式(今后简称矩形公式今后简称矩形公式)：(1.1)(1.2) 更一般地，我们可以在区间更一般地，我们可以在区间a，b上适当选取某些节点上适当选取某些节点 xk ，然后然后用用 f (xk )加权平均得到平均高度加权平均得到平均高度 f ()的近似值的近似值，这样构造出的，这样构造出的求积公式具有下列形式求积公式具有下列形式式中式中 xk 称为称为求积节点求积节点；Ak 称为称为求积系数求积系数，亦称为伴随节点，亦称为伴随节点 xk 的的权权权权Ak 仅仅与节点仅仅与节点

5、xk 的选取有关，而不依赖于被积函数的选取有关，而不依赖于被积函数 f(x)的的具体形式具体形式 ban0kkkxfAdxxf)()(使积分公式具有通用性使积分公式具有通用性 这类数值积分方法通常称作能这类数值积分方法通常称作能机械求积机械求积， 其特点是将积分其特点是将积分求值问题归结为函数值的计算，这就避开了牛顿求值问题归结为函数值的计算，这就避开了牛顿莱布尼兹公莱布尼兹公式需要寻求原函数的困难式需要寻求原函数的困难(1.3)二、二、代数精度的概念代数精度的概念 数值求积方法是近似方法，为要保证精度，我们自然希望求积数值求积方法是近似方法，为要保证精度，我们自然希望求积公式能对公式能对“尽

6、可能多尽可能多”的函数准确地成立，这就提出了所谓代数精的函数准确地成立，这就提出了所谓代数精度的概念度的概念 定义定义 1 如果某个求积公式对于次数如果某个求积公式对于次数m的多项式均能准确地成的多项式均能准确地成立，但对于立，但对于m+1次多项式就不准确成立，则称该求积公式具有次多项式就不准确成立，则称该求积公式具有m次代次代数精度数精度(或代数精确度或代数精确度) 一般地，欲使求积公式一般地，欲使求积公式 具有具有m次代数次代数精度，只要令它对于精度，只要令它对于f (x) = 1，x，xm 都能准确成立，这就要求都能准确成立，这就要求 bankkkxfAxxf0)(d)( . )(11;

7、)(21;1122mmmkkkkkabmxAabxAabA例例1： 考察其代数精度。考察其代数精度。 f(x)abf(a)f(b)梯形公式梯形公式解：解：逐次检查公式是否精确成立逐次检查公式是否精确成立代入代入 f(x) = 1： baabdx111 2 ab=代入代入 f(x) = x ：=代入代入 f(x) = x2 ： 222abbadxx 2baab 3233abbadxx 222baab 代数精度代数精度 = 1)()(2)(bfafabdxxfba 例例2 试构造形如试构造形如 f(x)dx A0f(0)+ A1f(h)+ A2f(2h) 的数值的数值求积公式求积公式,使其代数精度

8、尽可能高使其代数精度尽可能高,并指出其代数精度的阶数并指出其代数精度的阶数.3h0解解: 令公式对令公式对 f(x)=1,x, x2 均准确成立均准确成立,则有则有3h=A0+ A1+ A2h2=0 + A1h+ A22h9h3=0 + A1h2+ A24h229故求积公式的形式为故求积公式的形式为解之得解之得 A0= h, A1=0, A2= h. 34 34 f(x)dx f(0) + f(2h)3h43h43h0由公式的构造知由公式的构造知,公式公式至少至少具有具有2次代数精度次代数精度; 而当而当f(x)=x3时时,公式的左边公式的左边= h4, 右边右边=18h4, 公式的左边公式的

9、左边 右边右边,说明说明此公式对此公式对 f(x)=x3不能准确成立不能准确成立.因此因此,公式只具有公式只具有2次代数次代数精度精度.814三、三、求积公式的收敛性与稳定性求积公式的收敛性与稳定性 定理定理3表明，只要求积系数表明，只要求积系数Ak0 (k0,1,n)，就能保证，就能保证计算的稳定性计算的稳定性 定义定义2 在求积公式在求积公式 中，若中，若 其中其中 ，则称求积公式是收敛的，则称求积公式是收敛的 由于计算由于计算 f (xk)可能有误差可能有误差,实际得到实际得到 定义定义3 对任给对任给 e e 0，若，若 (k=0,1, ,n), 就有就有 , 则称求积公式是稳定的则称

10、求积公式是稳定的. bankkkxfAxxf0)(d)(e e |)(|00nkkknkkkfAxfA)(max11 iinixxhkkfxf)(0 ，只只要要 .)(,kkkkfxff 即即 bankkkhnxxfxfAd)()(lim00 定理定理3 若求积公式若求积公式(13)中系数中系数Ak0 (k0，1，n)，则此求积公式是稳定的则此求积公式是稳定的近似近似计算计算 badxxfI)(思思路路利用利用插值多项式插值多项式 则积分易算。则积分易算。)()(xfxPn 在在a, b上取上取 a x0 x1 x=1,2； y=exp(1./x) I=1/2*(y(1)+y(2)*1I =

11、2.1835 R1=(2-1)3/12*8.1548R1 = 0.6796估计估计截断误差截断误差为为用用Simpson公式公式计算：计算：=2. 0263=198.4306890. 0)(max2880) 12()4(2152xfRx=0.06890 x=1:0.5:2; y=exp(1./x)y = 2.7183 1.9477 1.6487 I2=(2-1)/6*(y(1)+4*y(2)+y(3)I2 = 2.0263 f41=exp(1/x1)*(1./x18+12./x17+36./x16+24./x15)f41 = 198.4346 R2=(2-1)5/2880*f41R2 = 0.

12、06894.3 4.3 复合求积公式复合求积公式高次插值有高次插值有Runge 现象现象，故采用分段低次插值，故采用分段低次插值 分段低次合成的分段低次合成的 Newton-Cotes 复合复合求积公式。求积公式。一、复合梯形公式一、复合梯形公式:),., 0(,nkhkaxnabhk 在每个在每个 上用梯形公式：上用梯形公式：,1 kkxx 11)()(2)(2nkkbfxfafh bankkkxfxfhdxxf101)()(2)(=Tn),(),()(12)()(12)(122102103bafabhnfabhfhfRnkknkk /*中值定理中值定理*/1,., 0,)()(2)(111

13、 nkxfxfxxdxxfkkxxkkkk二、复合辛普森公式二、复合辛普森公式:),., 0(,nkhkaxnabhk )()(4)(6)(1211 kkkxxxfxfxfhdxxfkkkx21 kx1 kx44444 )()(2)(4)(6)(1010121 nknkkkbabfxfxfafhdxxf= Sn)(2180)4(4 fhabfR 注：注：为方便编程，可采用另一记法：令为方便编程，可采用另一记法：令 n = 2n 为偶数，为偶数， 这时这时 ，有，有hkaxhnabhk ,2 )()(2)(4)(3 koddkevenkknbfxfxfafhS三、收敛速度与误差估计：三、收敛速度

14、与误差估计：定义定义 若一个积分公式的误差满足若一个积分公式的误差满足 且且C 0，则则称该公式是称该公式是 p 阶收敛阶收敛的。的。 ChfRphlim0)(,)(,)(642hOChOShOTnnn例例4：计算计算dxx 10142 解：解： )1()(2)0(161718fxffTkk8kxk 其中其中= 3.138988494 )1()(2)(4)0(241oddeven4fxfxffSkk8kxk 其中其中= 3.141592502运算量基运算量基本相同本相同function t=rctrap(fun,a,b,n)%复化梯形公式%n等分%a,b区间的左右端点h=(b-a)/n;t=0

15、;for i=1:n-1 t=t+h*feval(fun,i*h+a);endt=t+0.5*h*(feval(fun,a+eps)+feval(fun,b)function f=fun4(x)f=4/(1+x2) t8=rctrap(fun4,0,1,8)function s=sptrap(fun,a,b,n)% n，对应的等分点为2nh=(b-a)/(2*n);s=0;for i=1:ns=s+feval(fun,(2*i-1)*h+a)*4;endfor i=1:n-1s=s+feval(fun,(2*i)*h+a)*2;ends=s+feval(fun,a+eps)+feval(fun

16、,b);s=s*h/3 s4=sptrap(fun4,0,1,2)问题问题: 给定精度给定精度 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 409 k = 9

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

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

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

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

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

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

23、重复同样的手续，用柯特斯法二分前后的两个积分值Cn与与C2n的线性组合的结果可得到的线性组合的结果可得到龙贝格龙贝格(Romberg)求积公式求积公式nnnnnCCCCR631636414114423233 我们在变步长的过程中运用加速公式，就能将粗糙的梯我们在变步长的过程中运用加速公式，就能将粗糙的梯形值形值Tn逐步加工成精度较高的辛普森值逐步加工成精度较高的辛普森值Sn 、柯特斯值、柯特斯值Cn和龙和龙贝格值贝格值Rn . Romberg 算法：算法： e e ? e e ? e e ? T1 =)0(0T T8 =)3(0T T4 =)2(0T T2 =)1(0T S1 =)0(1T R

24、1 =)0(3T S2 =)1(1T C1 =)0(2T C2 =)1(2T S4 =)2(1T一般有：一般有：nnnSTT 1442nnnCSS 144222nnnRCC 144323Romberg 序列序列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 区间等分数区间等分数 梯形序列梯形序列 辛普森序列辛普森序列 柯特斯序列柯特斯序列 龙贝格序列龙贝格序列 龙贝格求积算法可用下表来表示：龙贝格求积算法可用下

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

26、 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.356

27、 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。三、理查森三、理查森(Richardson)外

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

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

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

31、若若取取算算法法上上式式也也称称为为龙龙贝贝格格求求积积Tm Romberg 算法算法 可以证明，如果可以证明，如果 f (x) 充分光滑，那么充分光滑，那么T 数表每一列的数表每一列的元素及对角线元素均收敛到所求的积分值元素及对角线元素均收敛到所求的积分值 I ，即，即ITmITkmmkmk )()(lim)(lim，固固定定例用龙贝格算法计算积分例用龙贝格算法计算积分.102/3dxxI function s,n=rbg1(fun,a,b,eps) if nargineps) h=(b-a)/2(k-1); w=0; if(h=0) for i=1:(2(k-1)-1) w=w+f (a+

32、i*h); end t(k,1)=h*(f(a)/2+w+f(b)/2) for l=2:kfor i=1:(k-l+1) t(i,l)=(4(l-1)*t(i+1,l-1)-t(i,l-1)/(4(l-1)-1)endends=t(1,k); s0=(t(1,k-1); k=k+1; n=k; else s=s0; n=-k; end end function f=fun(x)f=x(3/2) R=rbg1(fun,0,1,1e-6)自适应积分方法自适应积分方法(补充内容补充内容) 复合求积方法通常适用于被积函数变化不太大的积分复合求积方法通常适用于被积函数变化不太大的积分,如果在求积区间中

33、被积函数变化很大如果在求积区间中被积函数变化很大,有的部分函数值变化有的部分函数值变化剧烈剧烈,另一部分变化平缓另一部分变化平缓.这时统一将区间等分用复合求积公这时统一将区间等分用复合求积公式计算积分工作量大式计算积分工作量大. 要达到误差要求对变化剧烈部分必须将区间细分要达到误差要求对变化剧烈部分必须将区间细分,而平而平缓部分则可用大步长缓部分则可用大步长.针对被积函数在区间不同情形采用不针对被积函数在区间不同情形采用不用的步长用的步长,使得在满足精度前提下积分计算工作量尽可能小使得在满足精度前提下积分计算工作量尽可能小. 在不同区间上预测被积函数变化的剧烈程度确定相应在不同区间上预测被积函

34、数变化的剧烈程度确定相应步长步长,这种方法称为自适应积分方法这种方法称为自适应积分方法设给定精度要求设给定精度要求, 0 e e计算积分计算积分 badxxffI)()(的近似值的近似值.先取步长先取步长h=b-a,应用辛普森公式有应用辛普森公式有),(),()2(180),()()()4(4bafhabbaSdxxffIba (5.1)其中其中).()2(4)(6),(bfbafafhbaS 若把区间若把区间a,b对分对分,步长步长,222abhh 在每个小区间上用辛普森在每个小区间上用辛普森公式公式,则得则得),(),()2(180),()()4(422bafhabbaSfI (5.2)其

35、中其中).()43(4)2(6),2(),2()4(4)(6)2,(),2()2,(),(222bfhafhafhbbaShafhafafhbaaSbbaSbaaSbaS 实际上实际上(5.2)式即为式即为),(),()4(180),()()4(42bafhabbaSfI 与与(5.1)式比较，若式比较，若)()4(xf在在(a,b)上变化不大上变化不大,可假定可假定)()()4()4( ff 从而可得从而可得).()2(180),(),(1516)4(42 fhabbaSbaS 与与(5.2)比较比较,则得则得|,|151| ),(),(|151| ),()(|2122SSbaSbaSbaS

36、fI 如如果果有有).,(),(221baSSbaSS 这里这里,15|21e e SS则可得到则可得到,| ),()(|2e e baSfI此时此时,可取可取S2(a,b)作为作为I(f)的近似的近似,则可达到给定的误差精度则可达到给定的误差精度.(5.3)若不等式若不等式(5.3)不成立不成立,则应分别对子区间则应分别对子区间,22,bbabaa 及及再用辛普森公式再用辛普森公式,此时步长此时步长h3=0.5h2,得到得到)2()2,(33bbaSbaaS，及及 只要分别考察只要分别考察2| )2,()(|3e e baaSfI及及2| ),2()(|3e e bbaSfI是否成立是否成立

37、.对满足要求的区间不再细分对满足要求的区间不再细分,对不满足要求的还要继续上述过程对不满足要求的还要继续上述过程.最后还要应用龙贝格法则求出相应区间的积分近似值最后还要应用龙贝格法则求出相应区间的积分近似值.例例7 计算积分计算积分,112 . 02dxx 若用复合辛普森法若用复合辛普森法(3.5)式计算结果见表式计算结果见表此处此处hn即为公式中的即为公式中的h,积分精确值为积分精确值为4).nhnSn|Sn-Sn-1|10.84.9481480.7611120.44.1870370.16281930.24.0242180.02205440.14.0021640.00201050.054.0

38、00154计算到计算到|Sn-Sn-1|0.02为止为止,此时此时I(f)的近似值的近似值S50.2,1=4.000154,若再用龙贝格法则得到若再用龙贝格法则得到00002. 4151 , 2 . 0455 SSSRS整个计算将区间整个计算将区间32等分等分,计算计算33个个f(x)的值的值.现在若用自适应积分法现在若用自适应积分法,当当h2=0.4时有时有S20.2,0.6=3.51851852,S2 0.6,1=0.66851852,由于由于S2= S20.2,1= S20.2,0.6+ S20.6,1=4.187037,|S1-S2|=0.761111大于允许误差大于允许误差0.02,

39、故对故对0.2,0.6及及 0.6,1两区间再用两区间再用h3=h2/2做积分做积分.先计算先计算0.6,1的积分的积分S3 0.6,0.8=0.41678477, S3 0.8,1=0.25002572.由于由于S20.6,1-( S30.6,0.8+ S30.8,1)=0.66851852-0.66681049=0.001708 小于允许误差小于允许误差0.01,故在故在0.6,1区间的积分值为区间的积分值为66669662. 0)66851852. 066681049. 0 0 1 , 6 . 0 RS再计算再计算0.2,0.6的积分的积分, S2 0.2,0.

40、6=3.51851852,而对而对h3=h2/2得得S3 0.2,0.4=2.52314815, S3 0.4,0.6=0.83425926.由于由于S20.2,0.6-( S30.2,0.4+ S30.4,0.6)=0. 161111大于允许大于允许误差误差0.01,因此还要分别计算因此还要分别计算0.2,0.4及及0.4,0.6 积分积分.当当h4=h3/2时可求得时可求得S4 0.4,0.5=0.50005144, S4 0.5,0.6=0.33334864.而而S30.4,0.6-( S40.4,0.5+ S40.5,0.6)=0.000859 小于允许小于允许误差误差0.01,故在故

41、在0.4,0.6区间的积分值为区间的积分值为.8333428. 06 . 0 , 4 . 0 RS而对而对0.2,0.4的积分的积分, S3 0.2,0.4- S4 0.2,0.4不小于不小于0.005,故故还要分别计算还要分别计算 0.2,0.3及及0.3,0.4的积分的积分,其中其中 S40.3,0.4=0.83356954,当当h5=h4/2可求得可求得S5 0.3,0.35=0.47620166, S5 0.35,0.4=0.35714758,且且S40.3,0.4-(S5 0.3,0.35+ S5 0.35,0.4)=0.000220小于允小于允许误差许误差0.0025,故有故有.8

42、3333492. 04 . 0 , 3 . 0 RS最后子区间最后子区间0.2,0.3的积分可检验出它的误差小于的积分可检验出它的误差小于0.0025,且可得且可得.666686. 13 . 0 , 2 . 0 RS将以上各区间的积分近似值相加可得到将以上各区间的积分近似值相加可得到.00005957. 41 , 6 . 06 . 0 , 4 . 04 . 0 , 3 . 03 . 0 , 2 . 0)( RSRSRSRSfI它一共只需要计算它一共只需要计算17个个f(x)的值的值. 在构造在构造Newton-Cotes公式公式时，限定用积分区间的时，限定用积分区间的等分点等分点作为求积节点作

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

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

45、般性,我们研究带权积分我们研究带权积分 (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就很难求解故一般不通过解方程就很难求解故一般不通过解方程

46、(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不是线性方程组，不是线性方程组，不易求解。不易求解。 而从研究而从研究高斯点的基本特性高斯点的基本特性来着手解决来着手解决Gauss 求积公式求积公式的构造问题

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

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

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

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

51、大有性质：任意次数不大于于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 四、常用的高斯型求积公式四、常用的高斯型求积公式Step 2：

52、求求 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 注：注：构造正交多项式也可以利用最小二乘数据拟合中介构造正交多项式也可以

53、利用最小二乘数据拟合中介绍过的递推式进行。绍过的递推式进行。 特殊正交多项式族：特殊正交多项式族： Legendre 多项式族：多项式族：1)( x 定义在定义在 1, 1上，上，以以 Pn+1 的根为节点的求积公式称为的根为节点的求积公式称为Gauss-Legendre 公式公式。注：注：一般一般a,b上的积分可化为上的积分可化为-1,1上特殊高斯公式进行计算。上特殊高斯公式进行计算。例例 用用3点点Gauss公式计算积分公式计算积分 解解 查表得查表得x1=-0.7745966692,x2=0,x3=0.7745966692, A1=A3=0.5555555556,A2=0.8888888

54、889, 所以有所以有 Gauss-Legendre求积公式的余项为求积公式的余项为 ) 1 , 1(,)() 12()!2() !(2)2(3412nnfnnnfR11.cos xdxI68300355. 1coscoscos332211xAxAxAI误差为 5347103492. 6)cos(7) ! 6(62fR实际上实际上,I=2sin1=1.68294197, 误差为误差为|R |=6.158 10-5 . 用用Simpson公式公式,则有则有I 1.69353487, 误差为误差为|R |=1.06 10-2 . 由于由于因此因此,a,b上权函数上权函数 (x)=1的的Gauss型

55、求积公式为型求积公式为batabbaxdttabbafabdxxf)2)()()22(2)(11baniiixabbafAabdxxf1)22(2)(求积误差可表示为求积误差可表示为),(,)() 12()!2() !()()2(3412bafnnnabfRnn例例 用用3点点Gauss公式计算积分公式计算积分结果远比结果远比Simpson公式的结果精确公式的结果精确.102.14dxxI 解解 这里这里Gauss点和积分系数与上例相同点和积分系数与上例相同,所以所以 312112102141068. 32/ )1(1421)2121(142114iiixAdttdxxI Chebyshev

56、多项式族：多项式族：211)(xx 定义在定义在 1, 1上，上，) arccos( cos)(xkxTk Tn+1 的根为的根为 2212cosnkxkk = 0, , n以此为节点构造公式以此为节点构造公式 1102)()(11nkkkxfAdxxfx称为称为 Gauss-Chebyshev 公式公式。注意到积分端点注意到积分端点 1 可能是积分可能是积分的的奇点奇点，用普通，用普通Newton-Cotes公公式在端点会出问题。而式在端点会出问题。而Gauss公公式可能避免此问题的发生。式可能避免此问题的发生。Gauss-Laguerre求积公式为求积公式为 求积公式的误差为求积公式的误差

57、为 由于由于 所以所以,对对0, + )上权函数上权函数 (x)=1的积分的积分,也可以构造类似的也可以构造类似的Gauss-Laguerre求积公式求积公式:01)()(niiixxfAdxxfe), 0(,)()!2() !()2(2nfnnfR00)()(dxxfeedxxfxx01)()(niixixfeAdxxfi4.6 4.6 数值微分数值微分 数值微分的数值微分的概念概念 数值微分的数值微分的计算方法计算方法 原始概念近似原始概念近似: :中点法及外推法中点法及外推法 函数近似函数近似: :插值型的求导公式插值型的求导公式 函数相互关系转化函数相互关系转化: :利用数值积分求导利

58、用数值积分求导 数值微分的数值微分的误差分析误差分析 泰勒展式估计泰勒展式估计 事后误差估计事后误差估计 基本关系转化基本关系转化 数值微分数值微分就是就是用函数值的线性组合近似函数在某点用函数值的线性组合近似函数在某点的导数值的导数值一、中点法和外推法一、中点法和外推法 按导数定义按导数定义 , 是差商是差商 当当 时的极限时的极限取取差商差商作为作为导数导数的近似值的近似值，建立简单的数值微分方法，建立简单的数值微分方法:)(0 xf hxfhxf)()(00 0h hxfhxfxf000 (8.1)向后差商近似导数向后差商近似导数(8.2)(8.3)中心差商近似导数中心差商近似导数 hh

59、xfxfxf 000 hhxfhxfxf2000 hhxfhxfhD200 容易看出，就精度而言，以（容易看出，就精度而言，以（8.3）式更为可取，称）式更为可取，称(8.4)为为 的的中点公式中点公式, 其中其中h为一增量，称为为一增量，称为步长步长 这种数值这种数值微分方法称为微分方法称为中点方法中点方法, 它是前两种方法的算术平均它是前两种方法的算术平均)(0 xf 分别将分别将在在 x=a 处做处做Taylor展开有展开有)(haf )(! 5)(! 4)(! 3)(! 2)()()() 5(5) 4(432afhafhafhafhafhafhaf代入代入D(h)得得 )(! 5)(!

60、 3)()()5(42afhafhafhD,6)()(2MhhDaf 其中其中)(maxxfMhax 现在来考虑中点公式现在来考虑中点公式 的截断误差，的截断误差，hhafhafhD2)()()( (8.5)所以所以截断误差截断误差 )(! 5)(! 3)()() 5(42afhafhafhD(8.6)从截断误差的角度来看，步长从截断误差的角度来看，步长h越小，计算结果越准确。且越小，计算结果越准确。且 所以所以, 在在实际计算时实际计算时，通常，通常采用采用二分步长二分步长及及误差事后估误差事后估计法计法, 在变步长的过程中实现步长的自动选择，在保证截断在变步长的过程中实现步长的自动选择，在