% Math 32a   Fall 2003  Richard Palais  Brandeis Univ.

%       Discussion of Project 2

% Let's first write the code for the trapezoidal approximation 
% to the integral of a function f over [a,b] without subdivision.

 trapezoid = 0.5*(b-a)*(feval(f,a) + feval(f,b));
 
%  We recall that a bound for the error is given by 

 ErrT = C1*(b-a)^3/12;
 
%  where C1 is the maximum of the absolute value of f'' on [a,b] .

%  Similarly, the Simpson's Rule approximation to the integral 
%  of f over the same interval, without subdividing, is  

Simpson = ((b-a)/6)*(feval(f,a) + 4 * feval(f,(a+b)/2) + feval(f,b));

%  and in this case, a bound for the error is given by 

 ErrS = C2*(b-a)^5/90;

%  where C1 is the maximum of the absolute value of f'''' on [a,b] .
 
% Lets try this out on a simple function, say sin:
 
 f = 'sin';
 
%  on the interval [0,pi/2]
 
 a = 0;
 b = pi/2;
 
 % Since -cos is an antiderivative for sin, the exact value 
 %  of the integral is
 
  integral = (-cos(pi/2) - ( - cos(0)))
 
% Also, clearly C1 and C2 = 1;

C1 = 1;
C2 = 1;

%  Now lets design the M-Files.  

%     First       TrapezoidalRule1.m
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
 function    ApproxInt = TrapezoidalRule1(f,a,b,n)
 
%  This computes the approximation to the integral 
%  of the function f from a to b by dividing the 
%  interval into n equal parts and applying the 
%  Trapezoidal Rule approximation to each part and
%  then summing.
 
 %  Richard Palais
 
 h = (b-a)/n;
 ApproxInt = 0;
 for k = [0:n-1] 
   a1 = a + k*h;
   b1 = a1 + h;
   ApproxInt = ApproxInt + trapezoid(f,a1,b1);
 end;
 
 
        %%%%%%% Auxiliary Function  %%%%%%%
		
 function      term = trapezoid(f,a,b)
 
     term  = 0.5*(b-a)*(feval(f,a) + feval(f,b));
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	 
%	 Next    SimpsonsRule1.m
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
  function    ApproxInt = SimpsonsRule1(f,a,b,n)
 
%  This computes the approximation to the integral 
%  of the function f from a to b by dividing the 
%  interval into n equal parts and applying the 
%  Simpson's Rule approximation to each part and
%  then summing.
 
 %  Richard Palais
 
 h = (b-a)/n;
 ApproxInt = 0;
 for k = [0:n-1] 
   a1 = a + k*h;
   b1 = a1 + h;
   term = Simpson(f,a1,b1);
   ApproxInt = ApproxInt + term;
 end;
 
      %%%%%%% Auxiliary Function  %%%%%%%
	  
 function      term = Simpson(f,a,b)
 
   term = ((b-a)/6)*(feval(f,a) + 4 * feval(f,(a+b)/2) + feval(f,b));
 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%	The second versions below avoid he auxiliary 
%	function calls, AND more importantly, they 
%	avoid re-evaluating the integrand twice for  
%	each interior endpoint of an interval
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
%	TrapezoidalRule2.m
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	function    ApproxInt = TrapezoidalRule2(f,a,b,n)
 
 %  This computes the approximation to the integral 
 %  of the function f from a to b by dividing the 
 %  interval into n equal parts and applying the 
 %  Trapezoidal Rule approximation to each part and
 %  then summing.
 
 %  Richard Palais
 
 h = (b-a)/n;
 sum = (feval(f,a) + feval(f,b))/2 ;
 x = a;
 for k = [1:n-1] 
     x = x + h;
     sum = sum + feval(f,x);
 end;
 ApproxInt = h * sum;
 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
%    SimpsonsRule2.m
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
	
	function    ApproxInt = SimpsonsRule2(f,a,b,n)

%  This computes the approximation to the integral 
%  of the function f from a to b by dividing the 
%  interval into n equal parts and applying the 
%  Simpson's Rule approximation to each part and
%  then summing. It gains efficiency by grouping 
%  terms so that the function is not evaluated
%  more than once at the same argument.

 %  Richard Palais
 
h = (b-a)/n;
x = a;
z = a + h/2;
sum = (feval(f,a) + feval(f,b)) + 4 * feval(f,z);
   for k = [1:n-1] 
       x = x + h;
       z = z + h;
       sum = sum + 2 * feval(f,x) + 4 * feval(f,z);
   end;
ApproxInt = (h/6) * sum;   

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%	The third versions below avoid the loop by  
%	vectorizing the sum, and they also avoid 
%	evaluating the integrand more than once for
%	a given value of the argument.  The code is by 
%   Nick Dufresne. (slightly modified).
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%	TrapezoidalRule3.m
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%nick dufresne 
%
%usage:
%      result = TrapezoidalRule3(f,a,b,n)
%
% f is an inline function that maps a real input to a vector 
% a is the starting point
% b is the ending point
% n is the number of sub-intervals into which [a,b] is divided
% result is the estimated value of the integral
%

function thesum = TrapezoidalRule (f,a,b,n)
delta = (b-a)/n;
arguments = [a:delta:b];

%calculate all the values of the function
%in summing all n parts of the trapezoidal method we count internal points
%twice so we will multiply by 2

values = 2*feval(f,arguments);

%since we have calculated the value of the endpoints twice we need to
%subtract the values and we need to multiply by 0.5*delta  = 1/2*(b-a)/n

thesum = 0.5*delta*(sum(values)-feval(f,a)-feval(f,b));
 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
%    SimpsonsRule3.m
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
%nick dufresne 
%
%usage:
%      result = SimpsonsRule3(f,a,b,n)
%
% f is an inline function that maps a real input to a vector 
% a is the starting point
% b is the ending point
% n is the number of sub-intervals into which [a,b] is divided
% result is the estimated value of the integral
%

function thesum = SimpsonsRule (f,a,b,n)
delta = (b-a)/n;
arguments = [a:delta:b];
midpointArgs = [a+delta/2:delta:b-delta/2];
%calculate all the values of the function.
%in summing all n parts of the trapezoidal method, 
%we need to count internal pointstwice so we multiply by 2.

values = 2*feval(f,arguments);

%in the formula all midpoint values are multiplied by 4

midpointValues = 4*feval(f,midpointArgs);

%since we have calculated the value of the endpoints twice we need to
%subtract the values and we need to multiply by (1/6)*delta  = 1/6*(b-a)/n

thesum = (1/6)*delta*(sum(midpointValues)+sum(values)-feval(f,a)-feval(f,b));
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 format long
 
  f = inline('4./(1+ x.^2)','x');  % must be vectorized for third versions!
 
  a = 0;
  
  b = 1;
  
  C1 = 12;
 
  C2 = 96;
 
%%%%%%%%

  n=10000;
  
tic, A = TrapezoidalRule1(f,a,b,n); toc  %elapsed_time = 9.210232

tic, A = TrapezoidalRule2(f,0,1,n); toc  %elapsed_time = 4.533945

tic, A = TrapezoidalRule3(f,0,1,n); toc  %elapsed_time = 0.009738  !!!!! 

 A

 ErrT = C1*(b-a)^2/(12 * n^2)

 Error  = abs(A-pi)


tic, A = SimpsonsRule1(f,a,b,n); toc

tic, A = SimpsonsRule2(f,a,b,n); toc

tic, A = SimpsonsRule3(f,a,b,n); toc

  A

 ErrT = C1*(b-a)^2/(12 * n^2)

 Error  = abs(A-pi)

%%%%%

n = 3;

tic, A = SimpsonsRule1(f,a,b,n); toc   %elapsed_time = 0.005829

tic, A = SimpsonsRule2(f,a,b,n); toc   %elapsed_time = 0.004387

tic, A = SimpsonsRule3(f,a,b,n); toc   %elapsed_time = 0.002445n


 A

 ErrS = C2*(b-a)^4/(90 * n^4)

 Error = abs(A-pi)
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%   Let's check how the errors vary as a function 
%%%   of the  number n of subdivisions.
%%%   First Trapezoidal:
%
n = 5;
A = TrapezoidalRule1(f,a,b,n);
Error5 = abs(A-pi)          %   0.00666653977880
n=10;
A = TrapezoidalRule1(f,a,b,n);
Error10 = abs(A-pi)         %   0.00166666468263
%
%  Since Trapezoidal is a quadratic method, the ratio
%  of Error5 to Error 10 should be about (10/5)^2 = 4
%
 ratio = Error5/Error10    %  3.99992862887449
%
%%%%   Now Simpson's 
%
n = 5;
A = SimpsonsRule1(f,a,b,n);
Error5 = abs(A-pi)         % 3.965057793209326e-08
n=10;
A = SimpsonsRule1(f,a,b,n);
Error10 = abs(A-pi)        %  6.200080449048073e-10
%
%  Since SimpsonsRule1 is a fourth order method, the ratio
%  of Error5 to Error 10 should be about (10/5)^4 = 16
%  BUT,IN FACT:
 ratio = Error5/Error10     %  63.95171523650308  ~ (10/5)^6
%
n = 20; 
A = SimpsonsRule1(f,a,b,n);
Error20 = abs(A-pi)         %   9.687362023669266e-12
%
 ratio = Error10/Error20    %  64.00174200055011 ~  (20/10)^6
 
%  So, the conclusion seems nearly inescapable that for 
%  the present function, Simpson's is behaving as a 
%  sixth order method !!!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Shootout:   Comparison of Nick Dufresne's vectorized Trapezoidal and 
%              Simpson's with non-vectorized versions.

  f = inline('4./(1+ x.^2)','x');   % must be vectorized for third versions!
 
  a = 0;
  
  b = 1;

n=100000;

tic, A = TrapezoidalRule3(f,a,b,n); toc    %   elapsed_time = 0.07617

Error  = abs(A-pi)                         %  1.664046678229170e-11

n=14;

tic, A = SimpsonsRule3(f,a,b,n); toc       %  elapsed_time = 0.002523

Error  = abs(A-pi)                         %  8.234657400407741e-11

n = 100000;

tic, A = TrapezoidalRule2(f,a,b,n); toc    %  elapsed_time = 49.330224

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
