% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%