1、第四章数值积分和数值微分 4 1引言我们知道 若函数f x 在区间 a b 上连续且其原函数为F x 则可用Newton Leibnitz公式 求定积分的值 Newton Leibnitz公式无论在理论上还是在解决实际问题上都起了很大作用 但它并不能完全解决定积分的计算问题 因为积分学涉及的实际问题极为广泛 而且极其复杂 在实际计算中经常遇到以下三种情况 1 被积函数f x 并不一定能够找到用初等函数的有限形式表示的原函数F x 例如 Newton Leibnitz公式就无能为力了 3 被积函数f x 没有具体的解析表达式 其函数关系由表格或图形表示 对于这些情况 要计算积分的准确值都是十分困
2、难的 由此可见 通过原函数来计算积分有它的局限性 因而研究一种新的积分方法来解决Newton Leibniz公式所不能或很难解决的积分问题 这时需要用数值解法来建立积分的近似计算方法 将积分区间细分 在每一个小区间内用简单函数代替复杂函数进行积分 这就是数值积分的思想 用代数插值多项式去代替被积函数f x 进行积分是本章讨论数值积分的主要内容 同样对于函数f x 的求导问题 因为在微分学中 函数f x 的导数是通过极限定义的 若函数是以表格形式给出 或函数的表达式过于复杂时 也需要研究其数值计算方法 这是本章介绍的另一个内容 数值微分 4 2数值积分概述4 2 1数值积分的基本思想积分值在几何
3、上可以解释为由x a x b y 0以及y f x 这四条边所围成的曲边梯形面积 如图4 1所示 而这个面积之所以难于计算是因为它有一条曲边y f x 图4 1数值积分的几何意义 建立数值积分公式的途径比较多 其中最常用的有两种 1 由积分中值定理可知 对于连续函数f x 在积分区间 a b 内存在一点 使得即所求的曲边梯形的面积恰好等于底为 b a 高为的矩形面积 但是点 的具体位置一般是未知的 因而的值也是未知的 称为f x 在区间 a b 上的平均高度 那么只要对平均高度提供一种算法 相应地就获得一种数值求积方法 三个求积分公式 梯形公式 中矩形公式 按照这种思想 可构造出一些求积分值的
4、近似公式 例如分别取和 则分别得到中矩形公式和梯形公式 Simpson公式 f 的近似值而获得的一种数值积分方法 中矩形公式把 a b 的中点处函数值作为平均高度f 的近似值而获得的一种数值积分方法 Simpson公式是以函数f x 在a b a b 2这三点的函数值f a f b 的加权平均值似值而获得的一种数值积分方法 作为平均高度f 的近 2 先用某个简单函数近似逼近f x 用代替原被积函数f x 即 以此构造数值算法 从数值计算的角度考虑 函数应对f x 有充分的逼近程度 并且容易计算其积分 由于多项式能很好地逼近连续函数 且又容易计算积分 因此将选取为插值多项式 这样f x 的积分就
5、可以用其插值多项式的积分来近似代替 4 2 2插值求积公式设已知f x 在节点有函数值 作n次拉格朗日插值多项式 其中 称为求积系数 给出如下定义 设插值求积公式的余项为 由插值余项定理得 其中 当f x 是次数不高于n的多项式时 有R f 0 求积公式 4 4 能成为准确的等式 由于闭区间 a b 上的连续函数可用多项式逼近 所以一个求积公式能对多大次数的多项式f x 成为准确等式 是衡量该公式的精确程度的重要指标 为此给出以下定义 定义4 2 代数精度 设求积公式 4 4 对于一切次数小于等于m的多项式 是准确的 而对于次数为m 1的多项式是不准确的 则称该求积公式具有m次代数精度 简称代
6、数精度 由定义可知 若求积公式 4 4 的代数精度为n 则求积系数应满足线性方程组 或 这是关于的线性方程组 其系数矩阵 是范得蒙矩阵 当互异时非奇异 故有唯一解 定理4 1n 1个节点的求积公式为插值型求积公式的充要条件是公式至少具有n次代数精度 证 充分性设n 1个节点的求积公式为插值型求积公式 求积系数为又当f x 为不高于n次的多项式时 f x P x 其余项R f 0 因而这时求积公式至少具有n次代数精度 定理4 1n 1个节点的求积公式为插值型求积公式的充要条件是公式至少具有n次代数精度 必要性若求积公式至少具有n次代数精度 则对n次多项式 精确成立 即 而 例4 1设积分区间 a
7、 b 为 0 2 取时时 分别用梯形和辛卜生公式 计算其积分结果并与准确值进行比较解 梯形公式和辛卜生的计算结果与准确值比较如下表所示 f x 1xx2x3x4ex准确值222 6746 406 389梯形公式计算值2248168 389辛卜生公式计算值222 6746 676 421 从表中可以看出 当f x 是时 辛卜生公式比梯形公式更精确 一般说来 代数精度越高 求积公式越精确 梯形公式和中矩形公式具有1次代数精度 辛卜生公式有3次代数精度 下面以梯形公式为例进行验证 取f x 1时 两端相等 取f x x时 所以梯形公式只有1次代数精度 两端相等 例4 2试确定一个至少具有2次代数精度
8、的公式 解 要使公式具有2次代数精度 则对f x 1 x x2求积公式准确成立 即得如下方程组 解之得 所求公式为 例4 3试确定求积系数A B C使具有最高的代数精度解 分别取f x 1 x x2使求积公式准确成立 即得如下方程组 由于n 1节点的插值求积公式至少有n次代数精度 所以构造求积公式后应该验算所构造求积公式的代数精度 例如插值求积公式 有三个节点至少有2次代数精度 是否有3次代数精度呢 将f x x3代入公式两端 左端和右端都等于 b4 a4 4 公式两端严格相等 再将f x x4代入公式两端 两端不相等 所以该求积公式具有3次代数精度 的代数精度可以验证 对于f x 1 x时公
9、式两端相等 再将f x x2代入公式左端 例4 4考察求积公式 两端不相等 所以该求积公式具有1次代数精度 三个节点不一定具有2次代数精度 因为不是插值型的 右端 例4 5给定求积公式如下 试证此求积公式是插值型的求积公式 证 设 则以这三点为插值节点的Lagrange插值基函数为 由插值型求积公式的定义知 所给的求积公式是插值型求积公式 插值型求积公式为 例4 6求证 不是插值型的 证明 设x0 1 x1 0 x2 1 A0 1 2 A1 1 A2 1 2则以这三点为插值节点的Lagrange插值基函数为 例4 7给定求积公式 试确定求积系数A 1 A0 A1 使其有尽可能高的代数精度 并指
10、出其代数精度 解 令求积公式对f x 1 x x2准确成立 则有 例4 7给定求积公式 试确定求积系数A 1 A0 A1 使其有尽可能高的代数精度 并指出其代数精度 解之得 其代数精度至少为2 将f x x3代入求积公式两端相等 而将将f x x4代入求积公式两端不相等 所以其代数精度为3次 例4 8确定求积公式 使其具有尽可能高的代数精度 解 不妨设a 0 b h b a h 设所求公式的代数精度为2 则当f x 1 x x2时公式变成等式 即 例4 8确定求积公式 使其具有尽可能高的代数精度 解 不妨设a 0 b h b a h 设所求公式的代数精度为2 则当f x 1 x x2时公式变成
11、等式 即 其中h b a 令f x x3代入上式 两端不等 说明求积公式只有2次代数精度 解之得 构造插值求积公式有如下特点 复杂函数f x 的积分转化为计算多项式的积分求积系数Ak只与积分区间及节点xk有关 而与被积函数f x 无关 可以不管f x 如何 预先算出Ak的值n 1个节点的插值求积公式至少具有n次代数精度求积系数之和可用此检验计算求积系数的正确性 例4 9求证当节点为n 1个时 插值求积系数之和为 1 在积分区间 a b 上选取节点xk 2 求出f xk 及利用或解关于Ak的线性方程组求出Ak 这样就得到了 3 利用f x xn 验算代数精度 构造插值求积公式的步骤 例4 9对构
12、造一个至少有3次代数精度的求积公式 解 3次代数精度需4个节点 在 0 3 上取0 1 2 3四个节点构造求积公式 确定求积系数Ak k 0 1 2 3 利用求积系数公式 因为求积公式有4个节点 所以至少具有3次代数精度 只需将f x x4代入来验证其代数精度 将f x x4代入两端不相等 所以只有3次代数精度 4 3牛顿 柯特斯 Newton Cotes 求积公式在插值求积公式 中 当所取节点是等距时称为牛顿 柯特斯公式其中插值多项式求积系数 这里是插值基函数 即有 将积分区间 a b 划分为n等分 步长求积节点为为了计算系数Ak 由于 所以 作变量代换当时 有 于是可得 代入插值求积公式
13、4 4 有 称为牛顿 柯特斯求积公式 Ck称为柯特斯系数 容易验证 显然 Ck是不依赖于积分区间 a b 以及被积函数f x 的常数 要给出n 就可以算出柯特斯系数 譬如当n 1时 当n 2时 下表给出了n从1 8的柯特斯系数 我们看到 当n 8时 柯特斯系数有正有负 这时稳定性得不到保证 因此 实际计算不用高阶的牛顿 柯特斯公式 4 3 1偶阶求积公式的代数精度作为插值型的求积公式 n阶的牛顿 柯特斯公式至少具有n次的代数精度 定理1 实际的代数精度能否进一步提高呢 先看辛甫生公式 它是二阶牛顿 柯特斯公式 因此至少具有二次代数精度 进一步用f x x3进行检验 按辛甫生公式计算得 另一方面
14、直接求积得 这时有S I即辛甫生公式对次数不超过三次的多项式均能准确成立 又容易验证它对f x x4通常是不准确的 因此 辛甫生公式实际上具有三次代数精度 一般地 我们可以证明下述论断 定理当阶n为偶数时 牛顿 柯特斯公式至少有n 1次代数精度 证明只要验证 当n为偶数时 牛顿 柯特斯公式对f x x n 1 的余项为零 由于这里从而有 引进变换x a th 并注意到xj a jh 有 若n为偶数 则n 2为整数 再令t u n 2进一步有 是个奇函数 证毕 据此可以断定 因为被积函数 4 4几个低阶求积公式在牛顿 柯特斯求积公式中n 1 2 4时 就分别得到下面的梯形公式 辛卜生公式和柯特斯
15、公式 1 梯形公式当n 1时 牛顿 柯特斯公式就是梯形公式 定理4 2 梯形公式的误差 设f x 在 a b 上具有连续的二阶导数 则梯形公式的误差 余项 为 证 由插值型求积公式的余项 由于 x a x b 在 a b 中不变号 在 a b 上连续 根据高等数学中的积分中值定理 在 a b 上存在一点 使 因此 可知梯形公式的误差为 其中 2 辛卜生公式当n 2时 牛顿 柯特斯公式就是辛卜生公式 或称抛物线公式 定理4 3 辛卜生公式的误差 设在 a b 上具有连续的四阶导数 则辛卜生求积公式的误差为 3 柯特斯公式 当n 4时 牛顿 柯特斯公式为 定理4 4 柯特斯公式的误差 设在 a b
16、 上具有连续的6阶导数 则柯特斯求积公式的误差为 例4 11分别用梯形公式 辛卜生公式和柯特斯公式计算定积分的近似值 计算结果取5位有效数字 1 用梯形公式计算 2 用辛卜生公式 3 用柯特斯公式计算 系数为 积分的准确值为 可见 三个求积公式的精度逐渐提高 例4 12用辛卜生公式和柯特斯公式计算定积分 的近似值 并估计其误差 计算结果取5位小数 解 辛卜生公式 由于 知其误差为 由辛卜生公式余项 例4 12用辛卜生公式和柯特斯公式计算定积分 的近似值 并估计其误差 计算结果取5位小数 解 柯特斯公式 知其误差为 4 5复化求积公式由梯形 辛卜生和柯特斯求积公式余项可知 随着求积节点数的增多
17、对应公式的精度也会相应提高 但由于n 8时的牛顿 柯特斯求积公式开始出现负值的柯特斯系数 根据误差理论的分析研究 当积分公式出现负系数时 可能导致舍入误差增大 并且往往难以估计 因此不能用增加求积节点数的方法来提高计算精度 在实际应用中 通常将积分区间分成若干个小区间 在每个小区间上采用低阶求积公式 然后把所有小区间上的计算结果加起来得到整个区间上的求积公式 这就是复化求积公式的基本思想 常用的复化求积公式有复化梯形公式和复化辛卜生公式 4 5 1复化梯形公式及其误差将积分区间 a b 划分为n等分 步长h b a n求积节点为xk a kh k 0 1 n 在每个小区间上应用梯形公式 xk
18、xk 1 k 0 1 n 求出积分值Ik 然后将它们累加求和 用作为所求积分I的近似值 记 4 5 4 5 式称为复化梯形公式 当f x 在 a b 上有连续的二阶导数 在子区间上梯形公式的余项已知为 在 a b 上的余项 设在 a b 上连续 根据连续函数的介值定理知 存在 使 因此 余项 4 5 3复化辛卜生公式及其误差将积分区间 a b 划分为n等分 记子区间的中点为在每个小区间上应用辛卜生公式 则有 记 4 6 称为复化辛卜生公式 类似于复化梯形公式余项的讨论 复化辛卜生公式 4 6 的求积余项为 如果把每个子区间四等分 内分点依次记 同理可得复化柯特斯公式 求积余项为 复化求积公式的
19、余项表明 只要被积函数f x 所涉及的各阶导数在 a b 上连续 那么复化梯形公式 复化辛卜生公式与复化柯特斯公式所得近似值的余项和步长的关系依次为 因此当h 0 即n 时 都收敛于积分真值 且收敛速度一个比一个快 例4 13依次用n 8的复化梯形公式 n 4的复化辛卜生公式计算定积分 解 首先计算出所需各节点的函数值 n 8时 由复化梯形公式 4 5 可得如下计算公式 由复化辛卜生公式 4 6 可得如下计算公式 积分准确值I 0 9460831 这两种方法都需要提供9个点上的函数值 计算量基本相同 然而精度却差别较大 同积分的准确值 是指每一位数字都是有效数字的积分值 比较 复化梯形法只有两
20、位有效数字 T8 0 9456909 而复化辛卜生法却有六位有效数字 例4 14用复化梯形公式计算定积分才能使误差不超过 解 取 则 又区间长度b a 1 对复化梯形公式有余项 即 n 212 85 取n 213 即将区间 0 1 分为213等份时 用复化梯形公式计算误差不超过 问区间 0 1 应分多少等份 定义2如果一种复化求积公式 当时成立渐近关系式 定数 则称求积公式是p阶收敛的 在这种意义下 复化的梯形法 辛甫生法和柯特斯法分别具有二阶 四阶和六阶收敛性 而当h很小时 对于复化的梯形法 辛甫生法和柯特斯法分别有下列误差估计式 4 14 4 15 4 16 由此可见 若将步长h减半 即等
21、分数n加倍 则梯形法 辛甫生法与柯特斯法的误差分别减至原有误差的1 4 1 16与1 64 4 6龙贝格 Romberg 求积法复化求积方法对于提高计算精度是行之有效的方法 但复化公式的一个主要缺点在于要先估计出步长 若步长太大 则难以保证计算精度 若步长太小 则计算量太大 并且积累误差也会增大 在实际计算中通常采用变步长的方法 即把步长逐次分半 直至达到某种精度为止 4 6 1变步长的梯形公式变步长复化求积法的基本思想是在求积过程中 通过对计算结果精度的不断估计 逐步改变步长 逐次分半 直至满足精度要求为止 即按照给定的精度实现步长的自动选取 设将积分区间 a b n等分 即分成n个子区间
22、一共有n 1个节点 即x a kh k 0 1 n 步长 对于某个子区间 利用梯形公式计算积分近似值有 对整个区间 a b 有 将子区间再二等份 取其中点作新节点 此时区间数增加了一倍为2n 对某个子区间 利用复化梯形公式计算其积分近似值 对整个区间 a b 有 比较和有 4 7 4 7 式称为变步长梯形公式 当把积分区间分成n等份 用复化梯形公式计算积分I的近似值时 截断误差为 若把区间再分半为2n等份 计算出定积分的近似值 则截断误差为 当在区间 a b 上变化不大时 有 所以 可见 当步长二分后误差将减至1 4 将上式移项整理 可得验后误差估计式 上式说明 只要二等份前后两个积分值和相当
23、接近 就可以保证计算结果的误差很小 使接近于积分值I 4 6 2变步长的梯形求积算法实现 1 变步长的梯形求积法的计算步骤 变步长梯形求积法 它是以梯形求积公式为基础 逐步减少步长 按如下递推公式求二分后的梯形值 其中Tn和T2n分别代表二等分前后的积分值 如果 为给定的误差限 则T2n作为积分的近似值 否则继续进行二等分 即 转 再计算 直到满足所要求的精度为止 最终取二分后的积分值T2n作为所求的结果 例4 15用变步长梯形求积法计算定积分 所以有 然后将区间二等份 由于 故有 进一步二分求积区间 并计算新分点上的函数值 解 先对整个区间 0 1 用梯形公式 对于 有 这样不断二分下去 积
24、分的准确值为0 9460831 用变步长二分10次可得此结果 4 6 3龙贝格求积公式变步长梯形求积法算法简单 但精度较差 收敛速度较慢 但可以利用梯形法算法简单的优点 形成一个新算法 这就是龙贝格求积公式 龙贝格公式又称逐次分半加速法 根据积分区间分成n等份和2n等份时的误差估计式 4 8 可得 所以积分值的误差大致等于 如果用对进行修正时 与之和比更接近积分真值 所以可以将看成是对误差的一种补偿 因此可得到具有更好效果的式子 4 9 考察与n等份辛卜生公式之间的关系 将复化梯形公式 梯形变步长公式 代入 4 9 表达式得 故 这就是说 用梯形法二分前后两个积分值和作线性组合 结果却得到复化
25、辛卜生公式计算得到的积分值 再考察辛卜生法 其截断误差与成正比 因此 如果将步长折半 则误差减至 即有 由此可得 可以验证 上式右端的值其实等于Cn 就是说 用辛卜生公式二等份前后的两个积分值Sn和S2n作线性组合后 可得到柯特斯公式求得的积分值Cn 即有 4 11 用同样的方法 根据柯特斯公式的误差公式 可进一步导出龙贝格公式 4 12 在变步长的过程中运用 4 10 4 11 和 4 12 就能将粗糙的梯形值Tn逐步加工成精度较高的辛卜生值Sn 柯特斯值Cn和龙贝格值Rn或者说 将收敛缓慢的梯形值序列Tn加工成收敛迅速的龙贝格值序列Rn 这种加速方法称为龙贝格算法 龙贝格公式 见下表所示
26、4 6 4龙贝格求积法算法实现 1 龙贝格求积法计算步骤用梯形公式计算积分近似值按变步长梯形公式计算积分近似值将区间逐次分半 令区间长度 计算 按加速公式求加速值 梯形加速公式 辛卜生加速公式 龙贝格求积公式 精度控制 直到相邻两次积分值 其中 为允许的误差限 则终止计算并取Rn作为积分的近似值 否则将区间再对分 重复 的计算 直到满足精度要求为止 例4 16用龙贝格算法计算定积分要求相邻两次龙贝格值的偏差不超过 解 由题意 由于 于是有 4 8数值微分在微分学中 求函数f x 的导数通常是可以求得的 但有的比f x 复杂的多 另外 有时f x 仅由表格形式给出 则求也不容易 根据函数在若干个
27、点处的函数值去求该函数的导数的近似值称为数值微分 求数值导数也是实际问题经常遇到的 特别当该函数本身未知 但又需要对其求导数时 数值微分方法显得更为重要 4 8 1中点方法最简单的数值微分是用差商近似代替导数 即 同样 也可用向后差商近似代替导数 即 或中心差商的方法 即 可以看出中心差商是向前差商和向后差商的算术平均值 上述三种方法的截断误差分别为 和 如右图所示 前述三种导数的近似值分别表示弦线AB AC和BC的斜率 将这三条通过A点的弦的斜率与切线 利用中点公式计算导数 首先必须选取合适的步长 为此需要进行误差分析 分别将在x a处泰勒展开 有 代入 4 17 得 4 18 由此可知 从
28、截断误差的角度来看 步长越小 计算结果越准确 但从舍入误差角度 h越小 与越接近 直接相减会造成有效数字的严重损失 就舍入误差而言 步长是不宜太小的 怎样选择最佳步长 使截断误差与舍入误差之和最小呢 由式 4 18 知 当h适当小时 因而 对上式两边同乘以 并移项后得 由此可以看出 只要当二分前后的2个近似值G h 和很接近 就可以保证的截断误差很小 大致等于 所以比较二分前后所得的G h 和 若 则为所取的合适的步长 且 否则将步长再二等分 继续进行计算 4 8 2插值型求导公式函数f x 的定积分可以用插值多项式P x 的定积分来近似计算 同样 函数f x 的导数也可以用插值多项式P x
29、的导数来近似代替 即 4 19 这样建立的数值微分公式统称为插值型求导公式应当指出的是即使P x 与f x 处处相差不多 但与在某些点仍然可能出入很大 因而在使用插值求导公式时要注意误差的分析 由插值余项公式 得求导公式 4 19 的余项为 式中 在这一余项公式中 由于 和x是未知函数 因此无法对它的第二项作出估计 但在插值节点xk处 由于上式右端的第二项因式等于零 因而在插值节点处的导数余项为 4 20 下面给出实用的两点公式和三点公式 1 两点公式设已给出两个节点上的函数值 作线性插值 对上式两端求导 记 则有 注意到 于是有下列求导的两点公式 而利用余项公式知 带余项的两点公式是 2 三
30、点公式设已给出三个节点上的函数值 作二次插值 得 令 上式可表示为 两端对t求导 有 上式分别取t 0 1 2 得到三种三点公式 而带余项的三点公式如下 式中 截断误差是 用插值多项式P x 作为f x 的近似函数 还可以建立高阶数值求导公式 k 0 1 n 本章小结 本章介绍了积分和微分的数值计算方法 其基本原理主要是逼近论 即设法构造某个简单函数P x 近似表示f x 然后对P x 求积或求导得到f x 的积分或导数的近似值 基于插值原理 推导了数值积分和数值微分的基本公式 插值型求积公式介绍了牛顿 柯特斯公式和高斯公式两类 前者取等距节点 算法简单而容易编制程序 但是 由于在n 8时出现了负系数 从而影响稳定性和收敛性 因此实用的只是低阶公式 解决长区间与低阶公式的矛盾是使用复化求积公式 龙贝格算法是在区间逐次分半过程中 对用梯形法所获得的近似值进行多级 加工 从而获得高精度的积分近似值的一种方法 它具有自动选取步长且精度高 计算量小的特点 便于在计算机上使用 是数值积分中较好的方法 必须熟练地掌握 建立在代数精度概念上的待定系数法也是数值积分中的一般方法 按待定系数法确定的数值积分公式没有误差估计式 只能从代数精度出发 估计其精确程度 如果要求函数在某点的导数 而该函数仅仅是列表函数时 可使用插值求导公式去近似计算函数的导数 Thankyouverymuch
