数值计算方法之数值积分.ppt

第12章 数值积分与数值微分 设f x 是 a b 上连续可积的实函数 求f x 在 a b 上的数值积分也就是求f x 在 a b 上的定积分的数值解 即使我们能找到f x 的一个原函数F x 的解析形式 并利用牛顿 来布尼兹公式进行计算 在许多情况下这也是非常麻烦的 问题的关键在于 初等函数的原函数不一定还是初等函数 在这种情况下并不能利用牛顿 来布尼兹公式进行计算 本章的主要任务是寻找求数值积分的最有效的计算方法 最后建议的方法是变步长复化辛卜生加速算法 为此 需要经过一个曲折的分析过程 12 1求积公式与代数精度的概念 由定积分的定义可知 连续函数f x 在区间 a b 上的定积分近似值可以表示为 a b 内的一些点x0 x1 xn处的函数值f x0 f x1 f xn 的加权和或线性组合 即 1 其中w0 w1 wn仅与x0 x1 xn有关而与被积函数f x 无关 我们把这样的公式称为求积公式 也称为机械求积公式 1 术语和记号 为了计算f x 在区间 a b 上的定积分近似值 我们通常的做法是 把积分区间 a b 划分为n等分 记h b a n x0 a xk a kh k 0 1 2 n 称x0 x1 xn为 a b 的一个等份分划 假如x0 x1 xn为 a b 的一个等份分划那么求积公式 1 中的w0 w1 wn的选取仅仅只与n有关 从而可以简化对求积公式的研究 结论 只要给出了一个如何确定 1 式中的诸w0 w1 wn的机制 我们就可以得到相应的对任何被积函数都有效的计算定积分方法 2 求积公式的性质 微积分学中我们曾研究过 定积分保持函数的线性关系不变 它的含义是 若f x g x 都是 a b 上的可积函数 则对任意实数u v 我们有u f x v g x 也是 a b 上的可积函数 而且不难验证 求积公式也保持函数的线性关系不变 即 3 几种常见的求积公式 在后面的讨论中 我们将经常用到下面一些非常简单的求积公式 他们是中点公式 梯形公式和辛卜生公式 中点公式我是我们课程中强调的一个名词 与求积公式 1 对比分析 可以认为它是这样一种机制 把积分区间分为2等分 取w0 w2 0 w1 1所形成的求积公式 从几何上看 它实际上是取区间中点的函数值与区间长的积作为定积分值 类似于用中位线乘以高来计算梯形的面积 4 截断误差 在求积公式中 我们使用的是近似等号 这是因为 对于一般的被积函数来说 利用这些公式计算所得的结果除了舍入误差外 还有截断误差 因为定积分是用极限来定义的 有时为了进行误差分析 我们可以把上面的 1 式写成 1 其中R f 表示的就是截断误差 考察前面给出的三个求积公式 如果被积函数是线性函数 那么利用中点公式或梯形公式所得到的结果就是准确值 否则一般不是 对于一般的非线性函数 感觉上辛卜生公式更好一些 为了刻划求积公式对一般的被积函数的精确度 我们引进代数精度的概念 5 代数精度的概念 定义 一个求积公式如果对所有的次数不超过m的多项式严格相等 而对某些m 1次多项式不相等 则称该公式具有代数精度m 或该公式的代数精度为m 利用求积公式的线性性 我们不难证明下面的结论 定理 如果求积公式对1 x xm严格相等 而对xm 1不相等 则该公式的代数精度为m 作为课外练习 鼓励大家给出完整证明 6 基本结论 我们可以利用上面的定理所给出的方法证明辛卜生公式的代数精度是3 而中点公式和梯形公式的代数精度是1 现在我们可以对这三个公式作一个简单的评价 中点公式和梯形公式的代数精度虽然都是1 但中点公式只计算一个点的函数值 而梯形公式却要计算两个点处的函数值 所以中点公式优于梯形公式 与梯形公式相比 辛卜生公式只多计算一个点的函数值 但代数精度却增加到3 显然辛卜生公式更为优越 10 2牛顿 柯特斯求积公式 牛顿 柯特斯求积公式就是利用Lagrange插值多项式导出的求积公式 把一般的函数的积分转化为相应的插值多项式函数的积分也是我们学习插值法的基本目的之一 1 利用插值多项式近似替代被积函数 设f x 为被积函数 a b 为积分区间 x0 x1 xn为 a b 内的n 1个互异的点 记Ln x 为相应的拉格朗日插值多项式 那么我们有f x Ln x Rn x 2 利用插值多项式导出求积公式 提示 在上面给出的公式中 由于诸lk x 都是多项式函数 所以诸wk都可以精确地计算出来 从而可以得到一般性的求积公式 2 利用插值多项式导出求积公式 注释 回顾上一章关于多项式插值的结论 由于任意次数不超过n的多项式与它的任意n 1个基点的插值多项式恒等 再由求积公式的代数精度的定义 我们立即得到 由n 1个基点的拉格朗日插值多项式所形成的求积公式的代数精度至少式n 为此 我们上面的wk改写为w n k k 0 1 n 3 牛顿 柯特斯求积公式 牛顿 柯特斯求积公式就是利用等距基点的拉格朗日插值多项式导出的求积公式 将积分区间 a b 划分为n等分 记h b a n 取x0 a xk a kh k 0 1 n我们可以得到 3 牛顿 柯特斯求积公式 注释 在牛顿 柯特斯公式中 我们称c n k 为牛顿 柯特斯系数 一般可通过查表得到 崔国华教材p60列出了直到n 8的所有牛顿 柯特斯系数 应该说 实用意义不大 当n 8时牛顿 柯特斯公式并没有实际意义 实际上 我们通常只用到n 1 2 4的情形 相应的公式分别称为梯形公式 辛卜生公式 和柯特斯公式 对于现代的计算工具来说 有梯形公式和辛卜生公式也就够用了 利用牛顿 柯特斯系数 我们可以方便地写出牛顿 柯特斯求积公式 4 梯形公式 在牛顿 柯特斯求积公式中 如果我们取n 1 那么k可以取0和1 由此所形成的求积公式就是梯形公式 由于得到的结果是梯形面积公式 所以称梯形公式 5 辛卜生公式 在牛顿 柯特斯求积公式中 如果我们取n 2 那么k可以取0 1 2 由此所形成的求积公式就是辛卜生公式 6 柯特斯公式 作为课外作业 大家可以取n 4 相应地k可以取0 1 2 3和4 仿照上面的方式 可以得到 从而可进一步写出相应的求积公式 这就是柯特斯公式 在后面将要介绍的龙贝格求积算法中 我们将产生梯形序列 辛卜生序列 柯特斯序列和龙贝格序列 前三个序列都是基于牛顿 柯特斯公式产生的序列 而龙贝格序列则不是 10 3变步长复化梯形公式法 在上一节中我们介绍了牛顿 柯特斯公式以及它的特款 并且得到了牛顿 柯特斯公式的代数精度为n 接下来就是如何利用这个公式 一般说来 在一个较大的积分区间上利用较高阶的牛顿 柯特斯公式虽然可以得到较高的代数精度 但实际效果并不好 道理也不难理解 我们可以把积分区间划分为若干等分 在每个子区间上利用较低阶的求积公式 当然元可以是牛顿 柯特斯公式 也可以不是 并把这些积分值相加 按照这样的思路所得到的求积公式统称为复化型求积公式 1 复化中点公式 复化中点公式也许最不为人们所注意 以至在一般的教科书中还没有这个名称 我们在后面将会看到 对于求数值积分来说 它实际上是最有用的公式 把积分区间 a b 划分为n等分 记x0 x1 xn为等分点 记 xj 1 xj 为第j个子区间 zj为区间的中点 j 1 2 n 记h b a n 记Mn为所有子区间上利用中点公式所求得的积分值的和 那么我们有 1 1 复化中点公式 等价形式 现在我们把积分区间 a b 划分为2n等分 记y0 y1 y2n 1 y2n为等分点 此时我们可以把下标为奇数的点y2k 1 看成是区间 y2 k 1 y2k 的中点 k 1 2 n 如图 y0y1y2y3y4y2 n 1 y2n 1y2nx0 x1x2xn 1xn这样 我们也可以把复化中点公式写成 2 提示 这两种表示方法我们在后面都会用到 所以一定要记住 不管公式的形式如何 编写求Mn的程序还是相同的 1 复化中点公式 c语言代码 在实际应用中 n通常取2的某个整数幂2k 我们可以把计算结果作为数组M的第k个元素 程序段命名为 C语言源程序文件可以按如下的方式编写 include defineF X 4 0 1 0 X X staticdoublea 0 0 b 1 0 doubleGetM intk doublex y step intn j n 1 for j 0 j k j n 2 step b a n x a step 2 for j 0 j n j y F x x step return step y 2 复化梯形公式 把积分区间 a b 划分为n等分 记x0 x1 xn为等分点 记h b a n 则有记Ij j 1 2 n为f x 在第j个子区间 xj 1 xj 应用梯形公式所求得的积分值 则有 2 复化梯形公式 续 这就是我们的复化梯形公式 提示 虽然我们也可以直接编写求Tn的计算机程序 但是没有这个必要 3 变步长复化梯形公式 假设对某个n 我们利用复化梯形公式 也就是上面的 3 式 得到了Tn 如果它不满足我们的精度要求 那么我们可以把每个子区间再对分一次 这相当于把积分区间划分为2n等分 记y0 y1 y2 y2 n 1 y2n 1 y2n为等分点 记t b a 2n 则有再利用 3 式即得 3 变步长复化梯形公式 续 再利用上面的 2 式和 3 式 并注意到 3 式中的诸xk即为这里的y2k k 1 2 n 1 所以我们有 4 这就是变步长复化梯形公式 我们可以利用它形成一个自动调整精度的算法 4 变步长复化梯形公式方法 为了方便地利用 4 式形成一个算法 我们总是取n等于2的某个指数幂2k 并把Tn记为T k 由n 2k可得2n 2k 1 从而有T k 1 T K M K 2算法说明如下 T 0 b a f a f b 2 0 For K 0 k 20 k M k GetM k T k 1 M k T k 2 0 if fabs M k 1 M k EPS break 提示 如果取k 20 那么GetM k 中计算函数值的次数超过1000000次 所以我们一般不让k值超过20 12 4变步长复化辛卜生公式法 变步长复化梯形公式法的核心思想是 假设对于n 2K我们计算出来了T K 如果T k 不满足精度要求 则进一步把区间划分为2K 1等份 也就是在原来的每个子区间补上中点的函数值 并利用T k 1 T k M k 2来计算T k 1 现在我们也可以这样考虑 利用原来每个子区间的端点函数值f y2 k 1 和f x2k 以及后来补算的中点的函数值f y2k 1 采用辛卜生公式计算第k个子区间的积分值 再把他们相加结果计为S K 1 无论是理论分析还是直观考虑 这样做的效果应该更好一些 问题是能否也利用MK TK来简化我们的计算 1 复化辛卜生公式 假设x0 x1 x2 x2 k 1 x2k 1 x2k x2 n 1 x2n 1 x2n把积分区间 a b 划分为2n等分 我们可以认为其中的x0 x2 x2k x2n把区间划 a b 划分为n等份 并且x2k 1就是第k个子区间 x2 k 1 x2k 的中点 记Ik k 1 2 n为f x 在第k个子区间 x2 k 1 x2k 应用辛卜生公式所求得的积分值 则有这就是我们的复化辛卜生公式 虽然我们也可以直接编写求S2n的计算机程序 但是没有这个必要 2 关于复化辛卜生公式的结论 假设x0 x1 x2 x2 k 1 x2k 1 x2k x2 n 10 x2n 1 x2n把积分区间 a b 划分为2n等分 那么我们可以得到下面3个不同的积分值 利用x2k k 0 1 n这n 1个点处的函数值和复化梯形公式计算出Tn 利用xj j 0 1 2n这2n 1个点处的函数值和复化梯形公式计算出T2n 利用复化辛卜生公式计算出S2n 重要结论 对于上面的约定 我们有S2n 4T2n Tn 3 3 详细推导 4 变步长复化辛卜生公式加速算法 我们仍然取n为2的整数幂的形式 即n 2k 并记Sn为S k 那么我们有T k 1 M k T k 2S k 1 4T K 1 T K 3 0我们可以在变步长复化梯形公式法基础上稍加改进 即可得到一个相应的加速算法也就是简单地利用T K 1 和T K 的线性组合直接得到S k 1 额外的代码和计算量可忽略不计 5 算法说明 变步长复化辛卜生公式加速算法是一个实用的方法 而不是像前面介绍的方法只是过渡性的方法 所以要求大家都会编程 12 5龙贝格求积算法 前面我们介绍过变步长复化梯形公式法 利用它可以形成一个自动调整精度的算法 接下来又介绍了如何利用复化辛卜生公式实现变步长复化梯形公式法的加速 基本思想是利用T k 1 和T k 的线性组合来产生S k 1 我们称为进行一次外推 显然 它比起直接计算S k 1 优越的多 类似地 我们还可以利用S k 1 和S k 的线性组合来产生C k 1 称为柯特斯序列 利用C k 1 和C k 的线性组合来产生R k 1 称为龙贝格序列 显然这个过程还可以一直继续进行下去 这就是龙贝格算法的基本思想 1 统一的记号与公式 对于求f x 在区间 a b 上的定积分值 在龙贝格求积算法中 我们采用统一的记号它们之间的递推关系为 2 计算顺序 对于上面给出的外推序列 我们可按下图所示的顺序进行计算 如果有某个 即可停止进一步的计算 提示 为了便于编程 我们把每行的k值处理得相同 所以约定k m m 1 3 龙贝格算法的实现方法 在龙贝格算法中 我们实际的做法是把梯形序列外推3次 即m值取0 1 2 3 并把相应的序列称为梯形序列 辛卜生序列 柯特斯序列和龙贝格序列 并相应地记为T k S k C k R k 为了使我们的算法更为简明 我们再增加一个中点序列 记为M k 并且用一个专门的函数GetM intk 来计算 返回 M k 事实上 我们前面的变步长复化型辛卜生公式加速算法已经得到了中点序列 梯形序列和辛卜生序列 所以稍加改动即可得到龙贝格算法 4 实现龙贝格算法的代码 具体计算方法是 首先按公式组 1 设置M 0 T 0 S 0 C 0 R 0 接下来再对k 0 1 按公式组 2 计算M k T k 1 S k 1 C k 1 R k 1