880 likes | 1.05k Views
Chapter 7 数值积分与数值微分. 内容提纲( Outline). 求积公式的代数精度 插值型求积公式 复化求积法. Why do we do numerical integral?. 为什么要数值积分? 在微积分里,按 Newton-Leibniz 公式求定积分 要求被积函数 f ( x ) ☞ 有解析表达式; ☞ f ( x ) 的原函数 F ( x ) 为初等函数.. 问题 ☎ f ( x ) 没有解析表达式,只有数表形式 e.g. ☎ f ( x ) 有表达式,但原函数不是初等函数
E N D
内容提纲(Outline) • 求积公式的代数精度 • 插值型求积公式 • 复化求积法
Why do we do numerical integral? 为什么要数值积分? 在微积分里,按Newton-Leibniz公式求定积分 要求被积函数f(x) ☞ 有解析表达式; ☞f(x)的原函数F(x)为初等函数.
问题 ☎ f(x)没有解析表达式,只有数表形式 e.g. ☎ f(x)有表达式,但原函数不是初等函数 e.g. , 它们的原函数都不是初等函数.
求定积分就得通过近似计算-数值积分求得积分近似值求定积分就得通过近似计算-数值积分求得积分近似值 基本思想是对被积函数进行近似,给出数值积分,同时考虑近似精度。 下面首先给出代数精确度的概念
7.1 代数精确度 本章讨论的是形如 的定积分的数值计算,其中 为权函数, 要满足5.4节中所提的条件.
一般把积分区间n个点{xk}上的函数值f(xk)加权Ak的和 作为积分I(f)的近似, 即 或记 (2)
上式中xk,Ak分别称为求积节点、求积系数.求积系数与被积函数f(x)无关,而与求积节点、求积区间、权函数有关.称公式(2)为n点求积公式,有时也称上式中xk,Ak分别称为求积节点、求积系数.求积系数与被积函数f(x)无关,而与求积节点、求积区间、权函数有关.称公式(2)为n点求积公式,有时也称 为一个n点求积公式, 为求积公式的误差.用此公式)求积分近似值的计算称为数值积分或数值微分.
构造或确定一个求积公式,要讨论解决的问题有构造或确定一个求积公式,要讨论解决的问题有 (i) 确定求积系数Ak和求积节点n; (ii) 求积公式的误差估计和收敛性. 用什么标准来判定两个节点数相同的求积公式的“好”与“差”呢?通常用“代数精确度”的高低作为求积公式“好”与“差”的一个标准.在后面的讨论中我们将看到,节点相同的求积公式,代数精确度越高,求出的积分近似值精确度一般越好.下面给出代数精确度的定义.
定义1 若对任意的 ,求积公式(2)的误差都满足 ,则称该求积公式具有n次代数精确度. 定义1 若对任意的 ,求积公式(2)的误差都满足 ,则称该求积公式具有n次代数精确度. 验证一个求积公式所具有的代数精确度用定义1是极不方便的,为此给出另一个定义.
定义2 若对函数 , 求积公式(2)精确成立,即 而 , 则称其具有n次代数精确度. 因为函数组 是 的一组基函数,所以两个定义是等价的,但在具体应用时,定义2比定义1要方便的多.
例1 验证求积公式 例1 验证求积公式 具有3次代数精确度. 解: 当 而 有
(1)当 (2)当 (3)当
(1)当 故求积公式具有三次代数精确度.
7.2 插值型求积公式 这一节所讨论的求积公式,都是用在区间[a, b]上对被积函数f(x)作插值所得插值多项式Pn(x)代替被积函数f(x)导出的公式.这一类求积公式的求积节点 xk,就是对f(x)作插值时的插值节点,所以这类求积公式称为插值型求积公式. 为简便起见,这节讨论节点分布为等距并且权函数 时的插值型求积公式的构造等问题.
7.2.1 Newton-Cotes求积公式 一、公式的推导 设将积分区间[a,b]n等分,求积节点为 , 那么, 令x=a+th,则t=(x-a)/h,且由 可知 . 由Lagrange插值基函数有 而 ,所以
将n次Lagrange插值多项式Ln(x)代替被积函数f(x)得 记 称为Cotes求积系数.它与(3)式中的求积系数Ak相差一个常数b-a 即
把Ak代入到(3)式中,得到Newton-Cotes求积公式.例如把Ak代入到(3)式中,得到Newton-Cotes求积公式.例如 当n=4,5时,Newton-Cotes公式分别为 n=0,1,2三种情形,在讨论(3)式中的余项R(1, f)后再详细讨论.
二、误差估计 求积公式(3)计算出的积分I(f)的近似值In+1(f)的误差多大? 若被积函数 ,记 , 对n次Lagrange插值余项求积,可得n+1个节点的Newton-Cotes求积公式的误差估计式为 (5)
验证求积公式(3)的代数精确度,不用误差估计的(4)式,而用直接对插值余项求积的形式,即 (5) 由(5)式,显而易见,当 时,因 可知,R(1, f)=0,所以我们所n+1点的求积公式(3)至少具有n次的代数精确度.进一步可以证明,当n为偶数时,求积公式(3)的代数精确度可以达到n+1次.
三、几种常见的Newton-Cotes求积公式 对 n=0, 1, 2, 按公式(3)可以得出下面三种常见的Newton-Cotes求积公式. 1. n=0时的矩形求积公式 分别以积分区间[a, b]的左、右端点和区间中点,即x=a,b,(a+b)/2为求积节点得到: 左矩形求积公式: 右矩形求积公式: 中矩形求积公式: 三个求积公式的误差估计,可将函数f(x) 分别在 处展开到含f(x)的一阶导数的Taylor公式在区间[a, b]上积分 推得.
2. n=1时的梯形求积公式 按Cotes系数公式计算得 故求积系数A0, A1为 , 梯形求积公式为 记 (6)式的几何意义如图7-2所示(见p327) 容易验证公式(6)的代数精确度的次数为1. 考虑梯形求积公式(6)的误差估计R(1, f).假定 时,用推广的积分中植定理,将过(a, f(a)), (b, f(b))点的线性 插值的余项 在[a, b]上积分,可得 其中 也称为梯形求公式
3. n=2时的Simpson求积公式 按Cotes系数公式可以计算出 为此 , , 所以 (8) 公式(8) 称为Simpson求积公式.由7.1节例1可知 Simpson求积公式(8)具有3次的代数精确度. Simpson求积公式(8)的误差估计R(1, f)不能直接有插值 余项 利用推广的积分中值定理在 [a, b]上积分推出.原因是 在[a, b]上要变号.
7.3 Romberg积分 7.3.1 Richardson外推法 外推法是用精确度较低的近似公式组合成精确度较高的近似公式的一种方法. 设 是任意数, 是关于步长h逼近 的近似公式,它们的误差估计式为 (1) 这里, k1,k2,k3,… 是一组常数.
按(1)式,称 逼近 的误差为 .把h 的幂次称为误差的阶,例如, 称为二阶误差 我们希望找到一种简便的方法,用近似公式F(h)的组合,得到误差阶较高的近似公式 ,使 (2) 此时, 逼近 F*的误差为O(h2) 类似地,用 组合产生逼近F*的误差 为 O(h3) 的近似公式等.下面我们给出一种具体的组合方法.
把(1)式改写为 (3) 用h/2代替(3)式中的h,得 (4) 用2乘(4)式再减去(3)式,消去含h的项,得 (5) 令 ,且记
那么(5)式可写为 (6) 这里, 逼近 的误差为 再用 h/2 代替 h , 使(6)式变为 (7) 用4乘(7)式减去(6)式,消去含 的项,得 (8) 同样记
(8)式可以写为 (9) 这里 逼近 的误差为 还是用h/2代替h代入(9)式后,类似上述过 程,可以得到误差为 的 一般地,对 ,有逼近 的误差为 的递推公式 (10) 也称为关于步长h的外推公式. 表7-1列出了 时,按(10)式产生 的计算次序,表中各列左边黑体数字表示序号.
表7-1 例1 设 带余项的差分公式为
(11) 导出具有误差为 的外推公式. 解 令 用 h/2代替h,得 (12) 为消去含 的项,用4乘(12)式减去 (11)式,得
从而有 (13) 这里 这时, 逼近 的误差为 . 重复用h/2代替h并消去含 的项 ,得到逼近 的误差为 的 外推公式为
注意(14)式中第二项的分母为 而不 是(10)式中的 .这是由于(11)式中的余项 为关于 的幂次而不是关于h的幂次. 7.3.2 Romberg求积方法 Romberg求积方法是以复化梯形公式为基 础,应用Richardson外推法导出的数值求积方 法. 回忆7.2.1节的复化梯形公式,分别把积分区
间[a,b]分为1,2,4等分的结果列入表7-2. 表7-2 k 1 1 2 2 3 4
我们还可以进一步推导出它们的递推关系.由 可以化为 类似地,有 一般地,把区间[a,b]逐次分半k-1次 , 区间长度(步长)为 ,其中 .为
叙述方便起见,记 ,那么, (15) 或 (16) 从而有 (17) 其中 . 按外推法的思想,可以把(15)看成是关于
误差为 的一个近似公式.因此,复化梯形公 式的误差公式为 (18) 为消去 项,再取 代替(18)式中的 , 得 (19)
用4乘(19)式再减去(18)式,得 (20) 记 (21) 这是误差为 的外推公式. 重复上述过程,将区间逐次分半k-1次后,可
以得到误差为 的外推公式 (22) 当j=2时 (23) 当K=2时,有 这是n=2的复化Simpson公式的 .不难验证,对 一般的k, ,这里, 是 的复化
Simpson公式. 类似地,当j=3时, (24) 在实际计算中,经常直接应用(23)式和形式与(24) 式相类似的公式进行计算. 所谓Romberg求积方法,就是由上述两部分 组成.第一部分,对积分区间逐次分半k-1次,用复 化梯形求积公式(16)计算 ,第二部 分,用外推公式(22)计算
用Romberg求积方法计算 的计算值的 过程如下:首先,令k=1,区间长度 ,用梯形 求积公式计算 (表7-3中第一行);区间分半,令 K=2,区间长度 ,先按(16)式计算 ,再按 外推公式(22)式计算 (表7-3中第二行);再区 间分半,令k=3,区间长度 ,先按(16) 式计算 ,再按(22)式计算 (表7-3中第 三行)等等,逐次分半区间k次后的计算结果如表 7-3所示(见下页).
表7-3 : : : : ……. 表7-3中 的计算按行(k的序号)进行,每行第 1个元素 用复化梯形公式(16)计算,其他元 素 均按(22)式用 与 的组 合得到.在实际应用中,往往根据实际问题对计
算精确度的要求来确定区间逐次分半的次数. 常用不等式 (25) 作为达到精确度要求的判断准则,这里, 是给 定的一个小的正数. 例2用Romberg求积方法计算 (26) 的近似值,给定 解 首先令区间长度 ,用梯形求积公式计算
区间[0,1]分半,令区间长度 ,按16式 计算 再按(23)式计算 这时 未达到精确度要求. 为此,再将区间分半,令区间长度 按(16)式计算
按j=2和j=3的外推公式(23)和(24),分别用 和 的组合得到 以及用 和 的组合 得到 ,即 以及 这时, 已满足不等式(25)的要求. 作为积分(26)式 的近似,其误差为 .
下面给出用Romberg求积方法计算 近似值的计算步骤,用二维数组T的元素 存 放表7-3中的 . • 输入:积分区间端点 • 令 ,计算 • 令 • 令 ,计算 • for j=2,3,…,k 5.1 计算 end for (j)
if ,then goto 8 end if 7 令k=k+1, ,goto 4 8 输出 9 end
7.5 Gauss求积公式 7.5.1 引言 求积公式 (1) 当求积系数 、求积节点 都可以 自由选取时,其代数精确度最高可以达到多少 次? 下面的引理可以回答上述问题.
引理1当求积系数 和求积节点 都可以自由选取时,n点的求积公式(1)的代数精 确度最高可以达到2n-1次. 证 假设求积公式(1)具有m次代数精确度, 即对任意的m次代数多项式 求积公式(1)的精确成立.于是成立等式 即 若记 (2)
则(2)式成为 (3) 由于系数 的任意性,故使(3)式 成为恒等式的充要条件是 (4) (4)式的待定系数有2n个,所以确定待定系数的
独立条件至多给出2n个,从而可知m至多为 2n-1. 定义1 n点的求积公式(1)具有2n-1次代数 精确度(或称为具有最高的代数精确度)时,称 为Gauss型求积公式. Gauss型求积公式的求积节点 ,称为 Gauss点,它们可以通过求区间[a,b]上带权 的n次正交多项式 的n个根获得.所以先介 绍正交多项式及其性质.然后讨论Gauss型求 积公式的构造,等等.