龙贝格积分公式

龙贝格(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 。

龙贝格积分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第二题
求数值积分精确到。

解：该题可以用复合T型求积公式、复合Simpson求积公式、龙贝格（Romberg）公式来做。

复合T型求积公式、复合Simpson公式要自己设步长。

我采用的是龙贝格公式来算：
2.1用龙贝格公式进行计算
依题意得：
（1）a=-1，b=1;令，
求得：
（2）计算初始值：
（3）将积分区间[a,b]二分一次计算梯形值：
（4）鉴于理论分析表明梯形值T（h）具有如下误差渐进展开式
其中，诸，无关。

将步骤（2）、（3）的梯形值外推一次得逼近值.
（5）再将区间[a,b]二分二次计算梯形值：
(6将、外推一次得逼近值：
（7）进一步将、外推一次得逼近值：
（8）依据公式
得到：
重复上述步骤即可获得系列逼近值：
由于已达到预定精度，故取.。

实验题目2 Romberg 积分法摘要考虑积分()()b aI f f x dx =⎰欲求其近似值，可以采用如下公式： （复化）梯形公式 11[()()]2n ii 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 →∞=这是一个收敛速度更快的一个数值求积公式，我们称为龙贝格积分法。

该方法的计算可按下表进行0,0T 0,1T 0,2T … 0,m T 1,0T 1,1T … 1,1m T - 2,0T … 2,2m T - … … ,0m T很明显，龙贝格计算过程在计算机上实现时，只需开辟一个一维数组，即每次计算的结果,m k T ，可存放在0,k T 位置上，其最终结果,0m T 是存放在0,0T 位置上。

龙贝格（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 。

数值计算方法大作业注:由朱福利邹素云共同完成龙贝格求积公式【简介】龙贝格求积公式也称为逐次分半加速法。

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

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

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

【算法】对区间[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 .例：I = ∫0(4/1+X) dx解按上述五步计算，此处 f(x)=4/(1+x) 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)]}=3.131176471S2=1/3(4T4-T2)=3.141568628C1=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)]}=3.138988495S4=1/3(4T3-T4)=3.141592503C2=1/15(16S4-S2)=3.141594095R1=1/63(64C2-C1)=3.141585784把区间再二分，重复上述步骤算得T16=3.140941613 S8=3.141592652C4=3.141592662 R2=3.141592640程序步骤1.输入积分上线2.输入积分下线3.输入区间等分数4.输入要求的函数5.计算出所求函数的积分，分别是●复化梯形求积结果●辛普森求积结果●科特斯求积结果●龙贝格求积结果流程图源程序// MyT ask.cpp : Defines the entry point for the console application.//#include "stdafx.h"#include "math.h"#include<stdlib.h>#define IS_INT(x) ( (x) == int( (x) * 10 / 10 ))void fuHuaTiXing(int k, int a, int b);//复化梯形求积里面T1跟T2double f(double x);//原函数sinx / xdouble fuHuaTiXing(int k, double y,int a, int b);//复化梯形求积里面T1跟T2double leiJia(int a, int b, int i);//复化梯形求积里面的累加void simpson(int k);//k代表要计算的simpson的个数void cotes(int k);//k >=3void romberg(int k );double TTT[512] = {0.0};double TT[128] = {0.0};double S[64] = {0.0};double C[32] = {0.0};double R[16] = {0.0};int flag;int main(int argc, char* argv[]){int a,b,k;//a 积分下线b 积分上限k 等分数//k最大支持256等分int _k;double result;int loop_k;int g_quit = 0;int g_quit1 = 0;//欢迎界面printf("\t\t**************************************** ************\n");printf("\t\t**************************************** ************\n");printf(" 欢迎使用RomBerg方法计算积分\n");printf(" 本程序支持最小4等分最大256等分\n");printf("\t\t****************************************************\n");printf("\t\t**************************************** ************\n");printf("想继续运行吗，继续请按1，退出请按0\n");scanf("%d",&g_quit);if(1 == g_quit){while(1){printf("please input 3 numbers whitch means three constent.\n");printf("please input the first number whitch means 积分下限\n");getchar();scanf("%d",&a);printf("please input the second number whitch means 积分上限\n");getchar();scanf("%d",&b);printf("please input the third number whitch means 等分数\n\t\t\t\t\t注意:支持4等分到256等分\n \t\t\t\t\t注意:一定是2的整数幂\n");scanf("%d",&k); //if(4 == (int)k){_k = 2;}else if(8 == (int)k) {_k = 3;}else if(16 == (int)k) {_k = 4;}else if(32 == (int)k) {_k = 5;}else if(64 == (int)k) {_k = 6;}else if(128 == (int)k){_k = 7;}else if(256 == (int)k){_k = 8;}else{printf("您输入的等分区间不是2的整数幂,请重新输入!\n\n");getchar();getchar();getchar();continue;}aaaaa:printf("please select the function of f(x). \n");printf("if you want to chose function 1( pow(x,2)*sin(x) ) please input 1;\n");printf("if you want to chose function 2 ( sin(x)*cos(x) ) please input 2;\n");printf("if you want to chose function 3 ( sin(x)*sin(x)*x ) please input 3;\n");printf("if you want to chose function 4 ( sin(x)/x ) please input 4.\n");scanf("%d",&flag);if(flag==1 || flag==2 || flag==3 || flag==4 ){//null}else{printf("You have inputed wrong num! \nPlease try again!\n\n\n");goto aaaaa;}fuHuaTiXing(_k,a,b);simpson(_k);cotes(_k);romberg(_k);printf("output fuHuaTiXing(复化梯形) result.\n");for(loop_k=0;loop_k<=_k;loop_k++){printf("TT%d=%.14f\n",(int)pow(2,loop_k),TT[loop_k]);}printf("output simpson result.\n");for(loop_k=0;loop_k<=_k-1;loop_k++){printf("S%d= %.14f\n",(int)pow(2,loop_k),S[loop_k]);}printf("output cotes result.\n");for(loop_k=0;loop_k<=_k-2;loop_k++){printf("C%d = %.14f \n",(int)pow(2,loop_k),S[loop_k]);}printf("output romberg result.\n");for(loop_k=0;loop_k<=_k-3;loop_k++){printf("R%d= %.14f\n",(int)pow(2,loop_k),R[loop_k]);}printf("\n\n\n");printf("继续计算请按1，退出请按0\n");scanf("%d",&g_quit1);if(0 == g_quit1){exit(0);}}}else{exit(0);}return 0;}double f(double x)//原函数sinx/x {double result_f;switch(flag){case 1:result_f = pow(x,2)*sin(x);return (result_f);break;case 2:result_f = sin(x)*cos(x);return (result_f);break;case 3:result_f = sin(x)*sin(x)*x;return (result_f);break;case 4:result_f = sin(x);if(x != 0){result_f = result_f / x;return (result_f);}else return 1.0;break;default :break;}}void fuHuaTiXing(int k, int a, int b)//复化梯形求积里面T1跟T2{int loop_i;double result_b_a;result_b_a = double(b - a) / 2;for(loop_i = 0; loop_i<=k; loop_i++)//k >=3{int temp;if(loop_i == 0)//k means zhishu{TT[0] = result_b_a * ( f(a) + f(b) );//代表T1}else{temp = pow(2,loop_i);TT[loop_i] = TT[loop_i - 1] / 2 + ( (double)(b - a) / (temp ) ) * leiJia(a,b,loop_i);}}for(loop_i=0;loop_i<10;loop_i++){TTT[loop_i] = TT[loop_i];//把局部数组转存到全局数组中方便下面的访问}}double leiJia(int a, int b, int i)//复化梯形求积里面的累加{int loop_m ;int loop_i;double result_leijia = 0.0;int temp;double temp1,temp2;switch(i){case 1:loop_m = (int)pow(2,(i-1));break;case 2:loop_m = (int)pow(2,(i-1));break;case 3:loop_m = (int)pow(2,(i-1));break;case 4:loop_m = (int)pow(2,(i-1));break;case 5:loop_m = (int)pow(2,(i-1));break;case 6:loop_m = (int)pow(2,(i-1));break;default:break;}for(loop_i = 1; loop_i<=loop_m; loop_i++){temp = pow(2,i);temp1 = (double)(a + (2 * loop_i - 1) * ( (double)(b - a) / temp ) );//temp2 = f( temp1 );result_leijia = result_leijia + f( temp1 ); //第一次temp1=0.5}return (result_leijia);}//求simpsonvoid simpson(int k)//k代表要计算的ximpson的个数{int loop_i = 0, loop_j = 0;for(loop_i =0;loop_i<=k-1;loop_i++)S[loop_i] = ( 4 * TT[loop_i + 1] -TT[loop_i] ) / 3;}void cotes(int k){int loop_i;for(loop_i=0;loop_i<=k-2;loop_i++)C[loop_i] = ( 16 * S[loop_i+1] - S[loop_i] ) / 15;}void romberg(int k ){int loop_i;for(loop_i=0;loop_i<=k-3;loop_i++)R[loop_i] = ( pow(4,3) * C[loop_i +1] - C[loop_i] ) / (pow(4,3) - 1) ;}。

0.8 、1.用Romberg方法求解积分■. xdx，要求误差不超过a = 10解：Romberg.m 文件：function [I, step] = Romberg(f, a, b,EPS)% Romberg.m是用龙贝格公式求积分%f为被积函数% EPS为积分结果精度% a,b为积分区间的上下限% I为积分结果；step为积分的子区间数m = 1 k = 0 Er = 0.1 H =b-aS = zeros(1, 1)S(1, 1) = (H/2) *(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b)) while Er > EPSk = k + 1f1 = 0H = H/2for i = 1:mx = a +H*(2*i-1)f1 = f1 + subs(sym(f),fi ndsym(sym(f)),x)endS(k+1, 1) = S(k, 1)/2 + H*f1m = 2 * mfor n = 1:kS(k+1, n+1) = S(k+1, n) + (S(k+1, n)-S(k, n))/(4^ n-1) end Er = abs(S(k+1, n+1)-S(k, n))endI = S(k+1, k+1) step = k命令：clearclcformat shorta = 0;b = 0.8; EPS = 1e-2;[I, step] = Romberg('x A(1/2)', a, b, EPS)计算结果:Er =0.1000H =0.8000S =0.3578 f1 =H =0.40000 0.45660.4000 f1 =0.6325S =0.35780.4319S =0.35780.4319 Er =0.0988 f1 =0.20000.2000 f1 =0.4472 x =0.6000 f1 =1.2218S =0.3578 00.4319 0.45660.4603 0S =0.3578 00.4319 0.45660.4603 0.4698S =0.3578 0 00.4319 0.4566 00.4603 0.4698 0.4707 Er =0.0141 f1 =H =0.1000 x =0.1000 f1 =0.3162 x =0.3000 f1 =0.86400.5000f1=1.5711x =0.7000f1=2.4077S =0.3578 0 00.4319 0.4566 00.4603 0.4698 0.47070.47090 0 m=8S =0.3578 0 00.4319 0.4566 00.4603 0.4698 0.47070.4709 0.4745S =0.3578 0 00.4319 0.4566 0max X i 1 i 3 i * 3 5x-| 2x 2 X 3 12, X 1 4x 2 2x 3 20,试用Jacobi 迭代法求解此方程， x (0) 0,0,0 T ,2x 3x 2 10x 3 3. (k)X i 10 5 时终止迭代。