数值分析第4章数值积分与数值微分1详解

1 第第 4 章章 数值积分与数值微分数值积分与数值微分 1 数值积分的基本概念数值积分的基本概念 实际问题当中常常需要计算定积分 在微积分中 我们熟知 牛顿 莱布尼兹公式是计算定积 分的一种有效工具 在理论和实际计算上有很大作用 对定积分 若在区间 b a If x dx f x 上连续 且的原函数为 则可计算定积分 a b f x F x b a f x dxF bF a 似乎问题已经解决 其实不然 如 1 是由测量或数值计算给出数据表时 Newton Leibnitz 公式无法应用 f x 2 许多形式上很简单的函数 例如 2 322 sin1 1 sin cos ln x x f xxxxe xx 等等 它们的原函数不能用初等函数的有限形式表示 3 即使有些被积函数的原函数能通过初等函数的有限形式表示 但应用牛顿 莱布尼兹公式 计算 仍涉及大量的数值计算 还不如应用数值积分的方法来得方便 既节省工作量 又满足精度 的要求 例如下列积分 2 4 2 11211 lnarc 21 arc 21 14 2212 2 xx dxtgxtgxC xxx 对于上述这些情况 都要求建立定积分的近似计算方法 数值积分法 1 1 数值求积分的基本思想数值求积分的基本思想 根据以上所述 数值求积公式应该避免用原函数表示 而由被积函数的值决定 由积分中值定 理 对 存在 有 f xC a b a b b a f x dxba f 表明 定积分所表示的曲边梯形的面积等于底为而高为的矩形面积 图 4 1 问题在ba f 于点的具体位置一般是不知道的 因而难以准确算出 我们将称为区间上的平 f f a b 均高度 这样 只要对平均高度提供一种算法 相应地便获得一种数值求积分方法 f 如果我们用两端的算术平均作为平均高度的近似值 这样导出的求积公式 f 4 1 2 ba Tf af b 便是我们所熟悉的梯形公式梯形公式 图 4 2 而如果改用区间中点的 高度 近似地取代 2 ab c f c 平均高度 则可导出所谓中矩形公式中矩形公式 简称矩形公式矩形公式 f 4 2 2 ab Rba f 2 更一般地 我们可以在区间上适当选取某些节点 然后用加权平均得到平均高 a b k x k f x 度的近似值 这样构造出的求积公式具有下列形式 f a bab yf x yf x f yy xx 图 4 1 图 4 2 4 3 0 n b kk a k f x dxA f x 式中称为求积节点求积节点 成为求积系数求积系数 亦称伴随节点的权权 权仅仅与节点的选取有关 k x k A k x k A k x 而不依赖于被积函数的具体形式 f x 这类由积分区间上的某些点上处的函数值的线性组合作为定积分的近似值的求积公式通常称为 机械求积公式 机械求积公式 它避免了 Newton Leibnitz 公式寻求原函数的困难 对于求积公式 4 3 关键在于 确定节点和相应的系数 k x k A 1 2 代数精度的概念代数精度的概念 由 Weierstrass 定理可知 对闭区间上任意的连续函数 都可用多项式一致逼近 一般说来 多 项式的次数越高 逼近程度越好 这样 如果求积公式对阶多项式精确成立 那么求积公式的m 误差仅来源于阶多项式对连续函数的逼近误差 因此自然有如下的定义m 定义定义 4 1 如果某个求积公式对于次数不超过的多项式均准确地成立 但对于次多项式m1m 就不准确成立 则称该求积公式具有次代数精度 次代数精度 m 例例 1 判断求积公式 1 1 1 5 0 6 8 0 5 0 6 9 f x dxfff 的代数精度 解 记 1 1 1 5 0 6 8 0 5 0 6 9 I ff x dxI ffff 因为 1 1 1 1 2 1 585 2 9 IdxI 1 1 1 50 6805 0 6 0 9 I xxdxI x 0 3 1 222 1 12 50 68050 6 93 I xx dxI x 2 3 1 33333 1 1 5 0 6 05 0 6 0 9 I xx dxI x 0 1 444 1 12 50 36050 36 95 I xx dxI x 2 5 1 55555 1 1 5 0 6 05 0 6 0 9 I xx dxI x 0 1 66633 1 12 5 0 6 05 0 6 0 24 97 I xx dxI x 2 7 所以求积公式具有 5 次代数精度 1 3 插值型的求积公式插值型的求积公式 最直接自然的一种想法是用在上的插值多项式代替 由于代数多项式 f x a b n x f x 的原函数是容易求出的 我们以在上的积分值作为所求积分的近似值 即 n x a b I f b n a I fx dx 这样得到的求积分公式称为插值型求积公式 通常采用 Lagrange 插值 设上有个互异节点 的次 Lagrange 插值多项式为 a b1n 01 n xxx f xn 0 n nkk k Lxlx f x 其中 插值型求积公式为 0 n i k j ki j k xx lx xx 4 4 0 n b nkk a k I fLx dxA f x 其中 可看出 仅由积分区间与插值节点确定 与 0 1 b kk a Alx dx kn k A a b k x 被积函数的形式无关 求积公式 4 4 的截断误差为 f x 4 5 1 1 1 n bbb nnn aaa f Rff x dxLx dxx dx n 定义定义 4 2 求积公式 0 n b kk a k f x dxA f x 如其系数 则称此求积公式为插值型求积公式 b kk a Alx dx 定理定理 4 1 形如 4 3 的求积公式至少有次代数精度的充分必要条件是插值型的 n 证明证明 如果求积公式 4 3 是插值型的 由公式 4 5 可知 对于次数不超过的多项式 其n f x 余项等于零 因而这时求积公式至少具有次代数精度 R fn 反之 如果求积公式 4 3 至少具有次代数精度 那么对于插值基函数应准确成立 并n k lx 注意到 即有 kjjk lx 4 0 n b kj kjk a j lx dxA lxA 所以求积公式 4 3 是插值型的 1 4 求积公式的收敛性与稳定性求积公式的收敛性与稳定性 定义定义 4 3 在求积公式 4 3 中 若 0 0 lim n b kk an k h A f xf x dx 其中 则称求积公式 4 3 是收敛的 1 1 max ii i n hxx 实际使用任何求积公式时 除截断误差外 还有舍入误差 因此我们必须研究其数值稳定性 在求积公式 4 3 中 由于计算可能产生误差 实际得到 即 记 k f x k k f kkk f xf 00 nn nkknkk kk IfA f xIfA f 如果对任给正数 只要误差充分小就有0 k 4 6 0 n nnkkk k IfIfAf xf 它表明求积公式 4 3 计算是稳定的 由此给出 定义定义 4 4 对任给 若存在 只要就有 4 6 成立 0 0 0 1 kk f xfkn 则称求积公式 4 3 是稳定的 定理定理 4 2 若求积公式 4 3 中系数 则此求积公式是稳定的 若有正0 0 1 k Akn k A 有负 计算可能不稳定 证明证明 对任给 若取 对都有 则有0 ba 0 1 kn kk f xf 000 nnn nnkkkkkkk kkk IfIfAf xfAf xfA 注意对任何代数精度的求积公式均有0 0 1 1 n b kn a k AIdxba 可见时 有0 k A 00 nn nnkk kk IfIfAAba 由定义 4 4 可知求积公式 4 3 是稳定的 若有正有负时 假设 且 有 k A 0 kkk Af xf kk f xf 00 00 nn nnkkkkkk kk nn kk kk IfIfAf xfAf xf AAba 5 它表明初始数据的误差可能会引起计算结果误差的增大 即计算可能不稳定 2 Newton Cotes 公式公式 2 1 Cotes 系数系数 被积函数在积分区间内变化平缓 可用等距节点插值公式近似 将积分区间划分为等 a bn 分 步长 等距节点 此时求积公式 4 4 中的积分系数可得到 ba h n 0 1 k xakh kn 简化 00 nn bbb j kk aaa jj kj j kj k xx xajh Alx dxdxdx xxkj h 作变换 则有xath 000 000 1 1 n kn knnn nnn k jjj j kj kj k tj hhba Ahdttj dttj dt kj hk nkk nkn 令 0 0 1 n kn n n k j j k Ctj dt k nkn 则 求积公式 4 4 可简化为 n kk Aba C 4 7 0 n n kk k I fbaCf x 称为阶 Newton Cotes 公式 简记为 N C 公式 称为 Cotes 系数 n n k C 由的表达式可看出 它不但与被积函数无关 而且与积分区间也无关 因此可将 Cotes 系 n k C 数事先列成表格供查用 见表 4 1 N C 公式的截断误差为 4 8 1 2 1 0 00 1 1 nnnn bn n nj a jj fh Rfxx dxftj dt nn 时1n 4 9 11 222 ba I fbaf af bf af b 为梯形公式 时2n 4 10 141 4 662662 abbaab I fbaf aff bf aff b 为辛普生公式 时4n 6 4 11 3 7 32 12 32 7 90424 babababa I ff af aff af b 为 Cotes 公式 表 4 1 n n k C 1 1 2 1 2 2 1 6 4 6 1 6 3 1 8 3 8 3 8 1 8 4 7 90 16 45 2 15 16 45 7 90 5 19 288 25 96 25 144 25 144 25 96 19 288 6 41 840 9 35 9 280 34 105 9 280 9 35 41 840 7 751 17280 3577 17280 1323 17280 2989 17280 2989 17280 1323 17280 3577 17280 751 17280 8 989 28350 5888 28350 928 28350 10496 28350 4540 28350 10496 28350 928 28350 5888 28350 989 28350 从表 4 1 可看出 当时出现了负系数 由定理 4 2 可知 实际计算中将使舍入误差增大 8n 并且往往难以估计 从而 Newton Cotes 公式的收敛性和稳定性得不到保证 因此实际计算中不用 高阶 Newton Cotes 公式 2 2 偶阶求积公式的代数精度偶阶求积公式的代数精度 作为插值型的求积公式 阶的牛顿 柯特斯公式至少具有次的代数精度 求积公式的代数nn 精度能否进一步提高呢 定理定理 4 3 当阶为偶数时 公式 4 7 至少具有次代数精度 NC 1n 证明证明 我们只要验证 当为偶数时 N C 公式对的余项为零 按余项公式 4 8 n 1 n f xx 由于这里 从而 1 1 n fxn 0 n b nj a j Rfxx dx 引进变换 并注意到 有xath j xajh 2 0 0 n n n n j Rfhtj dt 当为偶数 则为整数 再令 进一步有n 2 n 2 n tu 2 22 22 02 22 0 2 nn nn nn nnn jjn n Rfhuj duhuj du 因为被积函数为奇函数 7 2 3 几种低阶求积公式的余项几种低阶求积公式的余项 梯形求积公式的余项为 2 b T a f RITxa xb dx 由于在上不变号 利用积分中值定理有 xa xb a b 4 12 3 212 b T a ff Rxa xb dxbaa b Simpson 公式的余项为 4 6 b S a ba RISf x dxf af cf b 这里 构造次数不超过 3 的多项式 使满足 2 ab c H x H af a H cf c H cfc H bf b 由于 Simpson 公式具有三次代数精度 它对于这样构造的三次式是准确的 即 H x 4 6 b a ba H x dxH aH cH b 所以 b S a Rf xH xdx 由第二章的例 6 可知 4 2 1 4 f xH xfxa xcxb 因在上保号 应用积分中值定理有 2 xa xcxb a b 4 13 4 2 4 4 1 4 1802 b s a Rfxa xcxb dx ba ba fa b 3 复化求积公式复化求积公式 前面导出的误差估计式表明 用 N C 公式计算积分近似值时 步长越小 截断误差越小 但 缩小步长等于增加节点数 亦即提高插值多项式的次数 Runge 现象表明 这样并不一定能提高精 度 理论上已经证明 当时 N C 公式所求得的近似值不一定收敛于积分的准确值 而且n 随着的增大 N C 公式是不稳定的 因此 实际中不采用高阶 N C 公式 为提高计算精度 可考n 虑对被积函数用分段低次多项式插值 由此导出复化求积公式 3 1 复化梯形公式复化梯形公式 8 将区间划分为等分 分点 在每个区间 a bn 0 1 k ba xakh hkn n 1 kk xx 上采用梯形公式 则得 0 1 1 kn 4 14 1 11 1 00 2 k k nn bx kkn ax kk h If x dxf x dxf xf xRf 记 4 15 11 1 01 2 22 nn nkkk kk hh Tf xf xf af xf b 称为复化梯形公式复化梯形公式 其余项为 32 11 1 00 1 1212 nn nnkkkkk kk hba h RfITffxx n 由于 且 2 f xCa b 1 0101 0 1 min max n kkk k nk n k fff n 所以存在使 a b 1 0 1 n k k ff n 于是复化梯形公式余项为 4 16 2 12 n ba Rfh f 从式 4 16 可以看出 余项误差是阶 所以当 有 2 h 2 f xCa b lim b n an Tf x dx 即复化梯形公式是收敛的 事实上只要 则可得收敛性 因为由 4 15 得 f xC a b 1 01 1 2 nn b nkk a kk baba Tf xf xf x dx n nn 所以复化梯形公式 4 15 收敛 此外 的求积系数为正 由定理 4 2 知复化梯形公式是稳定的 n T 3 2 复化辛普森公式复化辛普森公式 将区间划分为等分 在每个区间上采用辛普森公式 记则 a bn 1 kk xx 1 2 1 2 kk xxh 得 4 17 1 11 1 21 00 4 6 k k nn bx kkkn ax kk h If x dxf x dxf xf xf xRf 记 4 18 11 1 2 01 4 2 6 nn nkk kk h Sf af xf xf b 称为复化辛普森求积公式复化辛普森求积公式 其余项由 4 13 得 9 4 1 4 1 0 1 180 2 n nnkkkk k h RfIShfxx 于是当时 与复化梯形公式相似有 4 f xCa b 4 19 4 4 1802 nn bah RfISfa b 可以看出误差阶是 收敛性是显然的 事实上 只要 则有 4 h f xC a b 11 1 2 001 1 4 6 nnn nkkk kkk b a bababa Sf xf xf x nnn f x dx n 此外 由于中求积系数均为正数 故知复化辛普森求积公式计算稳定 n S 例例 2 根据函数表 4 1 表 4 1 k k x sin k k k x f x x k k x sin k k k x f x x 00155 80 9361556 11 80 997397863 40 9088516 21 40 989615877 80 8771925 33 80 9767267810 8414709 41 20 9588510 用复化梯形公式和复化辛普森公式计算的近似值 并估计误差 1 0 sin x Idx x 解 由复化梯形公式 7 1 1 0 1 2 0 945691 168 k k Ifff 由复化辛普森公式 34 11 121 0 1 2 4 0 946084 1648 kk kk Iffff 与准确值比较 显然用复化 Simpson 公式计算精度较高 0 9460831I 为了利用余项公式估计误差 要求的高阶导数 由于 sin x f x x 1 0 sin cos x f xxt dt x 所以有 11 00 cos cos 2 k kk k dk fxxt dttxtdt dx 于是 11 0001 1 max cos 21 kkk x k fxtxtdtt dt k 10 由复化梯形误差公式 4 16 得 2 2 88 01 111 max 0 000434 1212 83 x h RfITfx 由复化辛普森误差公式 4 19 得 4 6 44 111 0 271 10 180 85 RfIS 例例 3 若用复化求积分公式计算积分 1 0 x Ie dx 的近似值 要求计算结果有四位有效数字 应取多大 n 解 因为当时 有01x 1 0 31 x ee 于是 1 0 0 31 x e dx 要求计算结果有四位有效数字 即要求误差不超过 又因为 4 1 10 2 1 0 1 kx fxex 由式 4 16 得 2 24 11 10 12122 T h Rhf 即 开方得 因此若用复化梯形公式求积分 应等于 41 才能达到精度 4 1 10 6 n 40 8n n 若用复化 Simpson 公式 由式 4 19 44 4 4 4 1111 10 180 2180 16180 162 S hh Rf n 即得 故应取 1 62n 2n 4 龙贝格求积公式龙贝格求积公式 4 1 梯形公式的逐次分半算法梯形公式的逐次分半算法 如前所述 复化求积公式的截断误差随着步长的缩小而减少 而且如果被积函数的高阶导数容 易计算和估计时 由给定的精度可以预先确定步长 不过这样做常常是很困难的 一般不值得推崇 实际计算时 我们总是从某个步长出发计算近似值 若精度不够可将步长逐次分半以提高近似值 直到求得满足精度要求的近似值 设将区间分为等分 共有个分点 如果将求积区间再二分一次 则分点增至 a bn1n 个 我们将二分前后两个积分值联系起来加以考虑 注意到每个子区间经过二分只21n 1 kk xx 11 增加了一个分点 用复化梯形公式求得该子区间上的积分值为 11 2 1 2 kk k xxx 1 21 2 4 kkk h f xf xf x 注意 这里代表二分前的步长 将每个子区间上的积分值相加得 ba h n 11 211 2 00 42 nn nkkk kk hh Tf xf xf x 即 4 20 1 21 2 0 1 22 n nnk k h TTf x 这表明 将步长由缩小为时 等于的一半再加新增加节点处的函数值乘以当前步长 h 2 h 2n T n T 算法算法 4 1 1 输入 a b f x 2 置 0 1 2 ba mhTh f af b 3 置 对0F 1 1 2 2mk 21 FFf akh 4 0 1 2 TThF 5 若 输出 停机 否则 转 3 0 3TT IT 0 1 2 h mmh TT 4 2 李查逊 李查逊 Richardson 外推法 外推法 假设用某种数值方法求量的近似值 一般地 近似值是步长的函数 记为 相应的Ih 1 I h 误差为 4 21 12 112 k ppp k II hhhh 其中是与无关的常数 若用代替 4 21 中的 则 12 1 2 0 ik ippp hh h 得 4 22 12 1122 112 12 k kk ppp k pppppp k IIhhhh hhh 式 4 22 减去式 4 21 乘以 得 1 p 1 3321211 11 23 kk p ppppppppp k IIhII h hhh 取满足 以除上式两端 得 1 1 1 p 4 23 1 32 1 11 23 1 k p ppp k p IhI h Ib hb hb h 其中仍与无关 令 11 2 1 2 3 i ppp i bi h 12 1 1 11 2 1 p p IhI h Ih 由式 4 23 以作为的近似值 其误差至少为 因此收敛于的速度比 2 IhI 2 p O h 2 IhI 快 不断重复以上作法 可以得到一个函数序列 1 I h 4 24 1 1 11 2 3 1 m m p mm m p IhIh Ihm 以近似 误差为 随着的增大 收敛速度越来越快 这就是 m IhI m p m IIhO h m Richardson 外推法 4 3 龙贝格求积公式龙贝格求积公式 由前面知道 复化梯形公式的截断误差为 进一步分析 我们有如下欧拉 麦克劳林 2 O h Euler Maclaurin 公式 定理定理 4 4 设 则有 f xCa b 242 12 k k IT hhhh 其中系数与无关 1 2 k k h 把李查逊外推法与欧拉 麦克劳林公式相结合 可以得到求积公式的外推算法 特别地 在外 推算法式 4 24 中 取 并记 则有 1 2 2 k pk 0 T hT h 4 25 11 4 2 1 2 41 j mm m m h TTh Thm 经过次加速后 余项便取下列形式 1 2 m m 4 26 2 1 2 2 12 mm m ThIhh 上述处理方法通常称为李查逊 李查逊 Richardson 外推加速方法 外推加速方法 为研究 Romberg 求积方法的机器实现 引入记号 以表示二分次后求得的梯形值 且以 0 k Tk 表示序列的次加速值 则依以上递推公式得到 k m T 0 k Tm 1 11 41 1 2 4141 m kkk mmm mm ThTTk 称为龙贝格求积算法龙贝格求积算法 Romberg 公式的计算过程见下表 4 2 表 4 2 kh 0 k T 1 k T 2 k T 3 k T 4 k T 0ba 0 0 T 1 2 ba 1 0 T 0 1 T A 2 4 ba 2 0 T 1 1 T 0 2 T 13 3 8 ba 3 0 T 2 1 T 1 2 T 0 3 T 4 16 ba 4 0 T 3 1 T 2 2 T 1 3 T 0 4 T 算法算法 4 2 1 输入 a b f x 2 置 0 0 1 2 h hba Tf af bk 3 计算 1 2 1 00 1 11 22 k kk i ThThf aih 对1 jk 1 11 4 41 jkjkj jj kj j j TT T 4 输出停机 否则 返回 3 0 0 1kk TT 0 b k a Tf x dx 1 2 h h kk 例例 4 用 Romberg 算法计算积分 1 3 2 0 Ixdx 解 利用逐次分半算法 4 20 和 Romberg 算法 4 25 计算结果见表 4 3 0 0 1 0 1 0 500000 2 Tff 1 0 00 11 0 5 0 426777 22 TTf 2 1 00 113 0 25 0 407018 244 TTff 3 2 00 10 251357 0 401812 228888 TTffff 表 4 3 k 0 k T 1 k T 2 k T 3 k T 4 k T 5 k T 00 10 0 20 30 0 0 0 40 0 0 50 0 0 0 5 高斯求积公式高斯求积公式 14 5 1 一般理论一般理论 等距节点的插值型求积公式 虽然计算简单 使用方便 但是这种节点等距的规定却限制了求 积公式的代数精度 试想如果对节点不加限制 并适当选择求积系数 可能会提高求积公式的精度 Gauss 型求积公式的思想也正如此 亦即在节点数固定时 适当地选取节点与求积系数n k x 使求积分公式具有最高精度 k A 设有个互异节点的机械求积分公式1n 01 n xxx 2 27 0 n b kk a k x f x dxA f x 具有次代数精度 那么有取 式 1 精确成立 即m 0 1 l f xxlm 2 28 0 0 1 n b ll jj a j A xx x dx lm 式 2 构成阶的非线性方程组 且具有个未知数 所以当1m 22n 0 1 kk xAkn 给定后 只要 即时 方程组有解 这表明式个节点的求积 x 122mn 21mn 1n 公式的代数精度可达到 21n 另一方面 对式 1 不管如何选择与 最高精度不可能超过 事实上 对 k x k A21n 任意的互异节点 令 0 n kk x 2222 22101 nnn pxxxxxxxx 有 然而 22 0 0 n knk k A Px 21 0 b n a x Px dx 定义定义 4 如果求积分公式 4 27 具有次代数精度 则称这组节点为 Gauss 点点 相应公21n k x 式 4 27 称为带权的高斯求积公式 高斯求积公式 x 定理定理 5 插值型求积公式的节点是高斯点的充分必要条件是以这些节 01n axxxb 点为零点的多项式 101 nn xxxxxxx 与任何次数不超过的多项式带权正交 即n P x 4 29 1 0 b n a x P xx dx 证明证明 必要性 设 则 因此 如果是高斯点 n P xH 121 nn P xxH 01 n xxx 则式 1 对于精确成立 即有 1 n f xP xx 11 0 0 n b nkknk a k x P xx dxA P xx 故式 4 29 成立 再证充分性 对于 用除 记商为 余式为 即 21 n f xH 1 n x f x P x q x 其中 由式 4 29 可得 1 n f xP xxq x n P x q xH 4 30 bb aa x f x dxx q x dx 由于所给求积公式 4 27 是插值型的 它对于是精确成立的 即 n q xH 0 n b kk a k x f x dxA q x 15 再注意到 知 从而由 4 30 有 1 0 0 1 nk xkn 0 1 kk q xf xkn 0 n bb kk aa k x f x dxx q