Romberg(龙贝格)求积公式求解数值积分时的注意事项

《数值分析》第5版（李庆扬编著）的第四章课后习题第8-(2)题中，要求使用Romberg(龙贝格)求积公式求解f(x)=xsinx在区间[0,2pi]上的积分，要求误差小于10^(-5)。

针对此问题，套用计算公式求解即可。在第一步计算梯形公式时，出现了T0=pi*[f(0)+f(2pi)]。很显然，若按照sin(2pi)=sin(pi)=0计算，则Romberg(龙贝格)求积公式T0列前三个数值均会出现0，从而造成计算收敛速度较慢，需要迭代多次。（实际计算迭代次数大于10次时，依旧不能满足误差小于10^(-5)的条件）。

针对此问题，上述计算过程应对pi指定有效位数，化作有理数求解。本题中由于题目条件要求误差小于10^(-5)，故pi可以取6位小数或大于等于6位小数。本文中取7位小数，以pi=3.1415927为例进行计算，此时计算只需迭代5次，即可满足题目条件。下面给出对应的matlab程序

m_pi=3.1415927;   %取7位小数
t0=m_pi*(2*m_pi*sin(2*m_pi));   % T00   k=0
t0=vpa(t0,8)

N=2;    %k=1
j=1;
sum=0;
for i=1:1:N
    sum=sum+((j*m_pi/N)*(sin(j*m_pi/N)));
    j=j+2;
end
sum=sum*m_pi/(N);
t1=t0/2+sum;
t1=vpa(t1,8)   %t01

N=4;  %k=2
j=1;
sum=0;
for i=1:1:N
    sum=sum+((j*m_pi/N)*(sin(j*m_pi/N)));
    j=j+2;
end
sum=sum*m_pi/(N);
t2=t1/2+sum;
t2=vpa(t2,8)    %t02

N=8;  %k=3
j=1;
sum=0;
for i=1:1:N
    sum=sum+((j*m_pi/N)*(sin(j*m_pi/N)));
    j=j+2;
end
sum=sum*m_pi/(N);
t3=t2/2+sum;
t3=vpa(t3,8)    %t03

N=16;  %k=4
j=1;
sum=0;
for i=1:1:N
    sum=sum+((j*m_pi/N)*(sin(j*m_pi/N)));
    j=j+2;
end
sum=sum*m_pi/(N);
t4=t3/2+sum;
t4=vpa(t4,8)    %t04

N=32;  %k=5
j=1;
sum=0;
for i=1:1:N
    sum=sum+((j*m_pi/N)*(sin(j*m_pi/N)));
    j=j+2;
end
sum=sum*m_pi/(N);
t5=t4/2+sum;
t5=vpa(t5,8)    %t05

N=64;  %k=6
j=1;
sum=0;
for i=1:1:N
    sum=sum+((j*m_pi/N)*(sin(j*m_pi/N)));
    j=j+2;
end
sum=sum*m_pi/(N);
t6=t5/2+sum;
t6=vpa(t6,8)    %t06

%%  T1x
a0=(4/3)*(t1)-(1/3)*(t0);
a1=(4/3)*(t2)-(1/3)*(t1);
a2=(4/3)*(t3)-(1/3)*(t2);
a3=(4/3)*(t4)-(1/3)*(t3);
a4=(4/3)*(t5)-(1/3)*(t4);
a5=(4/3)*(t6)-(1/3)*(t5);
a0=vpa(a0,8)
a1=vpa(a1,8)
a2=vpa(a2,8)
a3=vpa(a3,8)
a4=vpa(a4,8)
a5=vpa(a5,8)


%%  T2x
b0=(16/15)*(a1)-(1/15)*(a0);
b1=(16/15)*(a2)-(1/15)*(a1);
b2=(16/15)*(a3)-(1/15)*(a2);
b3=(16/15)*(a4)-(1/15)*(a3);
b4=(16/15)*(a5)-(1/15)*(a4);
b0=vpa(b0,8)
b1=vpa(b1,8)
b2=vpa(b2,8)
b3=vpa(b3,8)
b4=vpa(b4,8)

%% T3x
c0=(64/63)*(b1)-(1/63)*(b0);
c1=(64/63)*(b2)-(1/63)*(b1);
c2=(64/63)*(b3)-(1/63)*(b2);
c3=(64/63)*(b4)-(1/63)*(b3);
c0=vpa(c0,8)
c1=vpa(c1,8)
c2=vpa(c2,8)
c3=vpa(c3,8)

%% T4x
d0=(256/255)*(c1)-(1/255)*(c0);
d1=(256/255)*(c2)-(1/255)*(c1);
d2=(256/255)*(c3)-(1/255)*(c2);
d0=vpa(d0,8)
d1=vpa(d1,8)
d2=vpa(d2,8)

%% T5x
e0=(1024/1023)*(d1)-(1/1023)*(d0);
e1=(1024/1023)*(d2)-(1/1023)*(d1);
e0=vpa(e0,8)
e1=vpa(e1,8)

%% T6x
f0=(4096/4095)*(e1)-(1/4095)*(e0);
f0=vpa(f0,8)

运行上述程序，会发现只需迭代5次，即可满足误差要求，计算结果为-6.2831853，计算表格如下表所示。

k	h	T0	T1	T2	T3	T4	T5
0	2pi	0.0000018322016
1	pi	-4.9348014	-6.5797359
2	pi/2	-5.9568328	-6.29751	-6.2786949
3	pi/4	-6.2022313	-6.2840308	-6.2831322	-6.2832026
4	pi/8	-6.2629859	-6.2832374	-6.2831845	-6.2831853	-6.2831852
5	pi/16	-6.2781379	-6.2831885	-6.2831853	-6.2831853	-6.2831853	-6.2831853

此外，本文给出计算Romberg(龙贝格)求积公式的C++程序，如下所示【标注：该程序数据类型需进一步完善--modified 2015-11-22】。

#include<iostream>
#include<math.h>
using namespace std;

# define Precision 0.00001//积分精度要求
# define e         2.71828183
#define  MAXRepeat 100  //最大允许重复

double function(double x)//被积函数
{
	double s;
	s = x* sin(x);
	return s;
}

double Romberg(double a, double b, double f(double x))
{
	int m, n, k;
	double y[MAXRepeat], h, ep, p, xk, s, q;
	h = b - a;
	y[0] = h*(f(a) + f(b)) / 2.0;//计算T`1`(h)=1/2(b-a)(f(a)+f(b));
	m = 1;
	n = 1;
	ep = Precision + 1;
	int num = 0;
	while ((ep >= Precision) && (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 = (y[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)
			q = (s*p - y[k - 1]) / (s - 1.0);//[pow(4,m)T`m`(h/2)-T`m`(h)]/[pow(4,m)-1],2m阶牛顿柯斯特公式，即龙贝格公式
			y[k - 1] = p;
			p = q;
		}
		ep = fabs(q - y[m - 1]);//前后两步计算结果比较求精度
		m = m + 1;
		y[m - 1] = q;
		n = n + n;   //  2 4 8 16 
		h = h / 2.0;//二倍分割区间
		num++;
	}
	cout <<"迭代次数为"<< num <<"次"<< endl;
	return q;
}
int   main()
{
	double a, b, Result;
	cout << "请输入积分下限:" << endl;
	cin >> a;
	cout << "请输入积分上限:" << endl;
	cin >> b;
	Result = Romberg(a, b, function);
	cout << "龙贝格积分结果：" << Result << endl;

	system("pause");

	return 0;
}

程序运行结果如下图所示。