4.4龙贝格求积公式
龙贝格求积公式
2
f( ) 2
T3 , 1 T3 , 2
T3 , 3
T4 , 1 T5 , 1
M
T4 , 2 T5 , 2
M
T4, 3 T5 , 3
M
T4 , 4 T5 , 4
M
用龙贝格方法计算积分的步骤为：
（1）：准备初值，先用梯形公式计算积分近
似 （值 2）：：T1按 变b 2步a[长f (a梯) 形f (公b)]式计算积分近似值：
4 0.918741799 0.916327874 0.916297224 0.916294351
5 0.916905342 0.916293190 0.916290077 0.916290776
T5 , 4 0.916290776
例3：取=0.00001,用龙贝格方法计算积分
1
I
4
dx
01 x2
解：由题意
f(x)=4/(1+x2) a=0 b=1 f(0)=4 f(1)=2 由梯形公式得
T1=1/2[f(0)+f(1)]=3 计算f(1/2)=16/5 用变步长梯形公式得
T2=1/2[T1+f(1/2)]=3.1 由加速公式得
S1=1/3(4T2-T1)=3.133333333
求出f(1/4) f(3/4) 进而求得 T4=1/2{T2+1/2[f(1/4)+f(3/4)]}
Simpson加速公式： Cn
42 S2n Sn 42 1
I
C2n
1 43 1(C2n
Cn )
43C2n Cn 43 1
Cotes加速公式：Rn
43C2n Cn 43 1
Romberg 值序列
4.4龙贝格求积公式
4 1 T1 ( k − 1) = T0 ( k ) − T0 ( k − 1) 3 3 16 1 T2 ( k − 1) = T1 ( k ) − T1 ( k − 1) 15 15 64 1 T3 ( k − 1) = T2 ( k ) − T2 ( k − 1) 63 63
k = 1 ,2 ,L
因此有如下递推公式 b−a [ f ( a ) + f (b )] T0 (0) = 2
1 T0 (k ) = T0 (k − 1) + hk 2
2 k −1 −1 j =0
∑ f (a + (2 j + 1)h )
k
k = 1, 2 ,L
上式称为递推的梯形公式
思考
递推梯形公式加上一个控制精度,即 可成为自动选取步长的复合梯形公式
带权）正交。 不大于n 不大于 的多项式 P(x) （带权）正交。
k =0
n
x 证明： 证明： “⇒”0 … xn 为 Gauss 点, 则公式 ∫a ρ( x) f ( x)dx ≈ ∑Ak f ( xk ) k=0 次代数精度。 至少有 2n+1 次代数精度。 对任意次数不大于 的多项式 不大于n 对任意次数不大于 ⇔ 求w(x) Pm(x)， Pm(x) w(x)的次数 ， 的次数 求 Gauss 点 不大于2n+1，则代入公式应精确成立： 精确成立： 不大于 ，则代入公式应精确成立 n 0 b ∫ ρ ( x ) Pm ( x ) w ( x )dx = ∑ Ak Pm ( x k ) w ( x k ) = 0
外推加速公式
由复合梯形公式的余项公式
I − T2 n 1 ≈ I − Tn 4 1 I − T2 n ≈ (T2 n − Tn ) 3 4 1 I ≈ T2 n − Tn 3 3 1 f ( x 1 )) − Tn j+ 3 2
龙贝格算法
龙贝格积分1. 算法原理采用复化求积公式计算时，为使截断误差不超过ε，需要估计被积函数高阶导数的最大值，从而确定把积分区间[]b a ,分成等长子区间的个数n 。
首先在整个区间[]b a ,上应用梯形公式，算出积分近似值T1；然后将[]b a ,分半，对 应用复化梯形公式算出T2；再将每个小区间分半，一般地，每次总是在前一次的基础上再将小区间分半，然后利用递推公式进行计算，直至相邻两个值之差小于允许误差为止。
实际计算中，常用ε≤-n n T T 2作为判别计算终止的条件。
若满足，则取n T f I 2][≈；否则将区间再分半进行计算，知道满足精度要求为止。
又经过推导可知，∑=-++=ni i i n n x x f h T T 112)2(221，在实际计算中，取kn 2=，则k a b h 2-=，112)1*2(2++--+=+k i i ab i a x x 。
所以，上式可以写为∑=++--+-+=+kk i k k ab i a f a b T T 211122)2)12((2211k开始计算时，取())()(21b f a f ab T +-=龙贝格算法是由递推算法得来的。
由梯形公式得出辛普森公式得出柯特斯公式最后得到龙贝格公式。
根据梯形法的误差公式，积分值n T 的截断误差大致与2h 成正比，因此步长减半后误差将减至四分之一，即有21114n n T T -≈-将上式移项整理，知2211()3n n n T T T -≈-由此可见，只要二分前后两个积分值n T 和2n T 相当接近，就可以保证计算保证结果计算结果2n T 的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计法。
按上式，积分值2n T 的误差大致等于21()3n n T T -，如果用这个误差值作为2n T 的一种补偿，可以期望，所得的()222141333n n n n n T T T T T T =+-=-应当是更好的结果。
龙贝格算法
2）理查森外推加速 从梯形公式出发， 将区间[a,b]逐次二分可以提高求积公式精度， 当[a,b] 分为 n 等份时，若记T������ = T ℎ ，当区间[a,b]分为 2n 等份时，则有 T2������ = T
ℎ 2
。再有泰勒公式展开为： T ℎ = ������ + ������1 ℎ2 + ������2 ℎ+ ⋯ +
建立一个命名为Romberg.m的function文件：
function[T]=Romberg(f,a,b,e) T=zeros(10,10); T(1,1)=(b-a)/2*(f(a)+f(b)); for k=2:10, sum=0; for i=1:2^(k-2), x=a+(2*i-1)*(b-a)/2^(k-1); sum=sum+f(x); end T(k,1)=T(k-1,1)/2+[(b-a)/2^(k-1)]*sum; for j=2:k, end %第一列用递推梯形公式 %定义龙贝格函数 %定义10阶的零元矩阵
1 3/2 ������ 0
������������。
算，并取������ ������, ������ ≈ ������ ；否则令 k=k+1，转（3）继续计算。 6）下图为我按照自己的算法所设计的示意表： 算法设计表：
k 1 2 3 4 …
h b-a (b-a)/2 (b-a)/4 (b-a)/8 …
������ ������, 1 T(1,1) T(2,1) T(3,1) T(4,1)
������ ������, 2
������ ������, 3
������ ������, 4
������ ������, 5
数值分析63 复化求积公式龙贝格求积公式讲解
只增加了一个分点
1 xk?1/ 2 ? 2 ( xk ? xk?1)
设hn=(b? a)/n, xk=a+kh n (k=0,1,? ,n),在[xk, xk+1] 上用梯形公式得
T1 ?
hn 2
?f
(
xk
)
?
f ? ( xk ? 1 )
复化求积的基本想法 :
将积分区间 [a, b]n等分, 步长
h?
b
? n
a
,
分点为
xk=a+kh (k=0,1,…,n) , 则由定积分性质知
? ? ? I ?
b
n?1
f ( x )dx ?
xk?1 f ( x )d x
a
k ? 0 xk
每个子区间 上的积分
?xk?1 f ( x )dx xk
用低阶求积公式 , 然后把所有区间的 计算结果求和 ,
注2: 同样也可用 | S4m-S2m |<ε 来控制计算的精度 . 这就是下面要介绍的 龙贝格求 积公式 .
6.4 龙贝格求积公式
6.4.1 梯形公式的递推化
复化求积方法可提高求积精度，实际计算时若
精度不够可将步长逐次分半 . 设将区间 [a, b]分为n等
分，共有 n+1个分点，如果将求积区间再分一次，则 分点增至 2n+1个，我们将二分 前后两个积分值 联系
果T8=0.9456909 只有2位有效数字，而应用复化辛普 森公式计算的结果 S4= 0.9460832 却有6位有效数字 .
注：为了利用余项公式估计误差，要求 f(x)=sin x/x 的高阶导数，由于
数值分析--第4章数值积分与数值微分[1]详解
![数值分析--第4章数值积分与数值微分[1]详解](https://img.taocdn.com/s3/m/9edd6ad82f60ddccdb38a082.png)
第4章 数值积分与数值微分1 数值积分的基本概念实际问题当中常常需要计算定积分。
在微积分中，我们熟知，牛顿-莱布尼兹公式是计算定积分的一种有效工具,在理论和实际计算上有很大作用。
对定积分()ba I f x dx =⎰，若()f x 在区间[,]ab 上连续，且()f x 的原函数为()F x ,则可计算定积分()()()baf x dx F b F a =-⎰似乎问题已经解决，其实不然。
如1）()f x 是由测量或数值计算给出数据表时,Newton-Leibnitz 公式无法应用。
2)许多形式上很简单的函数，例如222sin 1(),sin ,cos ,,ln x x f x x x e x x-= 等等，它们的原函数不能用初等函数的有限形式表示。
3）即使有些被积函数的原函数能通过初等函数的有限形式表示，但应用牛顿—莱布尼兹公式计算，仍涉及大量的数值计算，还不如应用数值积分的方法来得方便，既节省工作量，又满足精度的要求。
例如下列积分241arc 1)arc 1)1dx tg tg C x ⎡⎤=+++-+⎣⎦+⎰ 对于上述这些情况,都要求建立定积分的近似计算方法—-数值积分法。
1。
1 数值求积分的基本思想根据以上所述,数值求积公式应该避免用原函数表示，而由被积函数的值决定.由积分中值定理：对()[,]f x C a b ∈，存在[,]a b ξ∈，有()()()baf x dx b a f ξ=-⎰表明,定积分所表示的曲边梯形的面积等于底为b a -而高为()f ξ的矩形面积（图4-1）。
问题在于点ξ的具体位置一般是不知道的，因而难以准确算出()f ξ。
我们将()f ξ称为区间[,]a b 上的平均高度。
这样，只要对平均高度()f ξ提供一种算法，相应地便获得一种数值求积分方法.如果我们用两端的算术平均作为平均高度()f ξ的近似值，这样导出的求积公式[()()]2b aT f a f b -=+ (4—1） 便是我们所熟悉的梯形公式（图4-2）。
《数值计算方法》教学大纲
河北联合大学第2012-2013-1学期《数值计算方法》教学大纲依据我校章程，特制定了适合我校理工科各专业本科生的《数值计算方法》教学大纲。
一、课程计划课程名称：数值计算方法Numerical Calculation Methods开课单位：理学院课程类型：专业必修课开设学期：第五学期讲授学时：共15周，每周4学时，共60学时学时安排：课堂教学44学时+实验教学16学时适用专业：信科、数学、统计理科专业本科生教学方式：讲授（多媒体为主）+上机考核方式：闭卷40% +上机实验20%+课程报告20% +平时成绩10%学分：4学分与其它课程的联系预修课程：数学分析、高等代数、常微分方程、计算机高级语言等。
后继课程：偏微分方程数值解及其它专业课程。
二、课程介绍数值计算方法也称为数值分析，是研究用计算机求解各种数学问题的数值方法及其理论的一门学科。
随着计算科学与技术的进步和发展，科学计算已经与理论研究、科学实验并列成为进行科学活动的三大基本手段，作为一门综合性的新科学，科学计算已经成为了人们进行科学活动必不可少的科学方法和工具。
数值计算方法是科学计算的核心内容，它既有纯数学高度抽象性与严密科学性的特点，又有应用的广泛性与实际实验的高度技术性的特点，是一门与计算机使用密切结合的实用性很强的数学课程。
主要介绍数值计算的误差、插值法、函数逼近与曲线拟合、线性方程组迭代解法、数值积分与数值微分、非线性方程组解法、矩阵特征值与特征向量数值计算以及常微分方程数值解，并特别加强实验环节的训练以提高学生动手能力。
通过本课程的学习，不仅能使学生初步掌握数值计算方法的基本理论知识，了解算法设计及数学建模思想，而且能使学生具备一定的科学计算能力和分析与解决问题的能力，不仅为学习后继课程打下良好的理论基础，也为将来从事科学计算、计算机应用和科学研究等工作奠定必要的数学基础。
教学与实验教学课堂教学实验教学论文报告机动课内学时课外学时学时数44 16 8 2 60 10三、重点难点课程重点：理解各种常用数值计算方法的数学原理和理论分析过程，掌握各种数值计算方法的示范性上机程序，学会设计数值算法的基本思路、一般原理和各种数值算法的程序实现。
9-3-4s外推法与龙贝格求积公式
h2 = 1
⎧T1, 0 = ( h / 2 )[ f ( a ) + f ( b )], h = b − a ⎪ ⎨ T m , k +1 − T m , k , ⎪T m + 1, k = T m , k + 1 + 4m −1 ⎩
例题2 利用Romberg序列，近似计算 若T1,0=4, T2,0=5,求f (0.5).
1
3
x
cos xdx
使精度达到10-6.
引言
n+1个节点的插值求积公式
（一）定理：
n k =0
∫
b
a
f ( x) dx ≈ ∑ Ak f ( xk )
n+1个节点的插值型求积公式 ∫ ρ ( x ) f ( x ) dx ≈ a 代数精度最高不超过2n+1次。
b
∑A
k =0
n
k
f ( xk )
的
P261
Romberg方法计算
∫ f ( x)dx 步骤如下：
b a
（1）构造复合梯形序列{T1, k},k=0,1,…, 定义
T1,0= T1 =h0[f(a)+f(b)] /2
T1,k = T 2k , hk =
b−a 2k
h0=b-a
T1,1= T2 = (h1/2)[f(a)+f(b)+2f(a+h1)]=1/2[T1,0 +h0 f(a+h1)] , h1=(b-a)/2 T1,2 = T4 =(1/2){T1,1 +h1 [f(a+h2)+ f(a+3h2)] } , h2=(b-a)/4 ················
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
其中x
j
1 2
1 1 x j h a ( j )h 2 2
n 1 n 1 ba T2 n [ f ( a ) 2 f ( x j ) 2 f ( x 1 ) f (b)] j 4n j 1 j 0 2 n 1 ba b a n 1 [ f (a) 2 f ( x j ) f (b)] 2 f ( x 1 ) j 4n 4n j 0 j 1 2
--------(3)
ba ba n 1,2时， h0 b a, h1 hk k 2 2 则由(1)(2)(3)式,有
ba T1 [ f ( a ) f (b)] T0 (0) 2 1 ba 1 T2 T1 f ( a h) T0 (1) 1 T0 (0) h1 f (a h1 ) 2 2 2 2 若n 2k 记Tn T2k T0 (k ) k 0,1,2,
令
--------(7)
即 当然
Cn C2 k 1 T2 (k 1) C 2 n T2 (k )
--------(8)
同样由复合Cotes公式的余项
I C2 n 1 (C2 n Cn ) 63
64 1 64 1 得 I C2 n Cn T2 ( k ) T2 ( k 1) 63 63 63 63
令
64 1 T3 ( k 1) T2 ( k ) T2 ( k 1) 63 63
--------(9)
将 上 述 结 果 综 合 后
ba T0 (0 ) [ f ( a ) f (b )] 2 2 k 1 1 1 ba ba T0 (k ) T0 (k 1) k f (a (2 j 1) k ) 2 2 2 j 0
/* Romberg Integration */
4 2 0 1 xห้องสมุดไป่ตู้
1
dx
已知对于 = 106 须将区间对分 9 次，得到 T512 = 3.14159202
考察 由 I
4T2 n Tn 4 1 T2 n Tn 来计算 I 效果是否好些？ 41 3 3 Romberg 序列
4 1 T8 T4 = 3.141592502 = S 4 3 3 4T2 n Tn 42 S2n Sn Sn 一般有： Cn 2 41 4 1
43 C 2n C n Rn 3 4 1
Romberg
T1 = T0( 0 )
T2 = T0( 1 )
<?
算法： T4 = T0( 2 )
T8 = T0
(3)
S1 = T1( 0 ) S2 = T1
(1) (2)
<?
C1 = T2
(0) (1)
<?
………………
S n S2
k 1
--------(6)
即 当然
Sn S2 k 1 T1 ( k 1)
S 2 n T1 ( k )
--------(6)
因此由复合Simpson公式的余项
I S2 n 1 ( S2 n Sn ) 15
自己 证明
Cn
可得
16 1 16 1 I S 2 n S n T1 ( k ) T1 ( k 1) 15 15 15 15 16 1 T2 ( k 1) T1 ( k ) T1 ( k 1) Cn 15 15
b
--------(1)
如果将 a, b]分割为2n等份，而h (b a) /n不变, 则 [
n 1 n 1 ba T2 n [ f ( a ) 2 f ( x j ) 2 f ( x 1 ) f (b)] j 4n j 1 j 0 2 --------(2)
4 1 b a n 1 1 I ( Tn f ( x j 1 )) 3 Tn 3 2 2n j 0 2 1 4(b a ) n 1 Tn f (xj 1 ) 3 6n j 0 2
n 1 1 ba 4(b a) n 1 I [ ( f (a) f (b) 2 f ( x j )] f (xj 1 ) 3 2n 6n j 0 j 1 2 n 1 n 1 ba ( f ( a) f (b) 2 f ( x j ) 4 f ( x 1 )] j 6n j 1 j 0 2
Romberg算法的代 数精度为m的两倍 Romberg算法的收敛 阶高达m+1的两倍
T0 (0 ) T0 (1) T0 ( 2 ) T0 ( 3) T1 (0 ) T1 (1) T1 ( 2 ) T2 (0 ) T2 (1)
T3 (0 )
龙贝格积分
例：计算
I T2n 1 I Tn 4
外推 加速 公式
以上整个过程称为Romberg算法
其中外推加速公式可简化为
1 Tm ( k 1) m [ 4 m Tm 1 ( k ) Tm 1 ( k 1)] 4 1
--------(9)
并且m可以推广到m 1,2 ,
Romberg算法求解步骤
k 1, 2 ,
一、复合梯形公式的递推化
将定积分I f ( x)dx 的积分区间 a , b]分割为n等份 [ a ba x j a jh , j 0,1,, n 各节点为 h n 复合梯形(Trapz)公式为
n 1 ba Tn [ f ( a) 2 f ( x j ) f (b)] 2n j 1
1 h n1 Tn f ( x 1 ) j 2 2 j 0 2 1 b a n 1 1 Tn f (a ( j 2 )h) 2 2n j 0 1 b a n 1 ba Tn f (a (2 j 1) 2n ) 2 2 n j 0
S4 = T1
C2 = T2
R1 = T3
(0)
§3 Romberg Integration
理查德森外推法 /* Richardson’s extrapolation */ 利用低阶公式产生高精度的结果。
i 与 h 无关
设对于某一 h 0，有公式 T0(h) 近似计算某一未知值 I。由 Taylor展开得到： T0(h) I = 1 h + 2 h2 + 3 h3 + … h 1 (h) 2 (h)2 3 (h) 3 ... 现将 h 对分，得：T0 ( 2 ) I 2 2 2 Q：如何将公式精度由 O(h) 提高到 O(h2) ？ 2T0 ( h ) T0 ( h) 1 3 2 2 I 2 h 3 h 3 ... 21 2 4 2T0 ( h ) T0 ( h) 2 即： T1 ( h) I 1h2 2 h3 ... 2 1 2 2 T1 ( h ) T1 ( h) 2 T2 (h) I 1h3 2h4 ... 22 1 2 m Tm 1 ( h ) Tm 1 ( h) 2 Tm (h) I 1hm1 2hm2 ... 2m 1
4 1 T1 ( k 1) T0 ( k ) T0 ( k 1) 3 3 16 1 T2 ( k 1) T1 ( k ) T1 ( k 1) 15 15 64 1 T3 ( k 1) T2 ( k ) T2 ( k 1) 63 63
k 1, 2 ,
ba hk 1 k 1 2
ba x j a jhk 1 a j k 1 2
x
j
1 2
1 ba ba 1 x j hk 1 a ( j ) k 1 a ( 2 j 1) k 2 2 2 2
因此(1)(2)(3)式可化为如下递推公式 ba [ f ( a ) f (b )] T0 (0 ) 2 2k 1 1 T (k ) 1 T (k 1) h f (a (2 j 1)h ) 0 k k (4)------- 0 2 j 0 2k 1 1 1 ba ba T0 (k 1) k f (a (2 j 1) k ) 2 2 2 j 0
Sn
设n 2
k 1
复合Simpson公式
4 1 4 1 I T2 n Tn T0 ( k ) T0 ( k 1) 3 3 3 3 4 1 T1 ( k 1) T0 ( k ) T0 ( k 1) --------(5) 3 3
引入 T1 (k 1),令
上式称为递推的梯形公式
k 1, 2 ,
思考
递推梯形公式加上一个控制精度,即 可成为自动选取步长的复合梯形公式
二、外推加速公式
由复合梯形公式的余项公式
I T2 n 1 I Tn 4 1 I T2 n (T2 n Tn ) 3 4 1 I T2 n Tn 3 3
可得 由(3)式