Example 3.1

Numerical integration of int(sin(x)/(2*x^3),0,pi)

Contents

Initialization

clear N h X f I A B C err_A err_B df_a df_b C;
clear err_C;
close all;
format long;

N = 8;

h = (pi-1)/N;  % uniform spacing

X = [1:h:pi];
f = sin(X)./(2*X.^3);

Resolution

% A : Trapezoidal method
% B : Simpson
% C : End-corrected Trapezoidal
%
% err_A, err_B, err_C : % of error with 'exact' sol

I = 0.1985572988;
A = 0;         % trapezoidal
B = 0;         % quadrature
C = 0;         % correct end

%  Trapezoidal Rule
%--------------------
for i = [2:N]
    A = A + f(i);
end;
A = A + 1/2 * (f(N+1) + f(1));
A = A*h;
err_A = abs(I-A)/I * 100;

%  Simpson (quad)
%--------------------
for i = [1:N/2-1]
    B = B + 4*f(2*i) + 2*f(2*i+1);
end;
B = B + f(1) + f(N+1) + 4*f(N);
B = B*h/3;
err_B = abs(I-B)/I * 100;


%  Trap. End Correct
%--------------------

df_a = (cos(1)-3*sin(1))/2;
df_b = -1/(2*pi^3);

C = A - (h^2)/12 * (df_b - df_a);
err_C = abs(I-C)/I * 100;

Display results

disp (['N = ' int2str(N) '         Result' '    Error']);
disp(['Trapezoidal   ' num2str(A) '    ' num2str(err_A)]);
disp(['Simpson       ' num2str(B) '   ' num2str(err_B)]);
disp(['Trap end corr ' num2str(C) '   ' num2str(err_C)]);
disp('   ');

N = 8         Result    Error
Trapezoidal   0.2043    2.8943
Simpson       0.19883   0.1396
Trap end corr 0.19848   0.040948
   

Published with MATLAB® 7.5