Romberg算法
龙贝格算法
龙贝格积分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
5.3龙贝格(Romberg)算法
一、复合梯形公式的递推化
将定积分I = ∫ f ( x )dx的积分区间[ a , b ]分割为n等份 a b−a xk = a + jh , j = 0 ,1,L , n 各节点为 h= n 复合梯形(Trapz)公式为
n−1 b−a Tn = [ f ( a ) + 2 ∑ f ( x j ) + f (b )] 2n j =1
k −3
if(fabs(t3[k-3]-t2[k-2])<0.000000001) // 若 | R2 { printf("k=%3d,I=%12.9f\n",k,t3[k-3]); break; } } } double f(double x) { double y; y=sqrt(1.0+x*x); return y; }
13
其中外推加速公式可简化为
1 Tm ( k − 1) = m [ 4 m Tm − 1 ( k ) − Tm − 1 ( k − 1)] 4 −1
--------(11)
并且m可以推广到 m = 1,2 ,L
Romberg算法求解步骤
k = 1 , 2 ,L
Romberg算法的代 数精度为m的两倍 Romberg算法的收敛 阶高达m+1的两倍
第6章 数值积分与数值微分
龙贝格(Romberg)算法 算法 龙贝格
1
龙贝格(Romberg)算法 算法 龙贝格
综合前几节的内容,我们知道 梯形公式,Simpson公式,Cotes公式的代数精度分别为 1次,3次和5次 复合梯形、复合Simpson、复合Cotes公式的收敛阶分别为 2阶、4阶和6阶 无论从代数精度还是收敛速度,复合梯形公式都是较差的 有没有办法改善梯形公式呢?
7.3Romberg积分(共64张PPT)
第五页，共六十四页。
表7-1
O(h)
O(h2 )
O(h3 )
O(h4 )
例1 设
带余项的差分(chà fēn)公式为
第六页，共六十四页。
导出具有误差为 解令
用 h/2代替(dàitì)h,得
(11) 的外推公式.
h2
(12)
为消去含 h2的项,用4乘(12)式减去 (11)式,得
第七页，共六十四页。
第二页，共六十四页。
把(1)式改写(gǎixiě)为 (3)
用h/2代替(3)式中的h,得
(4) 用2乘(4)式再减去(3)式,消去含h的项,得
(5)
令
,且记
第三页，共六十四页。
那么(5)式可写为
(6) 这里, 逼近 的误差为
再用 h/2 代替(dàitì) h , 使(6)式变为
(7) 用4乘(7)式减去(6)式,消去含 的项,得
则所谓 gn (x)与 的根相同,即是指这两个正
交多项式的根有如下的关系.
第三十二页，共六十四页。
性质5 (1) 区间[a,b]上带权函数 (的x) 正交
多项式序列
与
对应元素之间
只相差一个比例常数.
(2)区间[a,b]上带权函数
首项 系数 (x)
(shǒu xiànɡ)
为1的
正交多项式序列
唯一.
常见的正交多项式有Legendre(勒让德)多
7.5.1 引言 求积公式
当求积系数
、求积节点
(1) 都可以
自由选取时,其代数精确度最高可以达到多少 次?
下面的引理可以回答上述(shàngshù)问题.
第二十四页，共六十四页。
引理1 当求积系数
Romberg龙贝格算法实验报告
Romberg龙贝格算法实验报告课程实验报告课程名称：专业班级： CS1306班学号： Uxx14967 姓名：段沛云指导教师：报告日期：计算机科学与技术学院目录1 实验目的 ........................................................12 实验原理 ........................................................13 算法设计与流程框图 (2)4 源程序 ........................................................ .. 45 程序运行 ........................................................76 结果分析 ........................................................77 实验体会 ........................................................71 实验目的掌握Romberg公式的用法，适用范围及精度，熟悉Romberg算法的流程，并能够设计算法计算积分31得到结果并输出。
1x2 实验原理2.1 取k=0,h=b-a,求T0=数)。
2.2 求梯形值T0(b-a),即按递推公式（4.1）计算T0。
k2h[f(a)+f(b)]，令1→k,(k记区间[a，b]的二分次22.3 求加速值，按公式(4.12)逐个求出T表的第k行其余各元素Tj(k-j)(j=1,2,….k)。
2.4 若|Tk+1-Tk|n-111T2n=[Tn+hn∑f(xi+)]22i=01Sn=T2n+(T2n-Tn)31Cn=S2n+(S2n-Sn)151Rn=C2n+(C2n-Cn)633 算法设计与流程框图算法设计：(先假定所求积分二分最大次数次数为20) 3.1 先求T[k][0] 3.2 再由公式T(k)m4m(k+1)1)=mTm-1-mTm(k-1(k=1,2,) 求T[i][j] 4-14-13.3 在求出的同时比较T[k][k]与T[k-1][k-1]的大小，如果二者之差的绝对值小于1e-5，就停止求T[k][k]；此时的k就是所求的二分次数，而此时的T[k][k]就是最终的结果 3.4 打印出所有的T[i][j]；程序流程图4 源程序#include #include #include #include int main(void) {float f(float(x)) {float y; y=1/x; return y; }floata,b,e,h,s,k,x,T1=0,T2=0,S1=0,S2=0,C1=0,C2=0,R1=0,R2=0; inti=0;printf("请输入积分下限 : "); scanf("%f",&a);printf("\n请输入积分上限 :"); scanf("%f",&b);printf("\n请输入允许误差 :"); scanf("%f",&e); k大学网=1; h=b-a;T1=h*(f(a)+f(b))/2;printf("____________________________________________\n"); printf("计算结果如下 : \n");printf("\nk T2 S2 C2 R2\n");printf("%d %10.7f %10.7f %10.7f %10.7f\n",i,T1,S1,C1,R1); do {x=a+h/2; s=0; while(x{ s=s+f(x); x=x+h; }T2=(T1+s*h)/2; S2=T2+(T2-T1)/3; if(k==1) {T1=T2; S1=S2; h=h/2; k=k+1; }else if(k==2) {C2=S2+(S2-S1)/15; C1=C2; T1=T2; S1=S2; h=h/2; k=k+1; }else if(k==3) {R2=C2+(C2-C1)/63; C2=S2+(S2-S1)/15; C1=C2; T1=T2; S1=S2; h=h/2; k=k+1; } else {C2=S2+(S2-S1)/15;R2=C2+(C2-C1)/63; if(fabs(R2-R1)printf("%d %10.7f %10.7f %10.7f %10.7f\n",i+1,T2,S2,C2,R2); break;} else { R1=R2; C1=C2; T1=T2; S1=S2; h=h/2; k=k+1; } } i++; printf("%d %10.7f %10.7f %10.7f %10.7f\n",i,T2,S2,C2,R2); } while(1); system("pause"); return 0; }5 程序运行6 结果分析如上所示的结果与课本中求得的结果完全一样，表明程序编写正确，且符合要求，事实上，只要再将所求值的精度设置得更小，则所求的结果将更加准确，最终将无限接近于标准值，由上表也可以看出用龙贝格积分法求函数的积分值在精度比较低的情况下就能求到很准确的值！7 实验体会本次实验较为简单，主要时间是耗费在循环判断上面，因为书上已经给了流程图，都是基本的C语言，难度不大。
数学实验题目2 Romberg积分法
数学实验题目2 Romberg 积分法摘要考虑积分()()b aI f f x dx =⎰欲求其近似值，可以采用如下公式：（复化）梯形公式 110[()()]2n i i i hT f x f x -+==+∑ 2()12b a E h f η-''=- [,]a b η∈ （复化）辛卜生公式 11102[()4()()]6n i i i i hS f x f x f x -++==++∑4(4)()1802b a h E f η-⎛⎫=- ⎪⎝⎭ [,]a b η∈ （复化）柯特斯公式 111042[7()32()12()90n i i i i hC f x f x f x -++==+++∑31432()7()]i i f xf x +++6(6)2()()9454b a h E f η-⎛⎫=- ⎪⎝⎭[,]a b η∈ 这里，梯形公式显得算法简单，具有如下递推关系121021()22n n n i i h T T f x -+==+∑因此，很容易实现从低阶的计算结果推算出高阶的近似值，而只需要花费较少的附加函数计算。
但是，由于梯形公式收敛阶较低，收敛速度缓慢。
所以，如何提高收敛速度，自然是人们极为关心的课题。
为此，记0,k T 为将区间[,]a b 进行2k等份的复化梯形积分结果，1,k T 为将区间[,]a b 进行2k等份的复化辛卜生积分结果，2,k T 为将区间[,]a b 进行2k等份的复化柯特斯积分结果。
根据李查逊（Richardson ）外推加速方法，可得到1,11,,0,1,2,40,1,2,41m m k m km k m k T T T m -+-=-⎛⎫=⎪=-⎝⎭可以证明，如果()f x 充分光滑，则有,lim ()m k k T I f →∞= （m 固定）,0lim ()m m T I f →∞=这是一个收敛速度更快的一个数值求积公式，我们称为龙贝格积分法。
最小二乘法和Romberg求积分
计算方法实验报告湖北师范学院计算机科学与技术学院一、实验题目最小二乘法和Romberg求积分二、实验内容1、最小二乘法拟合本实验是对于已知的m+1对不带权系数{}mi iw=的离散数据{},mi i ix y=，计00,min()max()i m i ma bi ix y≤≤≤≤==。
在连续的函数空间[]()(),0C a b f a f b⨯<中选定n+1个线性无关的基函数(){}nk kxσ=来构造拟合函数()()nk k kkx a xσσ==∑，求得当系数是最优化模型()()()201min I,,,,mn k i iia a a x y xσσ=⋅⋅⋅=-⎡⎤⎣⎦∑时的解***01,,,na a a⋅⋅⋅。
2、romberg求积分法Romberg算法是利用复合梯形公式，在对积分区间的步长逐次折半的过程中，求得积分I=∫a b f(x)dx的近似值序列{T2k}。
三、法思想描述1、最小二乘法拟合（1）、根据最小二乘法的拟合原理将问题简化为求由待拟合离散数据的正则方程组的求解问题；（2）、输入一散数据{},mi i ix y=，存放于数组[],[]X m Y m中，并给定拟合次数n根据待拟合的次数确定对应的系数矩阵为范德蒙行列式的线性方程组；（3）、由[][][][]TA m n A m n⨯得到正则方程组的系数矩阵[][]F n n，如下式，同时由[][][]TA m n Y m⨯得到值向量[]B n；(4)、用最列主元素法解出正则方程组的根，就得到“二”中描述的系数的最优解。
2、romberg求积分法Romberg算法是区间主次二分的过程中，对梯形值进行加权平均以获得准确程度较高a0a1a2…a n1x0x0 2…x0n1 x1x i n…x i n……………………1 x m x m2…x m n＝Y1Y2LLLY m的积分值的一种方法。
对于定积分I= 的复化梯形公式，其余项将积分区间[a，b]，逐次折半，假设,以保证复化梯形公式余项系数是非零的，则构成相应的外推法称为Romberg算法：下标m外推得到的第m个算法。
龙贝格(Romberg)算法的应用实验报告
Lab4 龙贝格（Romberg）算法的应用下面图1中的塑料雨蓬材料是由图2中所示的长方形平板塑料材料压制而成。
图1 图2已知图1的横截面曲线形状满足函数，则给定了雨蓬的长度后，要求需要平板原材料的长度。
函数接口定义：double Integral(double a, double b, double (*f)(double x, double y, double z), double TOL, double l, double t)在接口定义中：a、b分别为定积分的上、下界，f是积分核函数，其中x是积分哑元，y、z是本题目定义的特殊参数，分别对应中的l和t；TOL是要求积分达到的精度；l和t传入裁判输入的参数l和t的值。
另注意：的单位是厘米，输出的长度值要求以米为单位。
裁判程序样例如下：#include<stdio.h>#include<math.h>double f0( double x, double l, double t ){ /* 弧长积分的核函数*/return sqrt(1.0+l*l*t*t*cos(t*x)*cos(t*x));}double Integral(double a, double b, double (*f)(double x, double y, double z), double TOL, double l, double t);int main(){double a=0.0, b, TOL=0.005, l, t;while (scanf("%lf %lf %lf", &l, &b, &t) != EOF)printf("%.2f\n", Integral(a, b, f0, TOL, l, t));return 0;}裁判输入样例：2 100 1标准输出样例：1.68实验报告：1.求解步骤参照书上的龙贝格求积算法2.步骤①利用k=0;h=b-a;T[0][0]=(h/2)*(f(a,l,t)+f(b,l,t));求解T0(0)3.求解T0(0)到T0(k)的值由公式h=b−an 及h=b−a2k可得n=2k其中h为步长，n为二分次数又由递推公式：T2n=12T n+ℎ2∑f(xk+12)n−1k=0得T2k+1=12T2k+b−a2k+1∑f[a+b−a2k+1(2i−1)]2k−1i=1，k=0,1,2,3~其中xk+12= a+ℎ2(2i−1)。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Richardson 外推算法
Romberg 算法
记: T
(k) 0
T2k , T
(k) 1
S2k , T
(k) 2
C2k , T
(k) 3
R2k
T0( k ) : k 次等分后梯形公式计算所得的近似值 ( Tmk ) : m 次加速后所得的近似值
(0)
(1) (2) (3)
( Tmk ) ( ( 4 mTmk 1) Tmk )1 1 4m 1
ba h n
h0 b a
h1 ba 2
n=1
n=2 n=4
1 h1 1 T4 T2 f (a ih1 0.5h1 ) 2 2 i 0 1 h2 3 T8 T4 f ( a ih2 0.5h2 ) 2 2 i 0
ba h2 4
梯形法递推公式
Cha4.2.m
程序演示
Cha41.m: 复化N-C公式 Cha42.m: 梯形法的递推
龙贝格 (Romberg) 加速
梯形法递推公式算法简单，编程方便
但收敛速度较 慢 解决方法：龙贝格 (Romberg) 加速
复化积分公式的渐近状态
思想：利用余项公式与积分的定义
梯形公式：
I [ f ]=0.946083070367
1
0
sin( x ) dx x
T4 0.944513522
k
0 1 2 3 4 5 6 7 8 9 10
T2 0.939793285
T0(k)
0.920735492 0.939793285 0.944513522 0.945690864 0.945985030 0.946058561 0.946076943 0.946081539 0.946082687 0.946082975 0.946083046
2h 6 I Cn [ f (5) (b) f (5) (a)] 945 46
Romberg加速
梯形公式的加速
I T2 n 1 I Tn 4
辛普森公式
辛普森公式的加速
I S2n 1 I S n 16 I C2 n 1 I Cn 64
I
4 1 4 1 而 S n T2 n Tn I T2 n Tn 3 3 3 3 T2 n Tn 事后估计法 I T2 n 3
数值分析及计算软件
Chap 4 数值积分与数值微分
4.1 Romberg 算法
Romberg 算法
利用复化梯形公式、复化Simpson公式、复化 Cotes公式等计算定积分时，如何选取步长 h 太 大: 太 小: 计算精度难以保证
？
计算量增大
解决办法：变步长算法 通常采取将区间不断对分的方法，即取 n = 2k ， 反复采用复化求积公式，直到所得到的计算结果 满足指定的精度为止。
① T1 =T0
② T2 =T0
③ S1 =T1
(0) (1) (2)
Romberg 算法是收敛的
⑥ C1 =T2
⑨ C2 =T2
(0) (1) (0)
④ T4 =T0
⑦ T8 =T0
⑤ S2 =T1
⑧ S4 =T1
⑩ R1 =T3
举例
例：计算定积分
T1 0.920735492
梯形法的余项展开式的推导
T (h) I[ f ] 1h2 2h4 i h2i
必要性: Richardson外推法的基础
推导方法：Taylor展开方法
作业
1.证明：梯形公式的Romberg加速为辛普森公式 2. 教材 P 104/11
hk 1 ba 2 k 1
举例
计算精度满足 | T2n Tn | 107
I [ f ]=0.946083070367
例：用梯形法的递推公式计算定积分 解：
T (0)
1
0
sin( x ) dx , 要求 x
ba f ( a) f ( b) 0.920735492 2 1 (0) h0 0 T (1) T f (a ih0 0.5h0 ) 0.939793285 2 2 i 0 1 (1) h1 1 ( 2) T T f (a ih1 0.5h1 ) 0.944513522 2 2 i 0 1 ( 2) h2 1 ( 3) T T f ( a ih2 0.5h2 ) 0.945690864 2 2 i 0
定义
如果一种复化求积公式 I n 满足下列关系
I In lim p C , C 0 h 0 h
则称该求积公式是 p 阶收敛的. 有如下误差估计式
h2 I Tn [ f ' (b) f ' (a)] 12 h4 I Sn [ f ' ' ' (b) f ' ' ' (a)] 4 180 2
T (h) I[ f ] 1h2 2h4 i h2i
ba h n
其中系数
k , k 1,2,3,...
与h无关
梯形法的加速
T (h) I[ f ] 1h2 2h4 i h2i I[ f ] O( h2 )
梯形法递推公式
将 [a, b] 分成 n 等分 [xi , xi+1] ，xi a i h,
h ba n n1 n1 h f ( x ) f ( x ) h f (a) 2 f ( x ) f (b) Tn i i i 1 2 2 i 0 i 1
步长折半：[xi , xi+1/2] ， [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
1 S1 4T2 T1 3 0.946145882
1 S2 4T4 T2 3 0.946086934
1 C1 16 S2 S1 0.946083004 15
Cha43.m
举例
例：用 Romberg 算法计算定积分
满足
( ( ) | Tmk ) Tmk 1 | 107
I Tn 1 [ f ' (b) f ' (a)] 2 h 12
辛普森公式：
I Sn 1 [ f ' ' ' (b) f ' ' ' (a )] 4 4 h 180 2
Cotes公式：
I Cn 2 [ f (5) (b) f ( 5) (a)] h6 945 46
h h h h T I [ f ] 1 2 i 2 2 2 2
2 4 2i
4T ( h/2) T (h) 3I[ f ] ( 3 / 4)2h4 ( 15 ห้องสมุดไป่ตู้ 16)3h6
I [ f ]=0.4
1
0
x 3 dx , 要求计算精度
Cha44.m
解：逐步计算可得
k
0 1 2
T0(k)
0.50000000
T1(k)
T2(k)
T3(k)
T4(k)
T5(k)
0.42677670 0.40236893 0.40701811 0.40043192 0.40030278
3 3 3
0.40181246 0.40007725 0.40005361 0.40004965 0.40046340 0.40001371 0.40000948 0.40000878 0.40000862 0.40011767 0.40000243 0.40000168 0.40000155 0.40000152 0.40000152
hk 1 1 T2k T2k 1 2 2
2k 1 1
i 0
f ( a ihk 1 0.5hk 1 )
hk 1
ba k 1 2
k 1
记 T
(k)
T2k
2k 1 1
n2
T
(k)
1 ( k 1) hk 1 T 2 2
i 0
f (a ihk 1 0.5hk 1 )
n1
梯形法递推公式
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
ba T1 f (a) f (b) 2 h0 0 1 T2 T1 f (a ih0 0.5h0 ) 2 2 i 0
Cotes公式