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=0t0=vpa(t0,8)N=2;    %k=1j=1;sum=0;for i=1:1:N    sum=sum+((j*m_pi/N)*(sin(j*m_pi/N)));    j=j+2;endsum=sum*m_pi/(N);t1=t0/2+sum;t1=vpa(t1,8)   %t01N=4;  %k=2j=1;sum=0;for i=1:1:N    sum=sum+((j*m_pi/N)*(sin(j*m_pi/N)));    j=j+2;endsum=sum*m_pi/(N);t2=t1/2+sum;t2=vpa(t2,8)    %t02N=8;  %k=3j=1;sum=0;for i=1:1:N    sum=sum+((j*m_pi/N)*(sin(j*m_pi/N)));    j=j+2;endsum=sum*m_pi/(N);t3=t2/2+sum;t3=vpa(t3,8)    %t03N=16;  %k=4j=1;sum=0;for i=1:1:N    sum=sum+((j*m_pi/N)*(sin(j*m_pi/N)));    j=j+2;endsum=sum*m_pi/(N);t4=t3/2+sum;t4=vpa(t4,8)    %t04N=32;  %k=5j=1;sum=0;for i=1:1:N    sum=sum+((j*m_pi/N)*(sin(j*m_pi/N)));    j=j+2;endsum=sum*m_pi/(N);t5=t4/2+sum;t5=vpa(t5,8)    %t05N=64;  %k=6j=1;sum=0;for i=1:1:N    sum=sum+((j*m_pi/N)*(sin(j*m_pi/N)));    j=j+2;endsum=sum*m_pi/(N);t6=t5/2+sum;t6=vpa(t6,8)    %t06%%  T1xa0=(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)%%  T2xb0=(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)%% T3xc0=(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)%% T4xd0=(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)%% T5xe0=(1024/1023)*(d1)-(1/1023)*(d0);e1=(1024/1023)*(d2)-(1/1023)*(d1);e0=vpa(e0,8)e1=vpa(e1,8)%% T6xf0=(4096/4095)*(e1)-(1/4095)*(e0);f0=vpa(f0,8)

运行上述程序，会发现只需迭代5次，即可满足误差要求，计算结果为-6.2831853，计算表格如下表所示。khT0T1T2T3T4T502pi0.0000018322016     1pi-4.9348014-6.5797359    2pi/2-5.9568328-6.29751-6.2786949   3pi/4-6.2022313-6.2840308-6.2831322-6.2832026  4pi/8-6.2629859-6.2832374-6.2831845-6.2831853-6.2831852 5pi/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-1p = p + f(xk);      //计算∑f(xk+h/2),T}                   //   k=0p = (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;}

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