第4章 数值积分与数值微分
引言 数值求积的基本思想 微积分基本定理,对于积分 只要找到被积函数 的原函数 ,便有下列牛顿-莱 只要找到被积函数 的原函数 ,便有下列牛顿-莱 布尼兹(Newton-Leibniz)公式: 但对于下列情形:
(1)被积函数,诸如 等等,找不到用 初等函数表示的原函数; (2)当 是由测量或数值计算给出的一张数据表 时,牛顿-莱布尼兹公式也不能直接运用. 因此有必要研究积分的数值计算问题. 由积分中值定理知,在积分区间 内存在一点ξ, 成立
即:底为 而高为 的矩形面积恰等于所求 曲边梯形的面积 (图4-1). 图4-1
问题在于点ξ的具体位置一般是不知道的,因而难以 准确算出 的值. 将 称为区间 上的平均高度. 只要对平均高度 提供一种算法,相应地便 获得一种数值求积方法. 用两端点“高度” 与 的算术平均作为平均高度 的近似值,这样导出的求积公式 (1.1) 是梯形公式(几何意义参看图4-2).
图4-2 用区间中点 的“高度” 近似地取代平均 高度 ,则又可导出所谓中矩形公式(简称矩形公式) (1.2)
一般地,可以在区间 上适当选取某些节点 , 然后用 加权平均得到平均高度 的近似值, 这样构造出的求积公式具有下列形式: (1.3) 式中 称为求积节点; 称为求积系数, 亦称伴随节点 的权. 权 仅仅与节点 的选取有关, 而不依赖于被积函数 的具体形式.
这类数值积分方法通常称为机械求积,其特点是将积 分求值问题归结为函数值的计算,这就避开了牛顿-莱布尼 兹公式需要寻求原函数的困难.
代数精度的概念 数值求积是近似方法,为保证精度,自然希望求积公 式对尽可能多的函数准确成立. 定义1 如果某个求积公式对于次数不超过 的多项式 均能准确地成立,但对于 次多项式就不准确成立, 则称该求积公式具有 次代数精度. 梯形公式(1.1)和矩形公式(1.2)均具有一次代数精度.
欲使求积公式(1.3)具有 次代数精度,则只要令它 对 都准确成立,就得到 (1.4)
如果事先选定求积节点 ,譬如,以区间 的等 距分点作为节点,这时取 ,求解方程组(1.4)即可确 定求积系数 ,而使求积公式(1.3)至少具有 次代数精度. 构造形如(1.3)的求积公式,原则上是一个确定参数 和 的代数问题.
插值型的求积公式 设给定一组节点 且已知函数 在这些节点上的值, 做插值函数 . 取 作为积分 的近似值, 这样构造出的求积公式 (1.5)
称为是插值型的,式中求积系数 通过插值基函数 积分得出 (1.6) 由插值余项定理(第2章的定理2)即知,对于插值型的 求积公式(1.5),其余项 (1.7) 式中ξ与变量 有关,
当 是次数不超过 的多项式时,插值多项式就是 函数本身, 余项 为零, 所以这时插值型求积公式 至少具有 次代数精度. 反之,如果求积公式(1.5)至少具有 次代数精度,则 它必定是插值型的. 事实上,这时公式(1.5)对于插值基函数 应准确 成立,即有
注意到 上式右端实际上即等于 , 因而式(1.6) 成立. 这样,有 定理1 形如(1.5)的求积公式至少有 次代数精度的 充分必要条件是,它是插值型的.
求积公式的收敛性与稳定性 定义2 在求积公式(1.3)中,若 其中 则称求积公式(1.3)是收敛的. 在求积公式(1.3)中,由于计算 可能产生误差 , 实际得到将是 , 即 记
如果对任给小正数 只要误差 充分小就有 (1.8) 则表明求积公式(1.3)计算是稳定的, 由此给出: 定义3 对任给 若 只要 就有(1.8)成立,则称求积公式(1.3)是稳定的.
定理2 若求积公式(1.3)中系数 则此求积公式是稳定的. 证明 对任给 取 若对 都有 则当 时有
由定义3 ,知求积公式(1.3)是稳定的.
4.2 牛顿-柯特斯公式 柯特斯系数 设将积分区间 划分为 等分, 步长 选取等距节点 构造出的插值型求积公式 (2.1) 4.2 牛顿-柯特斯公式 柯特斯系数 设将积分区间 划分为 等分, 步长 选取等距节点 构造出的插值型求积公式 (2.1) 称为牛顿-柯特斯公式, 式中 称为柯特斯系数. 按(1.6)式,引进变换 则利用等距节点的 插值公式,有
柯特斯系数 (2.2) 当 时, 这时的求积公式就是梯形公式(1.1)
辛普森公式 当 时,按(2.2)式, 柯特斯系数为 相应的求积公式是辛普森(Simpson)公式 (2.3)
的牛顿-柯特斯公式称为柯特斯公式, 其形式是 (2.4) 这里 按(2.2)式,可构造柯特斯系数表.
从柯特斯系数表看到 时,柯特斯系数 出现 负值, 于是有 特别地,假定 且 则有
它表明初始数据误差将会引起计算结果误差增大,即计算 不稳定,故 的牛顿-柯特斯公式是不用的.
4.2.2 偶阶求积公式的代数精度 由定理1, 阶的牛顿-柯特斯公式至少具有 次的代数精度. 本节讨论代数精度的进一步提高问题. 4.2.2 偶阶求积公式的代数精度 由定理1, 阶的牛顿-柯特斯公式至少具有 次的代数精度. 本节讨论代数精度的进一步提高问题. 先看辛普森公式(2.3),它是二阶牛顿-柯特斯公式,因 此至少具有二次代数精度. 用 进行检验, 按辛普森公式计算得
另一方面,直接求积得 这时有 , 即辛普森公式对次数不超过三次的多项式 均能准确成立, 而它对 通常是不准确的, 因此, 辛普森公式实际上具有三次代数精度. 定理3 当阶 为偶数时,牛顿-柯特斯公式(2.1)至少 有 次代数精度.
证明 我们只要验证,当 为偶数时,牛顿-柯特斯 公式对 的余项为零. 由于这里 按余项公式(1.7) 有 引进变换 并注意到 有
若 为偶数,则 为整数, 再令 进一步有 因为被积函数 为奇函数,所以
4.2.3 几种低阶求积公式的余项 按余项公式(1.7),梯形公式(1.1)的余项 这里积分的核函数 在区间[a,b]上保号(非正), 4.2.3 几种低阶求积公式的余项 按余项公式(1.7),梯形公式(1.1)的余项 这里积分的核函数 在区间[a,b]上保号(非正), 应用积分中值定理,在 [a,b]内存在一点 使 (2.5)
为研究辛普森公式(2.3)的余项 Rs=I-S 构造次数 不超过3的多项式H3(x)满足 H3(c)=f(c), H3’(c)=f’(c) H3(a)=f(a), H3(b)=f(b) (2.6) 其中 辛普森公式具有三次代数精度,对于这样构造出的三次式H(x)应是准确的,即
由插值条件(2.6),上式右端实际上等于按辛普森公式(2.3) 求得的积分值 , 因此积分余项 对于多项式 ,其插值余项由第2章(5.11)得 故有
这时积分的核函数 在[a,b]上保号 (非正),再用积分中值定理有 (2.7) 类似的,对于柯特斯公式(2.4),结果如下:
4.3 复化求积公式 复化梯形公式 复化求积的基本思想是把积分区间分成若干子区间(通 常是等分),再在每个子区间上用低阶求积公式,目的是提 4.3 复化求积公式 复化求积的基本思想是把积分区间分成若干子区间(通 常是等分),再在每个子区间上用低阶求积公式,目的是提 高精度. 复化梯形公式 将区间[a,b]划分为 等分, 分点 在每个子区间[xk, xk+1] (k=0,1,…,n-1)上 采用梯形公式(1.1),则得
(3.1) 记 (3.2) 称为复化梯形公式.
由(2.5) ,其余项 由于 , 且 所以 使 于是复化梯形公式余项为
(3.3) 误差是h2阶, 且当 时有 即复化梯形公式是收敛的. 将Tn 改写为
只要 则当 时,上式右端括号内的两个 和式均收敛到积分 所以复化梯形公式(3.2)收敛. 此外, 的求积系数为正,由定理2知复化梯形公式是 稳定的.
4.3.2 复化辛普森求积公式 将区间[a,b]分为 等分, 在每个子区间[xk,xk+1]上 采用辛普森公式(2.3), 若记 则得 (3.4) 记 (3.5)
复化辛普森求积公式 称为复化辛普森求积公式. 由(2.7),其余项 于是当 时, 与复化梯形公式相似有 (3.6) 于是当 时, 与复化梯形公式相似有 (3.6) 误差阶为h4,显然是收敛的.
实际上,只要 则可得到收敛性, 即 此外,由于Sn中求积系数均为正数,故知复化辛普森 公式计算稳定.
例1 对于函数 , 给出 的函数表 (见表4-2), 试用复化梯形公式(3.2)及复化辛普森公式(3.5) 计算积分 并估计误差. 解 将积分区间[0,1]划分为8等分, 应用复化梯形法求得
而如果将[0,1]分为4等分,应用复化辛普森法有 以上得到的两个结果T8与S4,都需要提供9个点上的 函数值, 计算量基本相同,然而精度却差别很大. 同积分的准确值 比较, 复化梯形法的结果 只有两位有效数字. 接下来看误差估计 ,由于 所以有
于是 由(3.3)得复化梯形公式误差
对复化辛普森公式,由(3.6)得
4.4 龙贝格求积公式 梯形法的递推化 复化求积方法可提高求积精度,实际计算时若精度不 够可将步长逐次分半. 4.4 龙贝格求积公式 梯形法的递推化 复化求积方法可提高求积精度,实际计算时若精度不 够可将步长逐次分半. 设将区间 分为 等分,共有 个分点, 如果将求积区间再二分一次,则分点增至 个, 我们将二分前后两个积分值联系起来加以考察.
每个子区间[xk,xk+1],经过二分只增加了一个分点 用复化梯形公式求得该子区间上的积分值为 这里 代表二分前的步长. 将每个子区间上的积分值相加得
从而利用式(3.2)可导出下列递推公式 (4.1) 例2 计算积分值 解 先对整个区间[0,1]使用梯形公式. 对于函数 定义它在 的值 而 由梯形公式
将区间二等分,求出中点的函数值 利用递推公式(4.1),有 进一步二分求积区间,并计算新分点上的函数值 再利用式(4.1),有
这样不断二分下去,计算结果见下表. 它表明用复化梯形公式计算积分 要达到7位有效数 字的精度需要二分区间10次,即要有分点1025个,计算量 很大.
4.4.2 龙贝格算法 梯形法计算简单但收敛慢,本节讨论如何提高收敛速 度以节省计算量. 根据复化梯形公式的余项表达式(3.3) 假定 则有
移项整理,得 (4.2) 由此可见,只有二分前后的两个积分值 与 相 当接近,就可以保证计算结果 的误差很小. 这样直接用计算结果来估计误差的方法通常称作误差的事后估计法. 按式(4.2),积分近似值 的误差大致等于 若用这个误差值作为 的一种补偿,可以期望所得到的
(4.3) 可能是更好的结果. 直接验证知,按公式(4.3)组合得到的近似值恰为 即 (4.4) 也就是说,利用梯形法二分前后的两个积分值 与 , 按式(4.3)做线性组合,结果得到辛普森法的积分值
再考察辛普森法,按误差公式(3.6),其截断误差大致 与 成正比, 因此,若将步长折半则误差将减至原有误差 的1/16,即有 由此得到 不难直接验证,上式右端的值等于 , 为复化柯特斯公式, 它的精度为
就是说,用辛普森法二分前后两个积分值 与 , 可得到 的近似误差估计为 并且按上式做线性组合可得复化柯特斯的积分值 , 即有 (4.5) 重复同样的手续,依据柯特斯法的误差阶为 , 可进一步导出下列龙贝格(Romberg)公式: (4.6)
在变步长的过程中运用公式(4.4),(4.5)和(4.6), 就能将粗糙的梯形值 逐步加工成精度较高的辛普森值 柯特斯值 和龙贝格值 . 例3 用加速公式(4.4),(4.5)和(4.6)加工例2得到的 梯形值, 计算结果见下表( 代表二分次数):
可以看到,这里利用二分3次的数据(它们的精度都很差, 只有两三位有效数字),通过三次加速求得 这个结果的每一位数字都是有效数字,可见加速的效果是十分显著的.
4.4.3 理查森外推加速法 由梯形公式出发,将区间 逐次二分可提高求积 公式的精度,上述加速过程还可继续下去,设 若记 4.4.3 理查森外推加速法 由梯形公式出发,将区间 逐次二分可提高求积 公式的精度,上述加速过程还可继续下去,设 若记 当区间 分为 等分时, 则有 并且有 梯形公式余项可展成级数形式,即
定理4 设 则有 (4.7) 其中系数 与 无关. 定理4表明 是 阶, 在(4.7)中,若用 代替 ,有 (4.8) 若用4乘(4.8)式,减去(4.7)式再除3记之为 则得
(4.9) 这里 以及后面将出现的 均为与 无关的系数, 这样构造的 与积分值 近似的阶为 . 比较(4.9)与(4.4)可知,这样构造的序列 就是辛普森公式序列 又根据(4.9),有
若令 则又可进一步从余项展开式中消去 项,而有 这样构造出的 ,其实就是柯特斯公式序列. 它与积分值 的逼近阶为 如此继续下去,每加速一次,误差的量级便提高2阶.
一般地,若记 则有 (4.10) 经过 次加速后, 余项便取下列形式: (4.11) 上述处理方法通常称为理查森外推加速方法. 设以 表示二分 次后求得的梯形值,且以 表示 序列 的 次加速值,则依递推公式(4.10)可得
(4.12) 公式(4.12)也称为龙贝格求积算法. 计算过程: (1) 取 求 令 ( 记区间 的二分次数). (2) 求梯形值 即按递推公式(4.1)计算 (3) 求加速值,按公式(4.12)逐个求出T表(见表4-5)的 第 行其余各元素
(4) 若 (预先给定的精度),则终止计算, 并取 否则令 转(2)继续计算.
可以证明,如果 充分光滑,那么T表每一列的元素及对角线元素均收敛到所求的积分值 ,即 对于 不充分光滑的函数也可用龙贝格算法计算, 只是收敛慢一些,这时也可以直接使用复化辛普森公式计算.
例4 用龙贝格算法计算积分 解 在 上仅是一次连续可微, 用龙贝格算法计算结果见表4-6.
从表中看到用龙贝格算到 的精度与辛普森求积 精度相当. 这里 的精确值为