龙贝格求积公式
龙贝格求 积分
龙贝格(Romberg ）求积法1.算法理论Romberg 求积方法是以复化梯形公式为基础，应用Richardson 外推法导出的数值求积方法。
由复化梯形公式 )]()(2)([2222b f h a f a f h T +++=可以化为)]()]()([2[212112h a f h b f a f hT +++=＝)]([21211h a f h T ++一般地，把区间［a,b ］逐次分半k －1次，（k ＝1，2,……,n)区间长度（步长)为kk m a b h -=，其中mk ＝2k －1。
记k T ＝)1(k T由)1(k T ＝]))12(([21211)1(1∑=---++km j k k k h j a f h T 从而⎰badxx f )(＝)1(kT－)(''122k f h a b ξ- (1)按Richardson 外推思想，可将（1）看成关于k h ，误差为)(2k h O 的一个近似公式，因而，复化梯形公式的误差公式为⎰badxx f )(－)1(k T ＝......4221++kkh K h K ＝∑∞=12i i k i h K （2）取1+k h ＝k h 21有 ⎰ba dx x f )(－)1(1+k T ＝∑∞=+121221i ik ii hK （3）误差为)(2jh O 的误差公式 )(j kT＝)1(-j kT＋141)1(1)1(------j j k j k T T2。
误差及收敛性分析（1）误差，对复化梯形公式误差估计时，是估计出每个子区间上的误差,然后将n 个子区间上的误差相加作为整个积分区间上的误差。
(2）收敛性,记h x i =∆，由于∑=++=ni i i n x f x f h f T 01))]()([2)(＝))()((21101∑∑-==∆+∆n i ni i i i i x x f x x f上面两个累加式都是积分和，由于)(x f 在区间],[b a 上可积可知，只要],[b a 的分划的最大子区间的长度0→λ时,也即∞→n 时,它们的极限都等于积分值)(f I 。
Romberg求积公式
)
T3
(
h 2k
)
T4
(
h 2k
)
0 (1)T1 h
(3)T2 h
(6)T3 h (10)T4 h
1
(2)
T1
(
h 2
)
(5)
T2
(
h 2
)
(9)
T3
(
h 2
)
(14)
T4
(
h 2
)
2
(4)
T1
(
h 4
)
(8)
T2
(
h 4
)
(13)
T3
(
h 4
)
3
(7)
T1
(
h 8
)
(12)
T2
(
h 8
)
4
(11)
数值计算
Romberg求积公式
1.1外推法基本思想
以较小的计算量为代价,达到提高数值结 果的精度是外推法的中心思想.
设f (x) C 2k2[x h , x h ],用中心差商 22
f (x h) f (x h)
G(h)
2
2
h
逼近f ' (x)。
由Taylor展式起截断误差为
f '(x) G(h) 2h2 4h4 ... 2kh2k E(h)
4E1(h) E1(h) 4 1
,
r2 j
1 (4( j1) 3
1)2 j
( j 2,3,...k)。
此时f ' (x) G1(h) O(h4 )
同
理可
以有G1
(h)和G1
(
h 2
)构造G2
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
龙贝格求积公式
S2k
4T2k1 T2k 41
C2k
42 S2k1 S2k 42 1
R2k
4 C 3 2k 1
C2k
43 1
Romberg积分法的一般公式
Tm, j
4
T j 1 m,
j 1
Tm 1,
j 1
4 j1 1
其中
j 2,3, 4;m j
➢ Romberg积分思想
由上节分析知，用复化梯形公式计算积分值 I
T2n 的误差大约为：
1 3 (T2n Tn )
令
I
T2n
1 3 (T2n
Tn )
4T2n Tn 3
由复化梯形公式知
1
b a n1
T2n 2 Tn
2n
k0
f
(
x
k
1
)
2
4T2n Tn
3
1 3 Tn
2(b a) n1 3n k0
Tm,1 T2m1 (m 1) Tm,3 C2m3 (m 3)
Tm,2 S2m2 (m 2) Tm,4 R2m4 (m 4)
Romberg积分表
1
b a n1
T1 , 1
T2n 2 Tn
2n
k0
f
(
x
k
1
)
2
T2 , 1 T2 , 2
1 ba ba
T2,1 2 T1,1
=3.131176471 S2=1/3(4T4-T2)=3.141568628 C1=1/15(16S2-S1)=3.142117648 计算f(1/8) f(3/8) f(5/8) f(7/8)进而求得 T8=1/2{T4+1/4[f(1/8)+f(3/8)+f(5/8)+f(7/8)]}
龙贝格算法
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 的高阶导数，由于
6b复合求积公式龙贝格算法
步长折半：[xi , xi+1/2] ， [xi +1/2 , xi+1]
n1
xi xi +1/2 xi +1
h T2 n f ( xi ) f ( xi 1 2 ) f ( xi 1 2 ) f ( xi 1 ) i 0 4 n1 h f ( xi ) 2 f ( xi 1 2 ) f ( xi 1 ) i 0 4 h n1 h n1 1 h n1 f ( xi ) f ( xi 1 ) f ( xi 1 2 ) Tn f ( xi 1 2 ) 4 i 0 2 i 0 2 2 i 0 13
1 I T2 n (T2 n Tn ) 3
I Tn 4( I T2n )
3I 4T2n Tn
1 3
验后误差估计式 I T2 n (T2 n Tn )
当
T2n Tn 时，T2n即为所求的近似值。
1 (T2 n Tn ) 3
是T2n 的修正项，它与T2n 之和比T2n 、 Tn更接近与真值，即它是一种补偿。
|| T T2-T|< 2-T1|<
输出T2
16
举例
计算精度满足 | T2n Tn | 107
I [ f ]=0.946083070367
例：用梯形法的递推公式计算定积分 解：
1
0
sin( x ) dx , 要求 x
k
0 1 2 3 4 5 6 7 8 9 10
T (k)
梯形法递推公式
1 h n1 1 h n1 T2 n Tn f ( xi 1 2 ) Tn f ( a ih 0.5h) 2 2 i 0 2 2 i 0
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)。
// 龙贝格积分公式.cpp : 定义控制台应用程序的入口点。
//
#include"stdafx.h"
#include<iostream>
#include<math.h>
#define F(x) (exp(x)*cos (x)) //函数举例。
using namespace std;
//------------步长及4，16，64.......的实现----------
double pf (int i)
{
int s=1;
for (int j=0;j<i;j++)
s*=2;
return s/2;
}
//---------定义一个求t1,t2......的函数-------------
double gcc (double a, double b,double aa[][20],int i)
{
double s,h,x;
h=(b-a)/pf(i);
s=0;
x=a+h/2;
do
{
s+=F(x);
x+=h;
}
while (x<b);
aa[0][i]=aa[0][i-1]/2+h/2*s;
return 0;
}
//----------------------主函数---------------------------
int main()
{
double aa[20][20]={0},e,a,b,h;
int j,i,n;
cout <<"请输入积分区间：\na= ";
cin >>a;
cout <<"b= ";
cin >>b;
cout <<"请输入精度：e=";
cin >>e;
aa[0][0]=(b-a)*(F(a)+F(b))/2.0;
gcc(a,b,aa,1);
aa[1][0]=(4*aa[0][1]-aa[0][0])/3;
for (i=1;i<20;i++)
{
gcc(a,b,aa,i+1);//求下一个要用的t。
for ( n=i,j=1; n>0;n--,j++)//加速公式的实现。
aa[j][n]=(pf(2*j+1)*aa[j-1][n+1]-aa[j-1][n])/(pf(j*2+1)-1);
if (fabs(aa[i][1]-aa[i][0])<e)//判断是否达到精度。
break;
else
aa[i+1][0]=(pf(2*i+3)*aa[i][1]-aa[i][0])/(pf(2*i+3)-1);
}
cout <<"积分结果为："<<aa[i][1];
return 0;
}。