Romberg求积分公式

龙贝格(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数值积分实验报告一、实验目的1、加深外推法的原理理解, 掌握Romberg外推法的计算方法。

2、用matlab软件实现Romberg数值积分来计算题目的运算。

二、基本理论及背景1、理论推导：对区间[a， b]，令h=b-a构造梯形值序列{T2K}。

T1=h[f(a)+f(b)]/2把区间二等分，每个小区间长度为h/2=(b-a)/2，于是T2 =T1/2+[h/2]f(a+h/2）把区间四(2)等分，每个小区间长度为h/2 =（b-a）/4，于是T4 =T2/2+[h/2][f（a+h/4）+f（a+3h/4).....................把[a，b] 2等分，分点xi=a+（b-a）/ 2 ·i （i =0，1，2 · · · 2k）每个小区间长度为（b-a）/ 2 .2、参考Romberg数值积分，实现积分的数值求解，完成下列题目：三、算法设计及实现1、算法设计(a)function y=fun1(x)y=sqrt(4-(sin(x)).^2)；(b)function y=fun2 (x)y=log(1+2*x)/3;(c)function y=fun3(x)y=exp(-x).*sin(x.^2)四、实验步骤1、○1打开matlab软件，新建Romberg.m文件,在窗口中编辑Romber数值积分函数程序代码,并保存在指定的文件夹下,在Current Directory窗口右边点击《Browse For Folder》按钮指向Romberg.m文件；○2在Command Window中编辑相应要计算的题目的数值函数及相应的题目的表达式。

2、输出结果和初步分析说明（见附件一）。

五、使用说明实验结果分析1、在Command Window窗口中编辑要调用的函数名与指定的函数名字不同导致出现错误，通过改正与函数名相同即可。

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 ）求积法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++k k h 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 。

《MATLAB程序设计实践》课程考核1、编程实现以下科学计算算法，并举一例应用之。

“Romberg求积分公式”2、编程解决以下科学计算和工程实际问题。

1）、给定半径的为r，重量为Q的均质圆柱，轴心的初始速度为v0，初始角速度为w0且v0>r*w0，地面的摩擦系数为f，问经过多少时间后，圆柱将无滑动地滚动，求此时的圆柱轴心的速度。

2）、在一丘陵地带测量高程，x和y方向每隔100m测一个点，得高程数据如下，试拟合一曲面确定合适的模型，并由此找出最高点和该点的高程。

一、Romberg求积分公式1、算法说明：此算法可自动改变积分步长,使其相临两个值的绝对误差或相对误差小于预先设定的允许误差.Romberg加速法公式在等距节点的情况下，通过对求积区间（a，b）的逐次分半，由梯形公式出可逐次提高求积公式精度，这就是Romberg求积的基本思路，由于梯形公式余项只有精度，即，但当节点加密时可组合成其精度达到，如果再由与组合成则可使误差精度达到，于是依赖于x，若在上各阶导数存在，将展开，可将展成的幂级数形式，即，记的计算精度，可利用外推原理逐次消去式右端只要将步长h逐次分半，利用及组合消去，重复同一过程最后可得到递推公式，此时.说明用其误差阶为，这里表示m次加速。

计算时用序列表示区间分半次数，即具体计算公式为，就是Romberg求积方法。

2、程序代码：M文件1）、Romberg加速法function [s,n]=rbg1(a,b,eps)if nargin<3,eps=1e-6;ends=10;s0=0;k=2;t(1,1)=(b-a)*(f(a)+f(b))/2;while (abs(s-s0)>eps)h=(b-a)/2^(k-1);w=0;if (h~=0)for i=1:(2^(k-1)-1)w=w+f(a+i*h);endt(k,1)=h*(f(a)/2+w+f(b)/2);for l=2: kfor i=1:(k-l+1)t(i,l)=(4^(l-1)*t(i+1,l-1)-t(i,l-1))/(4^(l-1)-1);endends=t(1,k);s0=(t(1,k-1));k=k+1;n=k;else s=s0;n=-k;endend2）、改进的Romberg求积函数function [s,eer]=rbg2(a,b,eps)if nargin<3,eps=1e-6;endm=1;t(1,1)=(b-a)*(f(a)+f(b))/2;r(1,1)=0;while ((abs(t(1,m)-r(1,m))/2)>eps)c=0;m=m+1;for j=1:2^(m-1)c=c+f(a+(j-0.5)*(b-a)/2^(m-1));endr(m,1)=(b-a)*c/2^(m-1);for j=2:mfor k=1:(m-j+1)r(k,j)=r(k+1,j-1)+(r(k+1,j-1)-r(k,j-1))/(4^(j-1)-1);endt(1,j)=r(1,j-1)+2*(4^(j-2)-1)*(t(1,j-1)-r(1,j-1))/(4^(j-1)-1);endenderr=abs(t(1,m)-r(1,m))/2;s=t(1,m);3）定义f.m函数如下：function f=f(x);f=x.^3;4）运行命令及结果>> rbg1(0,2)>> rbg2(0,2)3、流程图二、圆柱体问题1、问题分析圆柱体水平方向受到地面的摩擦阻力f*Q，该摩擦力对轴心速度起减速作用，同时又产生一个力矩，对角速度起加速作用。

安徽中医药大学题目：FOnberg求积法c语言编程姓名：杨撞撞学号：13713042班级：13医软（1）班目录1简介2计算公式3算法描述4程序流程图5算法程序表示6算法结果截图1. 简介龙贝格求积公式也称为逐次分半加速法。

它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。

作为一种外推算法，它在不增加计算量的前提下提高了误差的精度.在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。

这样，前一次分割得到的函数值在分半以后仍可被利用，且易于编程。

2. 计算公式梯形公式----- 二复化辛普森公式------ A复化科特斯公式V龙贝格求积公式其对应的公式为：T2n=1/2（Tn+Hn）（梯形公式）Sn=4/（4-1）T2n-1（4-1）Tn （辛普森公式）Cn=4八2/（4八2-1）S2n-1/（4八2-1）Sn （柯特斯公式）Rn=4八3/（4八3-1）C2n-1/（4八3-1）Cn （龙贝格求积法公式）3. 算法描述3.1龙贝格算法基本描述先算出T0 （ 0）,从而计算出T0（1）,以此类推，直到计算出|T0（0）-Tn-1 （0）|<e即可利用加速推算公式推算出结果。

3.2龙贝格算法程序包步骤1.输入积分上限2 输入程序下限3输入区间等分数4输入要求的函数5计算出所求函数的积分，分别是：复化梯形求积结果辛普森求积结果柯特斯求积结果龙贝格求积结果4. 程序流程图例题：用Romberg 方法计算积分I= 01sin （x）/xdx 的相关算法流程图表示如下图开始 欢迎界面择函 择/——输入1——输入信息：积分:卜线，等分段积.4.x2*sinxs Ynx 2* cosx Yx *sin 2xYsi nx/xX合法？3 输入0Y*辛普林 求积科特斯 求积龙贝格 求积输出结 果5.算法程序表示#include<stdio.h> #include<math.h> #defi ne A(x)(si n(x)/x)//宏定义若干常用函数A,B,C,D,E,G#define B(x)(cos(x*x+2*x+1))#define C(x)(atan(sqrt(x*x+1)))#define D(x)(sqrt(exp(x)+sin(2*x)))#define E(x)(x*x*x+3*x*x+5)#define G(x)(log10(x)/pow(2,x))double t[20],s[20],c[20],r[20];// 定义全局数组double dh,fan,a,b,m; // 定义全局变量int jj=0;char hs;double F(double x) // 用switch 调用若干被积函数{switch(hs){case 'A':fan=A(x);break;case 'B':fan=B(x);break;case 'C':fan=C(x);break;case 'D':fan=D(x);break;case 'E':fan=E(x);break;case 'G':fan=G(x);break;default :printf(" 输入错误！");}return(fan);// 返回被积函数值} double H(int i) // 求和函数并返回和SUM{int j;double zh,SUM=0.0; II定义求和变量SUM并赋初值for(j=1;j<=pow(2,i-1);j++){ zh=(a+((2*j-1)*(b-a))Ipow(2,i));SUM=SUM+F(zh); II 调用F(x) 函数}SUM=(b-a)*SUMIpow(2,i);return(SUM);}double Txing(int k) II 梯形公式if(k==O) dh=t[jj]=((b-a)/2)*(F(a)+F(b));//分半次数为零时T 形公式求积 else {dh=0.5*Txi ng(k-1)+H(k); //Txing函数递归调用循环输出并返回dht[++jj]=dh; }m=pow(2,jj);prin tf("T[%0.0lf]=%0.7lf\t",m,t[jj]);//输出并返回dhreturn(dh); }double Simpson(int k) // 辛普森公式 错误！未找到引用源int i,j;b —俱 V 1.「 "‘+弓厂『[°十(雄-耳{廉G +甩)],Txin g(k); // 调用梯形公式for(i=0;i<=k-1;i++)s[i]=(4.0*t[i+1]-t[i])/3.0; II 森公式prin tf("\n");for(j=0;j<=k-1;j++) //{m=pow(2,j);prin tf("S[%0.0lf]=%0.7lf",m,s[j]);prin tf("\t");}return(s[k-1]); // 返回最后一个值}double Cotes(i nt k) //{int i,j;Simps on (k);for(i=0;i<=k-2;i++)c[i]=(16.0*s[i+1]-s[i])/15.0; //递推辛普循环输出s[k-1]科特斯公式递推科特斯公式prin tf("\n");for(j=0;jv=k-2;j++) II环输出{m=pow(2,j);prin tf("C[%0.0lf]=%0.7lf\t",m,c[j]);}return(c[k-2]); // 返回最后一个值}double Romberg(i nt k) //公式」-亠厂{int i,j; //斯公式Cotes(k);for(i=0;i<=k-3;i++) r[i]=(64.0*c[i+1]-c[i])/63.0; //隆贝格公式prin tf("\n");c[k-2]隆贝格调用科特递推for(j=0;j<=k-3;j++) // 环输出{m=pow(2,j); printf("R[%0.0lf]=%0.7lf\t",m,r[j]);printf("\n");}return(r[k-3]); // 返回最后一个值r[k-3] }main(){int k;char y;printf(" 请从以下公式中选择积分函数：\n"); printf("A: sin (x)/(x)\tB:cos(xA2+2x+1)\tC:ata n( sqrt(x A2+1))\n\n");tG:log10(x)/pow(2,x)\n\n");printf(" 请选择函数F(x)(大写)：\n");scanf("%c",&hs);。