Numerical integration (Trapezoid, Simpson, Gauss)

Hello again. I have to evaluate the integral of the function exp(x^2) from 0 to 2 using the following methods: composite trapezoid, composite Simpson's, Gauss' 2 and 3 points. The best result in my code is obtained by composite Simpson's rule. My question is: shouldn't Gauss' 3 points method give the best value for the integral ? My code:

clc, clear all, close all, format long e
% Numerical integration - Finding the mass of a rod of length l = 2 m and
% mass density of exp(x^2) kg/m where x is the distance frome the center of
% the rod

n = 10; % Partition
a = -1; % Lower bound
b = 1; % Upper bound
h = (b-a)/n;
x = [a:h:b]; % Domain
f = @(x) exp(x.^2); % function to integrate
N = length(x);
Q = integral(f,-1,1);

%--- Composite trapezoid rule
F = zeros(1,N);

for i = 1:N
    F(1,i) = f(x(i));
end
F = sum(F);
I_trap = (h/2)*(2*F - f(a) - f(b));

fprintf('\n\tComposite trapezoid gives:\n\n\t\t\t\t\t\tm = %f kg\n',I_trap)

	Composite trapezoid gives:

						m = 2.961309 kg

fprintf('\n\twith a deviation D = %f%%\n',abs(Q-I_trap)*100/Q)

	with a deviation D = 1.230835%

%--- Composite Simpson's rule

F1 = zeros(1,n/2-1);

for i = 1:n/2
    F1(1,i) = f(x(2*i));
end

F2 = zeros(1,n/2-1);

for j = 1:n/2-1
    F2(1,j) = f(x(2*j+1));
end

I_Simpson = (h/3)*(4*sum(F1) + 2*sum(F2) + f(a) + f(b));

fprintf('\n\tComposite Simpson''s  gives:\n\n\t\t\t\t\t\tm = %f kg\n',I_Simpson)

	Composite Simpson's  gives:

						m = 2.926204 kg

fprintf('\n\twith a deviation D = %f%%\n',abs(Q-I_Simpson)*100/Q)

	with a deviation D = 0.030780%

%---- Gauss

L = [0,sqrt(3)/3,sqrt(15)/5,sqrt(((15 + 2*sqrt(30)))/35)]; % Roots of Legendre polynomials
c = [2,1,8/9,5/9,(90 - 5*sqrt(30))/180,(90 + 5*sqrt(30))/180]; % Coefficients

%--- Two points Gauss

I_Gauss2 = f(L(2)) + f(-L(2));

fprintf('\n\tGauss'' 2 points gives:\n\n\t\t\t\t\t\tm = %f kg\n',I_Gauss2)

	Gauss' 2 points gives:

						m = 2.791225 kg

fprintf('\n\twith a deviation D = %f%%\n',abs(Q-I_Gauss2)*100/Q)

	with a deviation D = 4.583410%

%--- Three points Gauss

I_Gauss3 = (5/9)*f(L(3)) + (8/9)*f(0) + (5/9)*f(-L(3));

fprintf('\n\tGauss''s 3 points gives:\n\n\t\t\t\t\t\tm = %f kg\n',I_Gauss3)

	Gauss's 3 points gives:

						m = 2.913465 kg

fprintf('\n\twith a deviation D = %f%%\n',abs(Q-I_Gauss3)*100/Q)

	with a deviation D = 0.404681%

Answer 1

[John D'Errico]

"shouldn't Gauss' 3 points method give the best value for the integral ?"

Why? In fact, no. There is no reason for that to be true.

An n-node Gaussian rule will be order 2*n-1 as I recall. So a 2 node rule will be exact for a cubic polynomial on some interval. A 3 node rule will be exact for a degree 5 polynomial. All good. But, is your function a low order polynomial? OF COURSE NOT!

As such, none of these rules will be exact for the function exp(x^2). They will be approximations. Merely that. And yes, they may be decent approximations.

Look at the Simpson's rule. If it was only ONE panel of the Simpson's rule, then it would be comparable in some sense to the 3 point Gaussian rule, and then you would see it is less accurate. But by using multiple panels of the Simpson's rule (thus a composite Simpson's rule) you get a lower error. Somewhere in class you should have learned how the error decreases as a function of the step size.

So this is not at all unexpected.

Now, suppose you created a composite Gaussian 3 node rule? So by using multiple panels of a Gaussian rule? Here each panel would live on a much shorter interval. And now the composite Gaussian rule would indeed be more accurate than the Simpson's rule.

(Aside) Is that really a fully valid comparison? Sigh. Probably not. A composite Simpson's rule comes from the family of Newton-Cotes rules. They have the virtue that each panel shares some function values with the panels on either side. And that makes the rule a little more efficient, since what often matters is the total number of function evaluations.