%% ME 2016 Section
%% Fall 2007
%%
%% AUTHOR:  Chris Paredis
%%
%% PURPOSE: compare the accuracy and order of convergence for
%%          different numerical integration methods
%%
%% This script uses trapezoidal and Simpson integration to compute a known
%% integral for different numbers of integration segments.  The results are
%% plotted in a log-log plot so that they can be compared to theoretical
%% solutions
%%
%% ASSUMPTIONS:
%%         - none

%%
%% Problem formulation:
%%

% to avoid confusion with previously defined variables, it is sometimes
% good to clear the workspace before starting the real computations
% (be careful: this assumes of course that there is nothing in the
% workspace that you still need -- use this with caution)
clear all;
close all;

% define the integration problem to be solved
f = @(x)(cos(x).*cos(2*x));
lower_bound = 0;
upper_bound = 2;
exact_integ = sin(2)/2+sin(6)/6;

% consider the different number of segments -- making sure the number of
% segments is a whole number
% NOTE: logspace uses as its first two inputs the powers of 10 so that the
% range here extends from 10^1 to 10^6
segments = floor(logspace(1,6,20));

% assign memory for storing the errors
trapez_error = zeros(1,length(segments));
simpson_error = zeros(1,length(segments));

% do the same for the quad() errors
% note: we need to create a vector of desired accuracies
accuracies = logspace(-1,-16,20);
quad_error = zeros(1,length(accuracies));
quad_segments = zeros(1,length(accuracies));

%%
%% Problem solution:
%%

% compute the errors for trapezoidal and simpson
for i=1:length(segments)
    % perform the actual integration using two intetegration methods
    trapez_error(i) = abs(trapezoidal_integration(f, lower_bound, upper_bound, segments(i))-exact_integ);
    simpson_error(i) = abs(simpson_integration(f, lower_bound, upper_bound, segments(i))-exact_integ);
end

% compute the errors for quad()
for i=1:length(accuracies)
    [integral, f_evals] = quad(f, lower_bound, upper_bound, accuracies(i));
    quad_error(i) = abs(integral - exact_integ);
    quad_segments(i) = f_evals-1;
end

%%
%% Problem interpretation:
%%

% plot the error as a function of number of segments in a log-log plot.
loglog(segments,trapez_error,'x',segments,simpson_error,'o');
hold on
loglog(quad_segments, quad_error, 'r+');
legend('trapezoidal','simpson 1/3','quad()');
grid on
xlabel('number of segments')
ylabel('absolute error')

%%
% compute the order of convergence of each of the methods
% 1) trapezoidal:
coeffs = polyfit(log10(segments),log10(trapez_error),1);
fprintf('Trapezoidal order of convergence = %g\n',-coeffs(1));

% 2) simpson: because round-off is starting to kick in, we will only
% consider the first 12 points
range = 1:12;
coeffs = polyfit(log10(segments(range)),log10(simpson_error(range)),1);
fprintf('Simpson order of convergence = %g\n',-coeffs(1));

% 3) quad: because round-off is starting to kick in, we will only
% consider the first 16 points -- we cut off the first two also because
% they are identical to the third
range = 3:16;
coeffs = polyfit(log10(quad_segments(range)),log10(quad_error(range)),1);
fprintf('Quad order of convergence = %g\n',-coeffs(1));

%% ANSWERS TO THE INTERPRETATION QUESTIONS

%% Question 1: Describe the most important characteristics of the graph.
% The relation between the number of segments and the absolute error is
% more or less linear in a log-log scale.  This true for each of the three
% methods although the quad() performance has a fair amount of noise in it.
% Among the three methods, the error reduces most quickly for quad(),
% then Simpson, and finally trapezoidal.
% For Simpson and quad, the performance levels off at around 10^-15 and
% 10^-17, respectively.  After these points, the error slowly increases
% with the number of segments. 

%% Question 2: What is the order of convergence of all three integration
%% methods?
% The order of convergence of each of the methods can be computed by using
% a first order polynomial curve fit as applied to the log-log values (as
% is shown in the script above).  The order of convergence is the negative
% of the slope of this straightline fit.  The results of the script are:
% Trapezoidal order of convergence = 1.99959
% Simpson order of convergence = 3.93674
% Quad order of convergence = 6.1404

%% Question 3: Does the order of convergence match the theoretical results
%% that we derived in class?
% Yes, the results for trapezoidal and Simpson match pretty closely.  The
% theoretically predicted order of convergence for each is 2 and 4
% respectively.  The difference is due to the fact that the theoretical
% result really only applies in the limit when the number of segments
% approaches infinity.  The theoretical result also does not consider the
% round-off due to the use of floating point numbers.

%% Question 4: Why does the error for the Simpson rule seem to increase for
%% large numbers of segments?
% Since the algorithms are implemented using floating point computations,
% small round-off errors occur whenever we add or multiply two numbers.
% Although these errors are very small, when the algorithmic error becomes
% small, the round-of errors start to dominate. As a result, the error can
% never become smaller than the machine constant (about 10^-15 relative
% error -- Note: relative, not absolute). In addition, when a large number
% of arithmetic operations are performed, the total round-off error grows.
% As you can notice for Simpson, the growth rate is approximately linear in
% the number of arithmetic operations (slope of +1).

%% Question 5: Which method would you use and how many segments would you
%% select? Justify your answers.
% There are many correct answers for this question.  From the figure, it is
% clear though that using the Matlab quad() function typically gives the
% best performance.  For instance, to obtain an error of 10^-10, quad only
% needs about 130 function evaluations while Simpson requires about 380 and
% trapezoidal needs about 64,000!!  However, how many segments one should
% use depends on how accurately one wants to know the result.  If one wants
% to know the result as accurately as possible, one should have used quad()
% with approximately 3000 segments.  If one is satisfied with an error of
% 10^-10 then using quad with 130 function evaluations would be best.
% Notice that for very few function evaluations (or very inaccurate
% results), Simpson may actually work slightly better than quad().
% However, it is rather unlikely that anybody would be interested in
% results with such a poor accuracy.