%{
Integration: Trapezoid, Simpson and Gauss
Effect of number of intervals on relative error
AG: February 2022
%}
clc
tic
clf 
fun = @f;

% limits of integration and initial number of intervals
a =  0.0;
b =  2.0;
ntimes = 18;    % how many times to double intervals
eps = 1.0e-12;

% reserve arrays for plots
xerror = zeros(ntimes,1);
yerr1  = zeros(ntimes,1);
yerr2  = zeros(ntimes,1);
yerr3  = zeros(ntimes,1);
yerr4  = zeros(ntimes,1);

% Plot the function
figure (1)
%figure('DefaultAxesFontSize',14)
title('f(x)')
xlabel('x')
ylabel('f(x)')
xf = a: (b-a)/1000.: b;
yf = f(xf);
plot(xf,yf,'r','LineWidth',1.5)
grid;

% matlab integrals
I1 = integral(fun,a,b,'RelTol',0,'AbsTol',eps);
I2 =   quadgk(fun,a,b,'RelTol',0,'AbsTol',eps,'MaxIntervalCount',1024);
fprintf('quadgk = %12.8f \n',I2)

 if abs(I1-I2) >= eps
   fprintf(' integral and quadgk do not match')
 end

% calculate integrals using four methods
fprintf('          n    Trapezoid    Simpson      Gauss8       Gauss16 \n')
nint = 1;
for j = 1:ntimes

  nint = nint*2;
  [ITrap] = trapezoid(fun,a,b,nint);
  [ISimp] = simpson(fun,a,b,nint);
  [G8]    = gauss8(fun,a,b,nint);
  [G16]   = gauss16(fun,a,b,nint);
  fprintf(' %10i  %12.8f %12.8f %12.8f %12.8f\n',nint,ITrap,ISimp,G8,G16)

% for calculating relative errors
  xerror(j) = (nint);
  yerr1(j)  = ((abs(I1-ITrap))/I1);
  yerr2(j)  = ((abs(I1-ISimp))/I1);
  yerr3(j)  = ((abs(I1-G8))/I1);
  yerr4(j)  = ((abs(I1-G16))/I1);

end

% Plot errors
figure (2)
%figure('DefaultAxesFontSize',14)
loglog(xerror,yerr1,xerror,yerr2,xerror,yerr3,xerror,yerr4,'LineWidth',1.5)
legend('Trapezoid','Simpson','Gauss8','Gauss16')
% loglog(xerror,yerr1,xerror,yerr2,'LineWidth',1.5)
% legend('Trapezoid','Simpson')
% title('Relative error')
xlabel('n')
ylabel('relative error')

%print time of calculations
timeE = ceil(toc);
tminutes = int32(floor(timeE/60));
tseconds = timeE - tminutes*60;
fprintf('\n Elapsed time:  %2i minutes %2i seconds \n',tminutes,tseconds)


% Integrand ===========================
function f = f(x)
  %f = sin(x);
  f = exp(-x);
  %f = (x.*cos(10*x.^2))./(1.0+x.^2);
  %f = x.*sin(30*x).*cos(50*x);
  %f = 1.0 ./ (1.0- 0.998 * x.^2);
end
% ======================================

function [tr] = trapezoid(fun,a,b,nint)
  r = 0.0;
  dx = (b-a)/nint;
  for j=1:nint-1
    x = a + j*dx;
    r = r + fun(x);
  end
  r = (r + (fun(a)+fun(b))/2.0)*dx;
  tr = r;
end

function [sm] = simpson(fun,a,b,nint)
  if 2*floor(nint/2) ~= nint
    nint=nint+1;
  end
  s = 0.0;
  dx = (b-a)/nint;
  for i=2:2:nint-1
    x = a+i*dx;
    s = s + 2.0*fun(x) + 4.0*fun(x+dx);
  end
  s = (s + fun(a)+fun(b)+4.0*fun(a+dx) )*dx/3.0;
  sm = s;
end

function [G8] = gauss8(fun,a,b,nint)
n =4;
ti = [0.1834346424 0.5255324099 0.7966664774 0.9602898564];
ci = [0.3626837833 0.3137066458 0.2223810344 0.1012285362];
G8 = 0.0;
dx = (b-a)/nint;
for j=1:nint
  r = 0.0;
  ai = a + (j-1)*dx;
  bi = a + j*dx;
  m = (bi-ai)/2.0;
  c = (bi+ai)/2.0;
  for i=1:n
    r = r + ci(i)*(fun(m*(-1.0)*ti(i) + c) + f(m*ti(i) + c));
  end
  G8 = G8 + r*m;
end
end

function [G16] = gauss16(fun,a,b,nint)
n =8;
ti = [0.0950125098 0.2816035507 0.4580167776 0.6178762444 ...
      0.7554044083 0.8656312023 0.9445750230 0.9894009349]; 
ci = [0.1894506104 0.1826034150 0.1691565193 0.1495959888 ...
      0.1246289712 0.0951585116 0.0622535239 0.0271524594]; 
G16 = 0.0;
dx = (b-a)/nint;
for j=1:nint
  r = 0.0;
  ai = a + (j-1)*dx;
  bi = a + j*dx;
  m = (bi-ai)/2.0;
  c = (bi+ai)/2.0;
  for i=1:n
    r = r + ci(i)*(fun(m*(-1.0)*ti(i) + c) + f(m*ti(i) + c));
  end
  G16 = G16 + r*m;
end
end


% for the oscillating integral
% x = 2*pi;
% Iu = sin(80*x)/12800 - sin(20*x)/800 + (x*cos(20*x))/40 - (x*cos(80*x))/160;
% x = 0;
% Id = sin(80*x)/12800 - sin(20*x)/800 + (x*cos(20*x))/40 - (x*cos(80*x))/160;
% Ia = Iu-Id;
%fprintf('ratio %18.16f \n',I1/Ia)
%I1 = Ia;