% BRIEF RECAP of October 27th, 1998 NEEP271 lecture
% written by Fabian Waleffe 98/10/27
%
% BEGIN QUADRATURE = INTEGRATION
%
% example: probability to be within 2 std from mean for gaussian pdf 
%
%  prob(-2s,2s)= integral from -2s to 2s of p(x) dx
%
%  p(x) =1/sqrt(2 pi s^2) exp(-x^2/(2 s^2))   (mean 0, std s)
%
%  change variables: u=x/(sqrt(2) s) then
%
%  prob(-2s,2s) = 2/sqrt(pi) integral from 0 to sqrt(2) of exp(-u^2) du
%               = erf(sqrt(2)) = "error function"
%
% Matlab knows the error function: erf(x), see  "help erf"
%
% NUMERICAL INTEGRATION: Riemann sum, Trapezoidal rule, Midpoint rule
%
% Integral from x_1 to x_N :
%
% RS = SUM      f_i         Delta_i
%
% TR = SUM [f_i+f_(i+1)]/2  Delta_i 
%
% MR = SUM   f_(i+1/2)      Delta_i
%
%      where x_1, x_2, ..., x_N is a partition of the interval [x1 x_N]
%      Delta_i = x_(i+1)-x_i,  the SUM is from i=1 to N-1
%      f_i=f(x_i) and x_(i+1/2) = [x_i+x_(i+1)]/2

% EXAMPLE:
% 2/sqrt(pi) * integral from 0 to sqrt(2) of exp(-x^2) dx = erf(sqrt(2))

  Exact=erf(sqrt(2))

N=10;

  x=[0:1/N:1]'*sqrt(2); dx=sqrt(2)/N;  y=exp(-x.^2); 

  TR1= dx* (sum(y) -y(1)/2 -y(length(x))/2) 

% Let's repeat this for more segments:

% 20 segments
  N=20; x=[0:1/N:1]'*sqrt(2); dx=sqrt(2)/N;  y=exp(-x.^2); 
  TR2= dx* (sum(y) -y(1)/2 -y(length(x))/2) ; 

% 40 segments
  N=40; x=[0:1/N:1]'*sqrt(2); dx=sqrt(2)/N;  y=exp(-x.^2); 
  TR3= dx* (sum(y) -y(1)/2 -y(length(x))/2) ; 

% 80 segments
  N=80; x=[0:1/N:1]'*sqrt(2); dx=sqrt(2)/N;  y=exp(-x.^2); 
  TR4= dx* (sum(y) -y(1)/2 -y(length(x))/2) ; 

% put all these results in a vector "TR"
  TR=[TR1; TR2; TR3; TR4]

% compute relative error

  error=abs(2/sqrt(pi)*TR-Exact)/abs(Exact)

  loglog(sqrt(2)*[1/10;1/20;1/40;1/80]',error,'-o')
  grid, xlabel('dx'), ylabel('Relative error')
  text(0.1,1.4e-4,'What is the relation between error and dx?')

% That's simple enough but in fact TR is built-in:

  [trapz(x,y), TR4]

  help trapz

  help quad