第三章数值积分课件

一、数值求积的基本思想一、数值求积的基本思想 积分积分 只要找到被积函数只要找到被积函数 f(x)原函数原函数F(x)，便有，便有牛顿牛顿莱布尼兹莱布尼兹(NewtonLeibniz)公式公式实际困难实际困难：大量的被积函数（：大量的被积函数（,sin x2 等）等）,找不到用初等函找不到用初等函数表示的原函数数表示的原函数；另外；另外,f(x)是（测量或数值计算出的）一张数是（测量或数值计算出的）一张数据表时，据表时，牛顿牛顿莱布尼兹公式莱布尼兹公式也也不能直接运用不能直接运用。积分中值定理：在积分中值定理：在a,b内存在一点内存在一点 ，有，有 f()成立。成立。1 引言引言第三章第三章 数值积分数值积分(Numerical Integration)就是说就是说,底为底为b-a 而高为而高为f()的的矩形面积矩形面积恰等于所求恰等于所求曲边梯形的面积曲边梯形的面积.问题问题 在于点在于点的具体位置一般是不知道的，因而的具体位置一般是不知道的，因而 难以准确算出难以准确算出 f()的值的值我们将我们将f()称为区间称为区间a,b上的平均高度这样上的平均高度这样,只只要对要对平均高度平均高度f()提供一种算法提供一种算法,相应地便获得一种数值求积方法相应地便获得一种数值求积方法 如果用两端点的如果用两端点的“高度高度”f(a)与与f(b)的算术平均作为平均高度的算术平均作为平均高度f()的近似值，这样导出的求积公式的近似值，这样导出的求积公式:便是我们所熟悉的便是我们所熟悉的梯形公式梯形公式(trapezoidal rule).而如果改用区间中点而如果改用区间中点 的的“高度高度”f(c)近似地取代近似地取代平均高度平均高度f()，则又可导出所谓，则又可导出所谓中矩形公式中矩形公式(今后简称矩形公式今后简称矩形公式)：(1.1)(1.2)更一般地，我们可以在区间更一般地，我们可以在区间a，b上适当选取某些节点上适当选取某些节点 xk，然后然后用用 f(xk)加权平均得到平均高度加权平均得到平均高度 f()的近似值的近似值，这样构造出，这样构造出的求积公式具有下列形式的求积公式具有下列形式式中式中 xk 称为称为求积节点求积节点；Ak 称为称为求积系数求积系数，亦称为伴随节点，亦称为伴随节点 xk 的的权权权权Ak 仅仅与节点仅仅与节点xk 的选取有关，而不依赖于被积函数的选取有关，而不依赖于被积函数 f(x)的的具体形式具体形式使积分公式具有通用性使积分公式具有通用性 这类数值积分方法通常称作这类数值积分方法通常称作机械求积机械求积，其特点是将积分求其特点是将积分求值问题归结为函数值的计算，这就避开了牛顿值问题归结为函数值的计算，这就避开了牛顿莱布尼兹公式莱布尼兹公式需要寻求原函数的困难需要寻求原函数的困难(1.3)二、二、代数精度的概念代数精度的概念 数值求积方法是近似方法，为要保证精度，我们自然希望求积数值求积方法是近似方法，为要保证精度，我们自然希望求积公式能对公式能对“尽可能多尽可能多”的函数准确地成立，这就提出了所谓代数精的函数准确地成立，这就提出了所谓代数精度的概念度的概念 定义定义 1 如果某个求积公式对于次数如果某个求积公式对于次数m的多项式均能准确地成的多项式均能准确地成立，但对于立，但对于m+1次多项式就不一定准确，则称该求积公式具有次多项式就不一定准确，则称该求积公式具有m次代次代数精度数精度 一般地，欲使求积公式一般地，欲使求积公式 具有具有m次代数精次代数精度，只要令它对于度，只要令它对于f(x)=1，x，xm 都能准确成立，这就要求都能准确成立，这就要求例例1：考察其代数精度。考察其代数精度。f(x)abf(a)f(b)梯形公式梯形公式解：解：逐次检查公式是否精确成立逐次检查公式是否精确成立代入代入 P0=1：=代入代入 P1=x：=代入代入 P2=x2：代数精度代数精度=1例例2 试构造形如试构造形如 f(x)dx A0f(0)+A1f(h)+A2f(2h)的数值的数值求积公式求积公式,使其代数精度尽可能高使其代数精度尽可能高,并指出其代数精度的阶并指出其代数精度的阶数数.3h0解解:令公式对令公式对 f(x)=1,x,x2 均准确成立均准确成立,则有则有3h=A0+A1+A2h2=0+A1h+A22h9h3=0+A1h2+A24h229故求积公式的形式为故求积公式的形式为解之得解之得 A0=h,A1=0,A2=h.94 34 f(x)dx f(0)+f(2h)3h49h43h0由公式的构造知由公式的构造知,公式公式至少至少具有具有2次代数精度次代数精度;而当而当f(x)=x3时时,公式的左边公式的左边=h4,右边右边=18h4,公式的左边公式的左边 右边右边,说说明此公式对明此公式对 f(x)=x3不能准确成立不能准确成立.因此因此,公式只具有公式只具有2次代次代数精度数精度.814三、三、求积公式的收敛性与稳定性求积公式的收敛性与稳定性 定理定理3表明，只要求积系数表明，只要求积系数Ak0(k0,1,n)，就能保证就能保证计算的稳定性计算的稳定性 定义定义2 在求积公式在求积公式 中，若中，若 其中其中 ，则称求积公式是收敛的，则称求积公式是收敛的 由于计算由于计算 f(xk)可能有误差可能有误差,实际得到实际得到 定义定义3 对任给对任给 0，若若 (k=0,1,n),就有就有 ,则称求积则称求积公式是稳定的公式是稳定的.定理定理3 若求积公式若求积公式(13)中系数中系数Ak0(k0，1，n)，则此求积公式是稳定的则此求积公式是稳定的近似近似计算计算思思路路利用利用插值多项式插值多项式 则积分易算。则积分易算。在在a,b上取上取 a x0 x1 xn b，做做 f 的的 n 次插值次插值多项式多项式 ，即得到，即得到Ak由由 决定，决定，与与 无关。无关。节点节点 f(x)插值型积分公式插值型积分公式Interpolatory quadrature误差误差2 插值型的求积公式插值型的求积公式与与Newton-Cotes 公式公式关键是关键是f(x)如果求积公式是插值型的如果求积公式是插值型的,按余项式按余项式,对于次数对于次数 n的多项式的多项式 f(x)，其余项其余项R f 等于等于0，因而这时求积公式至少具有，因而这时求积公式至少具有n次代数精度次代数精度定理定理1：形如形如 的求积公式至少有的求积公式至少有 n 次代数精度次代数精度 该该公式为公式为插值型插值型（即：（即：）为便于计算，一般取为便于计算，一般取等距离节点等距离节点得到近似公式得到近似公式:一、一、Newton-Cotes 公式公式2、把把a,b二等分，作二等分，作2次插值，有次插值，有此公式称为此公式称为辛普森（辛普森（Simpson）公式）公式。1、对于对于a,b上上1次插值，有次插值，有此即此即梯形公式梯形公式。节点节点等距分布等距分布：令令Cotes系数系数注：注：Cotes 系数仅取决于系数仅取决于 n 和和 i，可查表得到。与可查表得到。与 f(x)及及区间区间a,b均无关。均无关。3、把把a,b n 等分，用插值等分，用插值Ln(x)近似近似 f(x)积分，有积分，有当当n=4时时,牛顿牛顿-柯特斯公式特别称作柯特斯公式特别称作柯特斯公式柯特斯公式,其形式为其形式为n=1:Trapezoidal Rule/*令令 x=a+th,h=b a,用用中值定理中值定理*/代数精度代数精度=1n=2:Simpsons Rule代数精度代数精度=3n=3:Simpsons 3/8-Rule,代数精度代数精度=3,n=4:Cotes Rule,代数精度代数精度=5,n 为为偶数阶偶数阶的的Newton-Cotes 公式至少有公式至少有 n+1 次代数精度。次代数精度。二、几种低阶求积公式的余项二、几种低阶求积公式的余项三、三、偶阶求积公式的代数精度偶阶求积公式的代数精度 作为插值型的求积公式，作为插值型的求积公式，n 阶的牛顿阶的牛顿-柯特斯公式至柯特斯公式至少具有少具有n 次的插值精度（定理次的插值精度（定理1）。实际的代数精度还可）。实际的代数精度还可进一步提高，一般地，可以证明下述定理进一步提高，一般地，可以证明下述定理:定理定理 2 当阶当阶 n 为偶数时，牛顿为偶数时，牛顿-柯特斯公式柯特斯公式至少有至少有 n+1 次代数精度次代数精度.注：注：由公式知，当由公式知，当n8时，柯特斯系数出现负值，这时，时，柯特斯系数出现负值，这时，初始数据误差将会引起计算结果误差增大，即计算不稳初始数据误差将会引起计算结果误差增大，即计算不稳定。因此，实际计算不用定。因此，实际计算不用n8的牛顿的牛顿-柯特斯公式柯特斯公式.估计估计截断误差截断误差为为解解 用用梯形公式梯形公式计算计算:=2.1835估计估计截断误差截断误差为为=0.6796用用Simpson公式公式计算：计算：=2.0263例例3 试分别使用梯形公式和试分别使用梯形公式和Simpson公式计算积分公式计算积分 的近的近似值，并估计截断误差似值，并估计截断误差.=198.43=0.068903 3 复化求积公式复化求积公式高次插值有高次插值有Runge 现象现象，故采用分段低次插值，故采用分段低次插值 分段低次合成的分段低次合成的 Newton-Cotes 复合复合求积公式。求积公式。一、复化梯形公式一、复化梯形公式:在每个在每个 上用梯形公式：上用梯形公式：=Tn/*中值定理中值定理*/二、复化辛普森公式二、复化辛普森公式:44444=Sn注：注：为方便编程，可采用另一记法：令为方便编程，可采用另一记法：令 n=2n 为偶数，为偶数，这时这时 ，有，有三、收敛速度与误差估计：三、收敛速度与误差估计：若若一一个个积积分分公公式式的的误误差差满满足足 且且C 0，则则称该公式是称该公式是 p 阶收敛阶收敛的。的。例例4：计算计算解：解：其中其中=3.138988494其中其中=3.141592502运算量基运算量基本相同本相同问题问题:给定精度给定精度 ，如何取，如何取 n?例如：要求例如：要求 ，如何判断，如何判断 n=？?上述上述例例4中若要求中若要求 ,则则即：取即：取 n=409通常采取将区间通常采取将区间不断对分不断对分的方法，即取的方法，即取 n=2k上述上述例例4中中2k 409 k=9 时，时，T512=3.14159202例例4中：中：S4=3.141592502注意到区间再次对分时注意到区间再次对分时可用来判断迭代可用来判断迭代是否停止。是否停止。(1)(2)(3)事后误差估计事后误差估计一、梯形法的递推化一、梯形法的递推化逐次分半法逐次分半法 上一节介绍的复化求积方法对提高精度是行之有效的，但上一节介绍的复化求积方法对提高精度是行之有效的，但在使用求积公式之前必须给出合适的步长，在使用求积公式之前必须给出合适的步长，步长步长取得取得太大精度太大精度难以保证难以保证，步长太小步长太小则会导致则会导致计算量计算量的的增加增加，而事先给出一个，而事先给出一个恰当的步长又往往是困难的恰当的步长又往往是困难的 实际计算中常常实际计算中常常采用变步长的计算方案采用变步长的计算方案，即在步长，即在步长逐次分逐次分半半(即步长二分即步长二分)的过程中，反复利用复化求积公式进行计算，的过程中，反复利用复化求积公式进行计算，直至所求得的积分值满足精度要求为止直至所求得的积分值满足精度要求为止 设将求积区间设将求积区间a，b分成分成n等分，则一共有等分，则一共有n+1个分点，按个分点，按梯形公式计算积分值梯形公式计算积分值Tn，需要提供需要提供n+1个函数值如果将求积个函数值如果将求积区间再二分一次，则分点增至区间再二分一次，则分点增至2n+1个，我们来个，我们来考察考察二分二分前后两前后两个积分值个积分值之间的之间的联系联系4 4 龙贝格求积公式龙贝格求积公式逐次分半逐次分半计算计算方案方案的实现的实现:注意到每个子区间注意到每个子区间xk，xk+1经过二分只增加了一个分经过二分只增加了一个分点点 xk+1/2(xk+xk+1)/2，用复化梯形公式求得该子区间上的用复化梯形公式求得该子区间上的积分值为积分值为这里这里 代表二分前的步长代表二分前的步长.将每个子区间上的积分值将每个子区间上的积分值相加得相加得 二、龙贝格算法二、龙贝格算法根据复化梯形公式的余项表达式根据复化梯形公式的余项表达式可见，可见，利用利用两种步长两种步长计算的结果能估计截断误差计算的结果能估计截断误差.若将该截断若将该截断误差加到计算结果中误差加到计算结果中，就得出就得出“改进的梯形求积公式改进的梯形求积公式”：事后误差事后误差估计估计例：例：计算计算已知对于已知对于 =10 6 须将区间对分须将区间对分 9 次，得到次，得到 T512=3.14159202由由 来计算来计算 I 效果是否好些效果是否好些？=3.141592502=S4改进梯形求积公式改进梯形求积公式的右边实际是的右边实际是这就是说用这就是说用梯形法二分前后的两个积分值梯形法二分前后的两个积分值Tn与与T2n的的线性组合线性组合的结果的结果得到得到复化辛普森法求积公式复化辛普森法求积公式 类似的情况，用辛普森法二分前后的两个积分值类似的情况，用辛普森法二分前后的两个积分值Sn与与S2n的线性组合的结果可得到的线性组合的结果可得到复化柯特斯求积公式复化柯特斯求积公式 重复同样的手续，用柯特斯法二分前后的两个积分值重复同样的手续，用柯特斯法二分前后的两个积分值Cn与与C2n的线性组合的结果可得到的线性组合的结果可得到龙贝格龙贝格(Romberg)求积公式求积公式 我们在变步长的过程中运用加速公式，就能将粗糙的梯我们在变步长的过程中运用加速公式，就能将粗糙的梯形值形值Tn逐步加工成精度较高的辛普森值逐步加工成精度较高的辛普森值Sn、柯特斯值柯特斯值Cn和龙和龙贝格值贝格值Rn.一般有：一般有：Romberg 算法：算法：?T1=)0(0T T8=)3(0T T4=)2(0T T2=)1(0T S1=)0(1T R1=)0(3T S2=)1(1T C1=)0(2T C2=)1(2T S4=)2(1TRomberg 序列序列0 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 区间等分数区间等分数 梯形序列梯形序列 辛普森序列辛普森序列 柯特斯序列柯特斯序列 龙贝格序列龙贝格序列 龙贝格求积算法可用下表来表示：龙贝格求积算法可用下表来表示：例例5 用龙贝格方法计算椭圆用龙贝格方法计算椭圆 x2/4+y2 l 的周长，使结果的周长，使结果具有五位有效数字具有五位有效数字 分析分析 为便于计算，先将椭圆方程采用参数形式表示为便于计算，先将椭圆方程采用参数形式表示,再根再根据弧长公式将椭圆周长用积分形式表示由于计算结果要求具据弧长公式将椭圆周长用积分形式表示由于计算结果要求具有五位有效数字，因此需要估计所求积分值有几位整数，从而有五位有效数字，因此需要估计所求积分值有几位整数，从而确定所求积分值的绝对误差限最后再应用龙贝格方法计算积确定所求积分值的绝对误差限最后再应用龙贝格方法计算积分分 解解 令令 x 2cos，y sin ，则椭圆的周长为则椭圆的周长为下表给出了用龙贝格方法计算积分下表给出了用龙贝格方法计算积分I=1+1+3sin2 dx 的过程的过程./200 1 2.356 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)外推加速法外推加速法 上面讨论说明由梯形公式出发上面讨论说明由梯形公式出发,将区间将区间a,b逐次二分逐次二分可提高求积公式的精度可提高求积公式的精度,上述加速过程还可继续下去上述加速过程还可继续下去.下面我们讨论其下面我们讨论其理论依据理论依据.若记若记Tn=T(h),当区间当区间a,b分为分为2n等分时等分时,有有 ,则则可见可见I=T(h)的误差为的误差为O(h2).若记若记 ，则，则 将梯形公式按余项展开将梯形公式按余项展开.由误差公式有由误差公式有显然显然T1(h)与与 I 近似的阶为近似的阶为O(h4).就是就是辛普森公式辛普森公式序列序列Sn,S2n,.这样构造的这样构造的 则又可进一步从余项中则又可进一步从余项中消去消去 h4 项，这样构造出的项，这样构造出的 ，其实就是，其实就是柯特斯公式柯特斯公式序序列，它与列，它与 I 的逼近阶为的逼近阶为O(h6).若令若令 ,一般地，若记一般地，若记T0(h)=T(h)，经过经过m(m=1,2,)次加速次加速后，则有后，则有如此继续下去，每加速一次，误差的量级便提高如此继续下去，每加速一次，误差的量级便提高2阶阶.Romberg 算法算法 可以证明，如果可以证明，如果 f(x)充分光滑，那么充分光滑，那么T 数表每一列的数表每一列的元素及对角线元素均收敛到所求的积分值元素及对角线元素均收敛到所求的积分值 I，即即 理查德森理查德森外推法外推法利用利用低低阶公式产生阶公式产生高高精度的结果。精度的结果。设设对于某一对于某一 h 0，有公式有公式 T0(h)近似计算某一未知值近似计算某一未知值 I。由。由Taylor展开得到：展开得到：T0(h)I=1 h+2 h2+3 h3+i 与与 h 无关无关现将现将 h 对分，得：对分，得：()()().)(3232222120+=hhhhIT Q：如何将公式精度由如何将公式精度由 O(h)提高到提高到 O(h2)？.432112)()(23322020 =hhIhTTh 即：即：HW:p.159 6,7,8.3)在构造在构造Newton-Cotes公式公式时，限定用积分区间的时，限定用积分区间的等分点等分点作为求积节点作为求积节点，这样做虽然使问题的处理过程得以简化，但，这样做虽然使问题的处理过程得以简化，但同时也同时也限制了精度限制了精度。求积公式含有求积公式含有2n+2个待定参数个待定参数xk、Ak(k0,1,n)若用若用待定系数法确定它们待定系数法确定它们,则最好需要则最好需要2n+2个独立的条件联立方个独立的条件联立方程组求解程组求解,从而易知求积公式的从而易知求积公式的最大代数精度最大代数精度可达到可达到2n+1次次.在节点数目固定为在节点数目固定为n 的条件下，能否通过的条件下，能否通过适当选取求积适当选取求积节点节点xk的位置以及相应的求积系数的位置以及相应的求积系数Ak，使求积公式，使求积公式具有尽可能高具有尽可能高(最高最高)的代数精度？的代数精度？这类求积公式称为这类求积公式称为高斯高斯(Gauss)求积公式求积公式。4 高斯求积公式高斯求积公式 将节点将节点 x0 xn 以及系数以及系数 A0 An 都都作为待定系数。作为待定系数。令令 f(x)=1,x,x2,x2n+1 代入可求解，得到的公式代入可求解，得到的公式具有具有2n+1 次代数精度。这样的节点称为次代数精度。这样的节点称为Gauss 点点，公式称为公式称为Gauss 型求积公式型求积公式。为使问题更具一般性为使问题更具一般性,我们研究带权积分我们研究带权积分 (x)为权函数为权函数,Ak(k0,1,n)为不依赖于为不依赖于f(x)的求积系的求积系数数,xk(k0,1,n)为求积节点为求积节点.要使要使(5.1)具有具有2n+1次代数精度，则需要满足次代数精度，则需要满足构造具有构造具有2n+1次代数精度的求积公式次代数精度的求积公式(5.1)从例中可看到求解非线性方程组从例中可看到求解非线性方程组(5.2)较复杂，通常较复杂，通常n2就很难求解故一般不通过解方程就很难求解故一般不通过解方程(5.2)求求 xk 及及 Ak(k0,1,n)例：例：求求 的的 2 点点 Gauss 公式。公式。解：解：设设 ，应有，应有 3 次代数精度。次代数精度。+101100)()()(xfAxfAdxxfx代入代入 f(x)=1,x,x2,x3不是线性方程组，不是线性方程组，不易求解。不易求解。而从研究而从研究高斯点的基本特性高斯点的基本特性来着手解决来着手解决Gauss 求积公式求积公式的构造问题的构造问题由插值型公式构由插值型公式构造知造知,关键求关键求xk，一、高斯点的基本特性一、高斯点的基本特性 证明：证明：“”x0 xn 为为 Gauss 点点,则公式则公式 至少有至少有 2n+1 次代数精度。次代数精度。对任意次数对任意次数不大于不大于n 的多项式的多项式 Pm(x)，Pm(x)w(x)的次的次数数不大于不大于2n+1，则代入公式应则代入公式应精确成立精确成立：0=0“”要证明要证明 x0 xn 为为 Gauss 点，即要证公式对任意点，即要证公式对任意次数次数不大于不大于2n+1 的多项式的多项式 Pm(x)精确成立，即证明：精确成立，即证明：设设0 x0 xn 为为 Gauss 点点 与与任任意意次次数数不大于不大于n 的多项式的多项式 P(x)（带权）正交带权）正交。定理定理求求 Gauss 点点 求求w(x)的的零点零点 Gauss 公式的余项：公式的余项：/*设设P为为f 的过的过x0 xn的插值多项式的插值多项式*/*只要只要P 的阶数不大于的阶数不大于2n+1，则下一步则下一步等式成立等式成立*/插值多项式插值多项式的余项的余项Q：什么样的什么样的插值多项式插值多项式在在 x0 xn 上上有有 2n+1 阶？阶？A：Hermite 多项式！多项式！满足满足二、高斯求积公式的余项二、高斯求积公式的余项三、高斯求积公式的稳定性与收敛性三、高斯求积公式的稳定性与收敛性 定理定理6 高斯求积公式高斯求积公式(5.1)的求积系数的求积系数 Ak(k0,1,n)全是正的全是正的 由本定理及定理由本定理及定理2，则得，则得 推论推论 高斯求积公式高斯求积公式(5.1)是稳定的是稳定的.定理定理7 设设 f(x)C a，b，则高斯求积公式则高斯求积公式(5.1)是收敛是收敛 的，即的，即 正交多项式族正交多项式族 0,1,n,有性质：任意次数不有性质：任意次数不大于大于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即：即：四、常用的高斯型求积公式四、常用的高斯型求积公式Step 2：求求 2=0 的的 2 个根，即为个根，即为 Gauss 点点 x0，x1Step 3：代入代入 f(x)=1,x 以求解以求解 A0，A1解解线性线性方程组，方程组，简单。简单。结果与前一方法相同：结果与前一方法相同：利用此公式计算利用此公式计算 的的值值注：注：构造正交多项式也可以利用最小二乘数据拟合中介绍构造正交多项式也可以利用最小二乘数据拟合中介绍过的递推式进行。过的递推式进行。特殊正交多项式族：特殊正交多项式族：Legendre 多项式族：多项式族：1)(x 定义在定义在 1,1上，上，满足：满足：由由 有递推有递推以以 Pn+1 的根为节点的求积公式称为的根为节点的求积公式称为Gauss-Legendre 公式公式。Chebyshev 多项式族：多项式族：211)(xx=定义在定义在 1,1上，上，Tn+1 的根的根为为k=0,n以此为节点构造公式以此为节点构造公式称为称为 Gauss-Chebyshev 公式公式。注意到积分端点注意到积分端点 1 可能是积分可能是积分的的奇点奇点，用普通，用普通Newton-Cotes公公式在端点会出问题。而式在端点会出问题。而Gauss公公式可能避免此问题的发生。式可能避免此问题的发生。其它公式见教材其它公式见教材p.144-148注：注：一般一般a,b上的积分可化为上的积分可化为-1,1上特殊高斯公式进行计算。上特殊高斯公式进行计算。5 5 数值微分数值微分数值微分的数值微分的概念概念数值微分的数值微分的计算方法计算方法原始概念近似原始概念近似:中点法及外推法中点法及外推法函数近似函数近似:插值型的求导公式插值型的求导公式函数相互关系转化函数相互关系转化:利用数值积分求导利用数值积分求导数值微分的数值微分的误差分析误差分析泰勒展式估计泰勒展式估计事后误差估计事后误差估计基本关系转化基本关系转化 数值微分数值微分就是就是用函数值的线性组合近似函数在某点的用函数值的线性组合近似函数在某点的导数值导数值一、中点法和外推法一、中点法和外推法 按导数定义按导数定义,是差商是差商 当当 时的时的极限极限取取差商差商作为作为导数导数的近似值的近似值，建立简单的数值微分方法，建立简单的数值微分方法:(6.1)向后差商近似导数向后差商近似导数(6.2)(6.3)中心差商近似导数中心差商近似导数容易看出，就精度而言，以（容易看出，就精度而言，以（6.3）式更为可取，称）式更为可取，称(6.4)为为 的的中点公式中点公式,其中其中h为一增量，称为为一增量，称为步长步长 这种数这种数值微分方法称为值微分方法称为中点方法中点方法,它是前两种方法的算术平均它是前两种方法的算术平均分别将分别将在在 x=a 处做处做Taylor展开有展开有代入代入D(h)得得其中其中现在来考虑中点公式现在来考虑中点公式 的截断误差，的截断误差，(6.5)所以所以截断误差截断误差(6.6)从截断误差的角度来看，步长从截断误差的角度来看，步长h越小，计算结果越准确。且越小，计算结果越准确。且 但从计算角度看，但从计算角度看，h 越小，越小，f(a+h)与与 f(a-h)越接近，直越接近，直接相减会造成有效数字的严重损失接相减会造成有效数字的严重损失(参看第参看第1章第章第4节节)。因此。因此,从舍入误差的角度来看，步长从舍入误差的角度来看，步长 h 不宜太小。不宜太小。所以所以,在在实际计算时实际计算时，通常，通常采用采用二分步长二分步长及及误差事后误差事后估计法估计法,在变步长的过程中实现步长的自动选择，在保证在变步长的过程中实现步长的自动选择，在保证截断误差满足的精度要求的前提下选取取尽可能大的步长。截断误差满足的精度要求的前提下选取取尽可能大的步长。kD(h)01239103.017652.791352.736442.722812.718282.71828 解解 这里采用的计算公式是这里采用的计算公式是 计算结果见表计算结果见表6.1，表中，表中 k 代表二分的次代表二分的次数，步长数，步长 。二分。二分 9 次得结果次得结果 D=2.71828，它的每一数字都是有效数字，它的每一数字都是有效数字(所所求导数的准确值为求导数的准确值为e=2.7182818)。M 例例6.1 用用变步长变步长的中点方法求的中点方法求 在在x=1处的导数值处的导数值,设取设取h=0.8起算。起算。表表6.1 计算结果计算结果 事后误差估计事后误差估计 我们看到，中点公式具有如下形式我们看到，中点公式具有如下形式 (6.7)式中的系数均与步长无关。式中的系数均与步长无关。若将步长二分，则有若将步长二分，则有 (6.8)若令若令 (6.10)则进一步消去则进一步消去D1(h)误差主项误差主项 ,有有 则可消去误差主项则可消去误差主项 ，得，得 取取(6.7)与与(6.8)加权平均加权平均(6.9)重复同样的手续重复同样的手续，再导出下列加速公式，再导出下列加速公式 这种加速过程还可继续下去。这种加速过程还可继续下去。这种加速方法通常称作这种加速方法通常称作Richardson外推加速法外推加速法。例例6.2 运用加速公式和加工例运用加速公式和加工例6.1的结果。的结果。表表6.2 Richardson外推加速法计算结果外推加速法计算结果h D(h)0.80.40.20.13.017652.791352.736442.722812.7159172.7181372.7182672.7182852.7182762.71828 解解 计算结果计算结果见表见表6.2。这里，。这里，加速的效果同样加速的效果同样是相当显著的。是相当显著的。最小最小,步长步长h不宜太大不宜太大,也不宜太小也不宜太小.其其最优步长最优步长应为应为 它表明它表明h 越小越小,舍入误差舍入误差越大越大,故它是病态的故它是病态的.用中点用中点公式计算公式计算的的误差上界误差上界为为要使误差要使误差E(h)则计算则计算 当当 f(a+h)及及 f(a-h)分别有舍入误差分别有舍入误差 1 及及 2 时时，若令若令的舍入误差上界为的舍入误差上界为 注注:中点心公式及其加速方法适合用表达式表示的函数。对中点心公式及其加速方法适合用表达式表示的函数。对于列表函数，则宜使用插值方法等导出数值求导公式。于列表函数，则宜使用插值方法等导出数值求导公式。二、插值型的求导公式二、插值型的求导公式 x x0 x1 x2 xn y y0 y1 y2 yn 对于列表函数对于列表函数 y=f(x)：插值多项式插值多项式 y=Pn(x)作为它的近似作为它的近似,我们取我们取统称统称插值型的求导公式插值型的求导公式的近似值，的近似值，作为作为建立的数值公式建立的数值公式 依据插值余项定理，依据插值余项定理，求导公式求导公式(6.11)的的余项余项为为式中式中(6.11)我们限定我们限定:求某个节点求某个节点 xk 上的导数值，上的导数值，上面的第二项变为零，上面的第二项变为零，这时有余项公式这时有余项公式(6.12)1 1两点公式两点公式 已给两节点已给两节点 x0,x1 上的函数值上的函数值 f(x0),f(x1)，做线性插值做线性插值记记 x1 x0=h，对上式两端求导，有，对上式两端求导，有于是有下列求导公式：于是有下列求导公式：而利用余项公式而利用余项公式(6.4)知，带知，带余项的两点公式余项的两点公式是：是：下面我们仅仅考察下面我们仅仅考察节点处的导数值节点处的导数值为简化讨论为简化讨论,假定假定所给的节点是等距的所给的节点是等距的 设已给出三节点设已给出三节点x0,xl=x0+h,x2=x0+2h上的函数值上的函数值,做二次插值做二次插值令令 x=x0+th,则则 2三点公式三点公式(6.13)这里撇号这里撇号()表示对表示对x求导求导上式分别取上式分别取 t=0,1,2,得到三种三点公式：得到三种三点公式：(6.14)而而带余项的三点求导公式带余项的三点求导公式如下：如下：公式公式(6.14)是我们所熟悉的是我们所熟悉的中点公式中点公式在三点公式中，它由于少用了一在三点公式中，它由于少用了一个函数值个函数值 f(x1)而引人注目而引人注目 例如，将式例如，将式(6.13)再对再对t求导一次，有求导一次，有于是有于是有 用插值多项式用插值多项式Pn(x)作为作为 f(x)的近似函数，还可以建立的近似函数，还可以建立高阶数值微分高阶数值微分公式：公式：而带余项的而带余项的二阶三点公式二阶三点公式如下：如下：(6.15)例例 已知函数已知函数y=ex 的下列数值的下列数值:x 2.5 2.6 2.7 2.8 2.9 y 12.1825 13.4637 14.8797 16.4446 18.1741试用二点、三点微分公式计算试用二点、三点微分公式计算x=2.7处的一阶、二阶导数值处的一阶、二阶导数值.本题没有明确指出用哪些点处的函数值来求本题没有明确指出用哪些点处的函数值来求 f(2.7)和和f(2.7).因此，随着步长因此，随着步长h不同不同,导数值有可能不同导数值有可能不同.另外另外,用两点函用两点函数值时，只能求一阶导数值数值时，只能求一阶导数值.解解:方法方法1:取取h=0.1时，两点公式有两种取法当时，两点公式有两种取法当 x0=2.6，x12.7时，时，f(2.7)f(2.7)-f(2.6)=14.8797-13.4637 =14.16000.110.61 当当x 0=2.7，x1=2.8时，时，f(2.7)f(2.8)-f(2.7)=15.64900.11 方法方法2:取取h=0.2时，两点公式有两种取法当时，两点公式有两种取法当 x0=2.5，x12.7时，时，f(2.7)f(2.7)-f(2.5)=13.4860当当 x0=2.7，x12.9时，时，f(2.7)f(2.9)-f(2.7)=16.47200.210.21三点公式取三点公式取x0=2.6，x1=2.7，x2=2.8，得，得 f(2.7)f(2.8)-f(2.6)=14.9045 f(2.7)f(2.8)-2f(2.7)+f(2.6)=14.89000.1212 0.1 1三点公式取三点公式取x0=2.5，x1=2.7，x2=2.9，得，得 f(2.7)f(2.9)-f(2.5)=14.9790 f(2.7)f(2.9)-2f(2.7)+f(2.5)=14.93000.2212 0.2 1 注注:f(2.7)和和 f(2.7)的真值都是的真值都是 14.87973，上面的计，上面的计算表明算表明:1）当使用两点公式时应取步长较小的函数值；）当使用两点公式时应取步长较小的函数值；2）一般情况下，同样步长的两点公式没有三点公式准确）一般情况下，同样步长的两点公式没有三点公式准确,步长越小越精确但如果高阶导数无界或舍入误差超过截断误步长越小越精确但如果高阶导数无界或舍入误差超过截断误差时差时,这个结论就不一定对了这个结论就不一定对了.附注附注:与积分相比，数值微分比较困难与积分相比，数值微分比较困难。积分描述了一个函数的整体或宏观性质，而微分则描述一个函积分描述了一个函数的整体或宏观性质，而微分则描述一个函数在一点处的斜率，这是函数的微观性质。因此积分对函数的形状数在一点处的斜率，这是函数的微观性质。因此积分对函数的形状在小范围内的改变不敏感。而微分却很敏感。在小范围内的改变不敏感。而微分却很敏感。一个函数小的变化，一个函数小的变化，容易产生相邻点的斜率的大的改变容易产生相邻点的斜率的大的改变。由于微分这个固有的困难，所以应尽可能避免数值微分，特别由于微分这个固有的困难，所以应尽可能避免数值微分，特别是对实验获得的数据进行微分。在这种情况下，最好用最小二乘曲是对实验获得的数据进行微分。在这种情况下，最好用最小二乘曲线拟合这种数据，然后对所得到的多项式进行微分。或用另一种方线拟合这种数据，然后对所得到的多项式进行微分。或用另一种方法，对该数据进行三次样条拟合，然后寻找该样条函数的微分；等法，对该数据进行三次样条拟合，然后寻找该样条函数的微分；等等。等。一般是先拟合（或逼近），再微分一般是先拟合（或逼近），再微分。三、利用数值积分求导三、利用数值积分求导 微分是积分的逆运算，因此可利用数值积分的方法来计算微分是积分的逆运算，因此可利用数值积分的方法来计算数值微分设数值微分设 f(x)是一个充分光滑的函数，设是一个充分光滑的函数，设则有：则有：对上式中积分用不同的求积公式对上式中积分用不同的求积公式,再从中求出再从中求出(x)就可得不就可得不同的微分公式同的微分公式.如如:对对用中矩形公式用中矩形公式则得则得从而得到从而得到中点微分公式中点微分公式若对若对(6.8)右端积分用辛普森求积公式，则有右端积分用辛普森求积公式，则有上式略去余项，并记上式略去余项，并记普森数值微分公式普森数值微分公式的近似值为的近似值为mk，则得到则得到辛辛这是关于这是关于m0，m1，mn这这 n+1 个未知量的个未知量的 n-1 个方程组，个方程组，可用求方程的方法求解。可用求方程的方法求解。