
function hw2_01

echo off, clc
clear all

  Nmin = 2;
  Nmax = 5;

% Exact value of the integral;
  Iex = - exp(-1.0) + 1.0;
%  Iex = 1 / (q + 1);
  
  i = 2;
  for n=Nmin:Nmax;
    Err(i,1) = n;
    Ig = Gauss_Integration( 0.0, 1.0, n );
    It = Trapezoidal_Integration( 0.0, 1.0, n );
    Err(i,2) = abs( (Ig - Iex) / Iex );
    Err(i,3) = abs( (It - Iex) / Iex );    
    i = i + 1;
  end;

% Plotting the result;
  figure(1); semilogy( Err(:,1), Err(:,2), 'x-', Err(:,1), Err(:,3), 'o-');
  legend('Gaussian', 'Trapezoidal'); axis( [ Nmin Nmax 1.0e-14 1.0] );
  




function [ I ] = Trapezoidal_Integration( a, b, N )

% Generating uniform grid;  
  x = linspace( a, b, N );
  h = ( b - a ) / ( N - 1 );
  
  
% Summing the nodes inside the domain;
  I = 0.0;
  for i=2:N-1
    I = I + fun( (i-1) * h + a );
  end;
  
% adding the endpoints (note different weights);
  I = I + 0.5 * ( fun(a) + fun(b) );
  
% rescaling;
  I = I * h;

%end;





function [ I ] = Gauss_Integration( a, b, N )
  
% rescaling the integration interval
  k1 = ( b - a ) / 2.0;
  k2 = ( b + a ) / 2.0;
  
% Determining the weights and the quadrature points
% depending on the order
  switch( N )
    case 1
     xi(1) = 0.000000000000000;
     w(1) = 2.000000000000000;

    case 2
     xi(1) = 0.577350269189626;
     xi(2) =-0.577350269189626;
     
     w(1) = 1.000000000000000;
     w(2) = 1.000000000000000;
     
    case 3
     xi(1) = 0.774596669241483;
     xi(2) =-0.774596669241483;
     xi(3) = 0.000000000000000;
     
     w(1) = 0.555555555555556;
     w(2) = 0.555555555555556;
     w(3) = 0.888888888888889;
  
    case 4
     xi(1) = 0.861136311594953;
     xi(2) =-0.861136311594053;
     xi(3) = 0.339981043584856;
     xi(4) =-0.339981043584856;
	     
     w(1) = 0.347854845137454;
     w(2) = 0.347854845137454;
     w(3) = 0.652145154862546;
     w(4) = 0.652145154862546;
     
    case 5
     xi(1) = 0.906179845938664;
     xi(2) =-0.906179845938664;
     xi(3) = 0.538469310105683;
     xi(4) =-0.538469310105683;
     xi(5) = 0.000000000000000;

     w(1) = 0.236926885056189;
     w(2) = 0.236926885056189;
     w(3) = 0.478628670499366;
     w(4) = 0.478628670499366;
     w(5) = 0.568888888888889;
     
   otherwise
    disp('Order of quadrature too high ...');

  end
  
% evaluating the integral
  I = 0.0;
  for i=1:N
    I = I + fun( k1 * xi(i) + k2 ) * w(i);
  end;
  
  I = I * ( b - a ) / 2.0;


% The function to be integrated
function [ y ] = fun( x )
%global q;
  
  y = exp( -x );
%  y = x.^q; 
