Points for attention in solving numerical integrals of Romberg (Romberg) quadrature formula

In the 5th edition of "Numerical Analysis" (Li Qingyang), the fourth chapter of Exercise 8th-(2) requires that the Romberg (Romberg) quadrature formula be used to solve the integral of f (x) =xsinx on the interval [0,2pi], and that the error is less than 10^ (-5).

To solve this problem, apply the calculation formula. T0=pi*[f (0) +f (2PI) was present in the first step when calculating the trapezoid formula. Obviously, if calculated according to sin (2pi) =sin (PI) =0, then the Romberg (Romberg) quadrature formula T0 the first three values of the column will appear 0, resulting in a slow calculation convergence, which needs to be iterated several times. (The actual calculation of the number of iterations is greater than 10 times, still can not meet the error less than 10^ (-5) conditions).

In order to solve this problem, the above calculation process should specify the effective digits of pi and be solved by rational number. In this topic because the question condition request error is less than 10^ (-5), therefore PI may take 6 decimal places or is greater than equals 6 decimal places. This article takes 7 decimal places, takes pi=3.1415927 as an example to calculate, at this time the computation only needs to iterate 5 times, can satisfy the question condition. The corresponding MATLAB program is given below

m_pi=3.1415927;   % 7 Decimal 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) Percent% 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) Percent% 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)

Run the above program, you will find that only 5 times to meet the error requirements, the result is 6.2831853, the calculation table as shown in the table below.

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

In addition, this paper gives the C + + program for Computing Romberg (Romberg) quadrature formula, as follows "annotation: The program data type needs to be further perfected--modified 2015-11-22".

#include <iostream> #include <math.h> using namespace std; # define Precision 0.00001//integral Precision Requirements # define E 2.71828183 #define MAXREPEAT 100//MAX Allow duplicate double function (double
	x)//integrand {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;//calculates 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);  Compute ∑f (XK+H/2), T}//K=0 P = (y[0] + h*p)/2.0;
		T ' m ' (H/2), variable step size trapezoid quadrature formula 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 ' (h)]
			/[pow (4,m) -1],2m-order Newton Christo formula, i.e. Romberg formula y[k-1] = p;
		p = q;
		EP = Fabs (Q-y[m-1]);//The result of the two-step calculation is more accurate m = m + 1; Y[m- 1] = q;   n = n + N;
	2 4 8 h = h/2.0;//twice times split interval num++;
	cout << "Iteration count" << num << "Times" << Endl;
return q;
	int main () {Double A, B, result;
	cout << "Please enter the lower limit of the points:" << Endl;
	Cin >> A;
	cout << "Please enter the integral limit:" << Endl;
	CIN >> B;
	result = Romberg (A, B, function);

	cout << "Romberg integration results:" << result << Endl;

	System ("pause");
return 0;
 }

The results of the program running are shown in the following figure.