数值积分及matlab实现PPT课件.ppt

数值积分与微分 2009 4 22 1 数值积分和数值微分 1引言我们知道 若函数f x 在区间 a b 上连续且其原函数为F x 则可用Newton Leibnitz公式 求得定积分 求定积分的值 Newton Leibnitz公式无论在理论上还是在解决实际问题上都起了很大作用 但它并不能完全解决定积分的计算问题 因为积分学涉及的实际问题极为广泛 而且极其复杂 在实际计算中经常遇到以下三种情况 2 1 被积函数f x 并不一定能够找到用初等函数的有限形式表示的原函数F x 例如 Newton Leibnitz公式就无能为力了 2 还有被积函数f x 的原函数能用初等函数表示 但表达式太复杂 例如函数 并不复杂 但积分后其表达式却很复杂 积分后其原函数F x 为 3 3 被积函数f x 没有具体的解析表达式 其函数关系由表格或图形表示 对于这些情况 要计算积分的准确值都是十分困难的 由此可见 通过原函数来计算积分有它的局限性 因而研究一种新的积分方法来解决Newton Leibniz公式所不能或很难解决的积分问题 这时需要用数值解法来建立积分的近似计算方法 将积分区间细分 在每一个小区间内用简单函数代替复杂函数进行积分 这就是数值积分的思想 用代数插值多项式去代替被积函数发f x 进行积分是本章讨论数值积分的主要内容 4 建立数值积分公式的途径比较多 其中最常用的有两种 1 由积分中值定理可知 对于连续函数f x 在积分区间 a b 内存在一点 使得即所求的曲边梯形的面积恰好等于底为 b a 高为的矩形面积 但是点 的具体位置一般是未知的 因而的值也是未知的 称为f x 在区间 a b 上的平均高度 那么只要对平均高度提供一种算法 相应地就获得一种数值求积方法 5 三个求积分公式 梯形公式 y f x y x a b y f x a b y x a b 2 中矩形公式 按照这种思想 可构造出一些求积分值的近似公式 例如分别取和 则分别得到中矩形公式和梯形公式 y f x a b a b 6 y f x y a b Simpson公式 a b 2 f 的近似值而获得的一种数值积分方法 中矩形公式把 a b 的中点处函数值作为平均高度f 的近似值而获得的一种数值积分方法 a b a b 2 在这三个公式中 梯形公式把f a f b 的加权平均值 作为平均高度 7 Simpson公式是以函数f x 在a b a b 2这三点的函数值f a f b 的加权平均值似值而获得的一种数值积分方法 作为平均高度f 的近 2 先用某个简单函数近似逼近f x 用代替原被积函数f x 即 以此构造数值算法 从数值计算的角度考虑 函数应对f x 有充分的逼近程度 并且容易计算其积分 由于多项式能很好地逼近连续函数 且又容易计算积分 因此将选取为插值多项式 这样f x 的积分就可以用其插值多项式的积分来近似代替 8 2 2插值求积公式设已知f x 在节点有函数值 作n次拉格朗日插值多项式 式中 这里 多项式P x 易于求积 所以可取作为的近似值 即 9 其中 称为求积系数 给出如下定义 定义1求积公式 其系数时 则称求积公式为插值求积公式 4 10 设插值求积公式的余项为 由插值余项定理得 其中 当f x 是次数不高于n的多项式时 有 0 求积公式 4 能成为准确的等式 由于闭区间 a b 上的连续函数可用多项式逼近 所以一个求积公式能对多大次数的多项式f x 成为准确等式 是衡量该公式的精确程度的重要指标 为此给出以下定义 11 定义2 代数精度 设求积公式 4 对于一切次数小于等于m的多项式 是准确的 而对于次数为m 1的多项式是不准确的 则称该求积公式具有m次代数精度 简称代数精度 或 12 定理1n 1个节点的求积公式为插值型求积公式的充要条件是公式至少具有n次代数精度 13 例1设积分区间 a b 为 0 2 取时时 分别用梯形和辛卜生公式 计算其积分结果并与准确值进行比较解 梯形公式和辛卜生的计算结果与准确值比较如下表所示 14 f x 1xx2x3x4ex准确值222 6746 406 389梯形公式计算值2248168 389辛卜生公式计算值222 6746 676 421 从表中可以看出 当f x 是时 辛卜生公式比梯形公式更精确 一般说来 代数精度越高 求积公式越精确 梯形公式和中矩形公式具有1次代数精度 辛卜生公式有3次代数精度 下面以梯形公式为例进行验证 15 取f x 1时 两端相等 取f x x时 取f x x2时 两端不相等 所以梯形公式只有1次代数精度 两端相等 16 构造插值求积公式有如下特点 复杂函数f x 的积分转化为计算多项式的积分求积系数Ak只与积分区间及节点xk有关 而与被积函数f x 无关 可以不管f x 如何 预先算出Ak的值n 1个节点的插值求积公式至少具有n次代数精度求积系数之和可用此检验计算求积系数的正确性 17 3牛顿 柯特斯 Newton Cotes 求积公式在插值求积公式 中 当所取节点是等距时称为牛顿 柯特斯公式其中插值多项式求积系数 这里是插值基函数 即有 18 将积分区间 a b 划分为n等分 步长求积节点为为了计算系数Ak 由于 所以 作变量代换当时 有 于是可得 19 k 0 1 n 代入插值求积公式 4 有 称为牛顿 柯特斯求积公式 Ck称为柯特斯系数 引进记号 k 0 1 n 则 20 容易验证 显然 Ck是不依赖于积分区间 a b 以及被积函数f x 的常数 只要给出n 就可以算出柯特斯系数 譬如当n 1时 21 当n 2时 22 2020 1 8 23 4几个低阶求积公式在牛顿 柯特斯求积公式中n 1 2 4时 就分别得到下面的梯形公式 辛卜生公式和柯特斯公式 1 梯形公式当n 1时 牛顿 柯特斯公式就是梯形公式 定理2 梯形公式的误差 设f x 在 a b 上具有连续的二阶导数 则梯形公式的误差 余项 为 24 2 辛卜生公式当n 2时 牛顿 柯特斯公式就是辛卜生公式 或称抛物线公式 定理3 辛卜生公式的误差 设在 a b 上具有连续的四阶导数 则辛卜生求积公式的误差为 定理证明从略 25 3 柯特斯公式 当n 4时 牛顿 柯特斯公式为 定理4 柯特斯公式的误差 设在 a b 上具有连续的6阶导数 则柯特斯求积公式的误差为 定理的证明从略 26 例11分别用梯形公式 辛卜生公式和柯特斯公式计算定积分的近似值 计算结果取5位有效数字 1 用梯形公式计算 2 用辛卜生公式 27 3 用柯特斯公式计算 系数为 积分的准确值为 可见 三个求积公式的精度逐渐提高 28 例12用辛卜生公式和柯特斯公式计算定积分 的近似值 并估计其误差 计算结果取5位小数 解 辛卜生公式 由于由辛卜生公式余项 知其误差为 29 例12用辛卜生公式和柯特斯公式计算定积分 的近似值 并估计其误差 计算结果取5位小数 解 柯特斯公式 知其误差为 30 例12用辛卜生公式和柯特斯公式计算定积分 的近似值 并估计其误差 计算结果取5位小数 该定积分的准确值 这个例子告诉我们 对于同一个积分 当n 2时 公式却是精确的 这是由于辛卜生公式具有三次代数精度 柯特斯公式具有五次代数精度 它们对被积函数为三次多项式当然是精确成立的 31 数值积分基本原理求解定积分的数值方法多种多样 如简单的梯形法 辛普生 Simpson 法 牛顿 柯特斯 Newton Cotes 法等都是经常采用的方法 它们的基本思想都是将整个积分区间 a b 分成n个子区间 xi xi 1 i 1 2 n 其中x1 a xn 1 b 这样求定积分问题就分解为求和问题 32 数值积分的实现方法1 变步长辛普生法基于变步长辛普生法 MATLAB给出了quad函数来求定积分 该函数的调用格式为 I n quad fname a b tol trace 其中fname是被积函数名 a和b分别是定积分的下限和上限 tol用来控制积分精度 缺省时取tol 0 001 trace控制是否展现积分过程 若取非0则展现积分过程 取0则不展现 缺省时取trace 0 返回参数I即定积分值 n为被积函数的调用次数 33 例1求定积分 1 建立被积函数文件fesin m functionf fesin x f exp 0 5 x sin x pi 6 2 调用数值积分函数quad求定积分 S n quad fesin 0 3 pi S 0 9008n 77 34 2 牛顿 柯特斯法基于牛顿 柯特斯法 MATLAB给出了quad8函数来求定积分 该函数的调用格式为 I n quad8 fname a b tol trace 其中参数的含义和quad函数相似 只是tol的缺省值取10 6 该函数可以更精确地求出定积分的值 且一般情况下函数调用的步数明显小于quad函数 从而保证能以更高的效率求出所需的定积分值 35 例2求定积分 1 被积函数文件fx m functionf fx x f x sin x 1 cos x cos x 2 调用函数quad8求定积分 I quad8 fx 0 pi I 2 4674 36 例3分别用quad函数和quad8函数求定积分的近似值 并在相同的积分精度下 比较函数的调用次数 调用函数quad求定积分 formatlong fx inline exp x I n quad fx 1 2 5 1e 10 I 0 28579444254766n 65 37 调用函数quad8求定积分 formatlong fx inline exp x I n quad8 fx 1 2 5 1e 10 I 0 28579444254754n 33 38 3 被积函数由一个表格定义在MATLAB中 对由表格形式定义的函数关系的求定积分问题用trapz X Y 函数 其中向量X Y定义函数关系Y f X 例4用trapz函数计算定积分 命令如下 X 1 0 01 2 5 Y exp X 生成函数关系数据向量trapz X Y ans 0 28579682416393 39 1 3二重定积分的数值求解使用MATLAB提供的dblquad函数就可以直接求出上述二重定积分的数值解 该函数的调用格式为 I dblquad f a b c d tol trace 该函数求f x y 在 a b c d 区域上的二重定积分 参数tol trace的用法与函数quad完全相同 40 例5计算二重定积分 1 建立一个函数文件fxy m functionf fxy x y globalki ki ki 1 ki用于统计被积函数的调用次数f exp x 2 2 sin x 2 y 2 调用dblquad函数求解 globalki ki 0 I dblquad fxy 2 2 1 1 kiI 1 57449318974494ki 1038 41 2数值微分2 1数值差分与差商2 2数值微分的实现在MATLAB中 没有直接提供求数值导数的函数 只有计算向前差分的函数diff 其调用格式为 DX diff X 计算向量X的向前差分 DX i X i 1 X i i 1 2 n 1 DX diff X n 计算X的n阶向前差分 例如 diff X 2 diff diff X DX diff A n dim 计算矩阵A的n阶差分 dim 1时 缺省状态 按列计算差分 dim 2 按行计算差分 42 例6生成以向量V 1 2 3 4 5 6 为基础的范得蒙矩阵 按列进行差分运算 命令如下 V vander 1 6 DV diff V 计算V的一阶差分 43 例7用不同的方法求函数f x 的数值导数 并在同一个坐标系中做出f x 的图像 程序如下 f inline sqrt x 3 2 x 2 x 12 x 5