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 =+-=-应当是更好的结果。

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 积分法摘要考虑积分()()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求积分二、实验内容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个算法。

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)。

利用Matlab 实现Romberg 数值积分算法一、内容摘要针对于某些多项式积分，利用Newton —Leibniz 积分公式求解时有困难，可以采用数值积分的方法，求解指定精度的近似解，本文利用Matlab 中的.m 文件编写了复化梯形公式与Romberg 的数值积分算法的程序，求解多项式的数值积分，比较两者的收敛速度。

二、数值积分公式1.复化梯形公式求解数值积分的基础是将区间一等分时的Newton —Cotes求积公式:I =(x)[f(a)f(b)]2bab af dx -≈+⎰ 其几何意义是，利用区间端点的函数值、与端点构成的梯形面积来近似(x)f 在区间[a,b]上的积分值，截断误差为：3"(b a)()12f η-- (a,b)η∈ 具有一次的代数精度，很明显，这样的近似求解精度很难满足计算的要求，因而，可以采用将积分区间不停地对分，当区间足够小的时候，利用梯形公式求解每一个小区间的积分近似值，然后将所有的区间加起来，作为被求函数的积分，可以根据计算精度的要求，划分对分的区间个数，得到复化梯形公式：I =11(b a)(b a)(x)dx [f(a)f(b)2(a )]2n bak k f f n n -=--≈+++∑⎰其截断误差为：2"(b a)h ()12R f η--=(a,b)η∈ 数值积分算法使用复化的梯形公式计算的数值积分，其收敛速度比减慢，为此，采用Romberg 数值积分。

其思想主要是，根据I 的近似值2n T 加上I 与2n T 的近似误差，作为新的I 的近视，反复迭代，求出满足计算精度的近似解。

用2n T 近似I 所产生的误差可用下式进行估算：12221()3n n n I T T T -∆=-=-新的I 的近似值：122n n j T T -=∆+ j =（0 1 2 ….） Romberg 数值积分算法计算顺序i=0 (1) 002Ti=1 (2) 102T(3) 012Ti=2 (4) 202T (5) 112T(6) 022Ti=3(7)302T (8) 212T (9) 122T(10) 032T i=4 (11)402T(12) 312T(13) 222T(14) 132T…………其中，第一列是二阶收敛的，第二列是四阶收敛的，第三列是六阶收敛的，第四列是八阶收敛的，即Romberg 序列。