09-2龙贝格算法

2013-2014（1）专业课程实践论文
题目：龙贝格算法

一、算法理论
龙贝格求积公式也称为逐次分半加速法。它是在梯形公式、辛卜生公式和
柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。 作为一
种外推算法, 它在不增加计算量的前提下提高了误差的精度。
在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法
进行。这样，前一次分割得到的函数值在分半以后仍可被利用,且易于编程。

具体实现为：对区间,ab，令hab构造梯形值序列2kT。

1
[()()]/2Thfafb
把区间二等分，每个小区间长度为 /2()/2hba，于是

21
/2/2(/2)TThfah

把区间四(2)等分，每个小区间长度为/2()/4hba，于是

42
/2(/2)(/4)(3/4)TThfahfah
.....................

把,ab 2等分，分点()/(2*)ixabai 0,1,2....2ik每个小
区间长度为()/2ba
二、算法框图
开始
输入a, b ,e
b-a=> h
h/2[f(a)+f(b)]=>T

0 => s
a + h/2 => x

X < b?
t/2 + h/2 *s =>T
It1-t2I < ε？
输出 T2

结束

h/2 => h
T2 => T1

S +f(x) =>s
X + h =>x
三、算法程序
#include
#include

using namespace std;
const int MAXRepeat = 100; //最大允许重复
double function(double x)//被积函数,根据自己的需要手工输入
{
double s;
s = 1.0 / (1 + x);//此为被积函数，积哪个就换哪个
return s;
}
void Romberg(double a, double b, double epsion, double f(double x))
{
int m = 1;
int n = 1;
int k;
double h;
double ep;
double p;
double xk;
double s;
double q;
double T[MAXRepeat];

h = b - a;
T[0] = 0.5 * h * (f(a) + f(b));
ep = epsion + 1.0;
while ((ep >= epsion) && (m < MAXRepeat))
{
p = 0.0;
for (k = 0; k < n; k++)
{
xk = a + (k + 0.5) * h; // n-1
p = p + f(xk); //计算∑f(xk+h/2),T
} // k=0
p = (T[0] + h * p) / 2.0; //T`m`(h/2),变步长梯形求积公式
s = 1.0;
for (k = 1; k <= m; k++)
{
s = 4.0 * s;
//[pow(4,m)T`m`(h/2)-T`m`(h)]/[pow(4,m)-1],2m阶牛顿柯斯特公式，即龙
贝格公式
q = (s * p - T[k - 1]) / (s - 1.0);
T[k-1] = p;
p = q;
}
ep = fabs(q - T[m - 1]);
m++;
T[m - 1] = q;
n++; // 2 4 8 16
h /= 2.0;
}

for (int i = 0; i < m; i++)
{
int j;
if (!(i % j))
{
cout< }
else
{
cout< }
j++;
}
}

int main()
{
double a;
double b;
double epsion;
cout<<"Please input the lower limit: ";
cin>>a;
cout<<"Please input the upper limit: ";
cin>>b;
cout<<"Please input the precision : ";
cin>>epsion;
Romberg( a, b, epsion, function);
return 0;
}
四、算法实现
例1．
求积分101/(1)xdx在精度为0.01下的积分值。

解：把上页代码编译运行后，输入被积函数的上下线和被要求精度，即可求得
如下结果

例2．
求积分3221/(1)xdx在精度为0.02下的积分值。

解：在代码中标注被积函数的地方把原被积函数换成21/(1)x，然后重新编译
运行