|
《数值分析》第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，计算表格如下表所示。
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
#include
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 |