计算方法引论-第五章

┆┆
┆
T2(0) T2(1) ┆ ┆ ┆
T3(0) ┆ ┆ ┆ ┆
误差 O(h2) O(h4) O(h6) O(h8) 系数 4,-1 42,-1 43,-1 44,-1
T (k) m
4 T m (k 1) m1
T (k) m1
4m 1
计算方法引论( 第三版)
5.18
徐萃薇、孙绳武 高教2007
Romberg求积公式(续)
5.14
徐萃薇、孙绳武 高教2007
事后误差估计 (续)
• 进一步的结果
– 逐次分半梯形法序列改进所得是逐次分 半抛物线法序列
– 逐次分半抛物线法序列又可改进得另一 个序列
– 这一过程可继续下去 – 改进也称外推,其依据是Euler-Maclaurin
公式
计算方法引论( 第三版)
5.15
徐萃薇、孙绳武 高教2007
计算方法引论( 第三版)
5.1
徐萃薇、孙绳武 高教2007
数值积分公式
• 数值积分公式
• 需求
b
n
f (x)dx
a
Ak f (xk )
k 0
– 原函数不一定能用初等函数表示 – 原函数虽能用初等函数表示,但表达式太复杂
不便计算
– 函数由表格给出
• 构造
– 可利用近似函数的积分作积分的近似，例如 由Lagrange插值多项式、Newton插值多项式 代入.
这里xk,Ak待定.乃有：
A1 A2 An 0
A1 x1 A2 x2 An xn 1
A1 x12
A2 x22
An
x
2 n
2,
A1 x1m A2 x2m An xnm m
k
b x k dx, k 0,1,m
a
计算方法引论( 第三版)
5.22
徐萃薇、孙绳武 高教2007
•例
– 梯形公式具有1次代数精确度
– 抛物线公式具有3次代数精确度
– 一般的Newton-Cotes公式的代数精确度
m=n,当n=2k+1 m=n+1,当n=2k
计算方法引论( 第三版)
5.7
徐萃薇、孙绳武 高教2007
代数精确度的确定
• 一个求积公式的代数精确度m可以从它们的截断
误差推出 ,也可按定义依次将1,x,x2,x3,…代入求积 公式从而检查出准确成立的最高次数.
n=n+n,s=0；
for i=1:2:n s=s+f(a+i*h)；
end
T=T0/2+h*s if | T-T0|<ε,break,end h=h/2；T0=T； end
计算方法引论( 第三版)
5.12
徐萃薇、孙绳武 高教2007
逐次分半梯形法(续)
•例
– 用逐次分半梯形法计算
–
I
结果
1
0
sin x
– 利用正交多项式可方便地给出结果
计算方法引论( 第三版)
5.23
徐萃薇、孙绳武 高教2007
• 说明
– 每列是一个逐次分半法序列,m=0是逐次分半梯形 公式的值以后是外推(改进)一次、外推
(改进)二次…的序列.
– 外推用公式,系数列于前列底.
– 计算可在算出梯形公式的新值后就外推至最右,即 按m+k=0,1,…自左而右计算.表中每行每列均收敛. 但若Euler-Maclaurin公式中某项为零或因光滑性不 够而中断时将不再有上述效果.
事后误差估计
• 判据 | T-T0|<ε
– 由I-Th≈Ch2及I-Th/2≈C(h/2)2相减得
Th/2-Th≈(4-1)C(h/2)2
– 估计误差 – 改进结果
I
Th / 2
Th / 2 Th 4 1
I 4Th / 2 Th 4 1
– 故ε应为预定误差容许值的三倍
计算方法引论( 第三版)
1 0.9397933 0.9460869 0.9460831
2 0.9445135 0.9460833
3 0.9456909 误差 O(h2) 系数 4,-1
O(h4) 42,-1
O(h6) 43,-1
O(h8) 44,-1
计算方法引论( 第三版)
5.20
徐萃薇、孙绳武 高教2007
算法
• 可在逐次分半梯形公式算法上增加外推並
– 梯形公式具有1次代数精确度(二法皆可) – 抛物线公式具有3次代数精确度(依定义) – 一般的Newton-Cotes公式的代数精确度(依定
义) m=n,当n=2k+1 m=n+1,当n=2k
计算方法引论( 第三版)
5.8
徐萃薇、孙绳武 高教2007
复化梯形公式
• 复化梯形公式
– 积分区间分成若干小区间在每个小区间 上用梯形公式即得复化梯形公式.
k 0
– 系数Ak有表可查. 从误差考虑应当采用系 数皆正的那些公式
计算方法引论( 第三版)
5.6
徐萃薇、孙绳武 高教2007
求积公式代数精确度
• 定义
– 一个求积公式对一切次数不超过m的多项式是 准确的,而有大于m的多项式不准确,则称该求积 公式具有m次代数精确度,或该求积公式的代数 精确度是m.
xdx
0.946083070
nT
n
T
1 0.92073549240395 2 0.93979328480618 4 0.94451352166539 8 0.94569086358270 16 0.94598502993439 32 0.94605856096277 64 0.94607694306006
Gauss型求积公式(续)
• 最高代数精确度的插值求积公式
– 当m=2n-1时方程数等于未知数个数.可望有 解(此前节点指定,系数待定时,m=n-1有解).求 得n个节点2n-1次代数精确度的插值求积公 式.
– 一般地
b
n
w(x) f (x)dx
a
Ak f (xk )
k 1
(其中权函数w(x)≥0,与各次多项式乘积的积 分存在,a,b皆可取∞)也能得到n个节点2n-1次 代数精确度的插值求积公式 .
– 误差
R(
f
)
h5 90
(
f
(4) (1 )
f
(4) (2 )
f
(4) (n
))
ba 180
h4
f
(4) ()
计算方法引论( 第三版)
5.10
徐萃薇、孙绳武 高教2007
逐次分半梯形法
• 复化梯形公式Tn与T2n的关系
令xk=a+kh,k=0,1, …,2n,h=(b-a)/(2n)可得
T2n
1 2 Tn
存储
n=1,h=(b-a)/2；T(1)=h(f(a)+f(b)) for k=2:kmax
n=n+n,s=0； for i=1:2:n
s=s+f(a+i*h)； end T(k)=T(k-1)/2+h*s; M=1; t=T(1); for j=k-1:-1:1
M=4*M;T(j)=(M*T(j+1)- T(j))/(M-1); end h=h/2；T if abs(T(1)-t)<ε,break,end end
计算 方法 引论： 数值 分析
误差 插值法与数值微分 数据拟合法 快速傅氏变换 数值积分
第五章 数值积分
• 插值求积公式 • Newton-Cotes公式 • 求积公式代数精确度 • 复合公式 • 逐次分半梯形法 • Richardson外推与 Romberg求积公式 • Gauss型求积公式
都能提高h2 这一过程可用Romberg求积公式实现
计算方法引论( 第三版)
5.17
徐萃薇、孙绳武 高教2007
Romberg求积公式
• Romberg求积公式(逐次分半外推法)
k m=0 m=1 m=2 m=3 …
0
T0(0)
T1(0)
1
T0(1)
T1(1)
2
T0(2)
T1(2)
3
T0(3) ┆
计算方法引论( 第三版)
5.21
徐萃薇、孙绳武 高教2007
Gauss型求积公式
• 问题
– n次插值构造的插值求 积公式至少有n次代数 精确度.那么，最高能 达多少次，又如何达 到.
– 考虑
b
n
f (x)dx
a
Ak f (xk )
– 代入1,x,x2,…,kxm1 看左右
相等最高的m是多少.
– 公式:令xk=a+kh,k=0,1,2, …,n,h=(b-a)/n
b a
f
(x)dx
h 2
(
f
(x0 )
2
f
( x1 )
2
f
(xn1 )
f
(xn ))
– 误差
R(
f
)
h3 12
(
f
(1 )
f
(2 )
f
(n ))
b a h2 f ()
12
பைடு நூலகம்
计算方法引论( 第三版)
5.9
徐萃薇、孙绳武 高教2007
• 称Ak为求积公式系数,R(f)为其截断断误差 • 易见对次数不超过n的多项式R(f)=0
计算方法引论( 第三版)
5.3
徐萃薇、孙绳武 高教2007
Newton-Cotes公式
• 梯形公式
– 用一次插值构造的插值求积公式称梯形 公式.几何上就是用梯形面积逼近曲边 梯形的面积
– 公式: 令 h=b-a
b
h
a f (x)dx 2 ( f (x0 ) f (x1))
– 误差
R( f ) b f (x )(x a)(x b)dx h3 f (), a b
a2
12
– 对次数不超过一的多项式梯形公式是准 确的
计算方法引论( 第三版)
5.4
徐萃薇、孙绳武 高教2007
Newton-Cotes公式(续)
计算方法引论( 第三版)
5.2
徐萃薇、孙绳武 高教2007
插值求积公式
• 将Lagrange插值公式代入积分:
n
f (x) f (xk )lk (x) Rn (x) k 0
b
n
f (x)dx
a
Ak f (xk ) R( f )
k 0
其中Ak,R(f)分别是lk(x),Rn(f)的积分.这样得 到的数值积分公式称插值求积公式.
– 停止判据:|Tm(0)-Tm-1(0)|<ε(容许误差)
计算方法引论( 第三版)
5.19
徐萃薇、孙绳武 高教2007
Romberg求积公式例
• Romberg求积公式例(I同逐次分半梯形法 例)
k
m=0
m=1
m=2
m=3
0 0.9207355 0.9461459 0.9460830 0.9460831
Euler-Maclaurin公式
• Euler-Maclaurin公式:
I-Th=a1h2+ a2h4 + a3h6 +… a1,a2,a3…与h无关,与f有关
• 用于事后误差估计及改进结果推导中
– 逐次分半梯形法序列改进一次误差是 O(h4)
– 外推(改进)两次所得序列的误差是O(h6) – 继续外推(改进)…每次外推序列的误差
h(
f
(x1)
f
(x3)
f
( x2 n1 ))
• 据此依次计算T1,T2, T4…
– T1,梯形公式 – 重复进行:区间分半,积分值之半加上新点上函
数值之和与h之积.
• 计算可用误差控制並限定分半次数
计算方法引论( 第三版)
5.11
徐萃薇、孙绳武 高教2007
逐次分半梯形法(续)
•算法
n=1,h=(b-a)/2；T0=h(f(a)+f(b)) for k=1:kmax
复化抛物线公式
• 复化抛物线公式
– 积分区间分成2n个小区间在每两个小区间 上用抛物线公式即得复化抛物线公式.
– 公式:令xk=a+kh,k=0,1, …,2n,h=(b-a)/(2n)
b a
f
(x)dx
h 3
(
f
(x0 )
4
f
(x1)
2
f
(x2 )
2 f (x2n2 ) 4 f (x2n1) f (x2n ))
5.5
徐萃薇、孙绳武 高教2007
Newton-Cotes公式(续)
• 一般的Newton-Cotes公式
– 取xk=a+kh,k=0,1,2, …,n,h=(b-a)/n,作n次 插值,构造的插值求积公式通称Newton-
Cotes公式亦称等距基点插值求积公式
b
n
f (x)dx
a
Ak f (xk )
都能提高h2
计算方法引论( 第三版)
5.16
徐萃薇、孙绳武 高教2007
Richardson外推
• 根据Euler-Maclaurin公式,可作Th和Th/2的线
性组合得到误差O(h4)的值称Richardson外 推(前称改进结果)于是
– 逐次分半梯形法序列外推(改进)一次误 差是O(h4)
– 外推(改进)两次误差是O(h6) – 继续外推(改进)…每次外推序列的误差
•抛物线公式
– xk=a+kh,k=0,1,2,h=(b-a)/2,作二次插值构 造的插值求积公式称抛物线公式.
b
a
f
(x)dx
h 3
(
f
(x0 )
4
f
(x1)
f
( x2 ))
– 误差 R( f ) h5 f (4) (), a b
90
– 对次数不超过三的多项式抛物线公式是 准确的
计算方法引论( 第三版)
128 256 512 1024 2048 4096
0.94608153854315 0.94608268741135 0.94608297462823 0.94608304643245 0.94608306438350 0.94608306887126
计算方法引论( 第三版)
5.13
徐萃薇、孙绳武 高教2007