function richardson
% richardson: Demonstrate Richardson's extrapolation.
%
%   richardson()
%
% Demonstrates the application of Richardson's extrapolation to improve
% the accuracy of the standard first order accurate forward finite
% difference formula for the first derivative.
%
% Essentially a script file (only a function to allow subfunctions).
%
% No Parameters.
%
% Ian Mitchell for CS 303, 3/11/07.

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Base stepsize (this should ideally be exactly represented in floating point).
% If the base stepsize is modified, the Richardson's extrapolation
% formula will probably be incorrect.
h_base = 0.5;
h_const = 1.0;

% stepsize powers (integers).
h_powers = -2:15;
num_powers = length(h_powers);

% Actual stepsizes.
h_values = h_const * h_base .^ h_powers;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% We will use 2 schemes: the standard first order accurate forward finite
% difference formula for the first derivative, and the derived
% Richardson's extrapolation formula.
num_schemes = 2;

% Some constants to make indexing more obvious.
i_standard = 1;
i_richardson = 2;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Function whose derivative we want.
diff_func = inline('atan(x)', 'x');

% Sample point.
x_0 = 1;

% True derivative.
soln_func = inline('1 ./ (1 + x.^2)', 'x');

% True solution.
soln_value = feval(soln_func, x_0);

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% For recording the values of each discretization and scheme.
approxs = zeros(num_powers, num_schemes);

% For recording the errors for each discretization and scheme.
errs = zeros(num_powers, num_schemes);

% Approximate the solution with the standard scheme.
for j = 1 : num_powers

  % What is the standard approximation?
  approxs(j, i_standard) = standard_approx(diff_func, x_0, h_values(j));
  errs(j, i_standard) = abs(soln_value - approxs(j, i_standard));
  
end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Approximate the solution with Richardson's extrapolation.  Since this
% formula requires the standard solution with h and h/2, we cannot create
% an approximation for the final stepsize.
for j = 1 : num_powers - 1

  approxs(j, i_richardson) = richardson_approx(approxs(j, i_standard), ...
                                               approxs(j+1, i_standard));
  errs(j, i_richardson) = abs(soln_value - approxs(j, i_richardson));
  
end

% We don't have a valid value for Richardson's formula for the final
% stepsize.  To flag this fact, set them to Not-a-Number.
approxs(end, i_richardson) = NaN;
errs(end, i_richardson) = NaN;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Plot empirical convergence curve, and some lines of fixed slope.
figure;
handle_standard = loglog(h_values, errs(:, i_standard), 'bo--');
hold on;
handle_richardson = loglog(h_values(1:end-1), ...
                           errs(1:end-1, i_richardson), 'b+-');

% Lines of various slopes for comparison.
loglog([ max(h_values), max(h_values) * 0.01 ], ...
       [ errs(1,1), errs(1,1) * 0.01 ], 'k-');
loglog([ max(h_values), max(h_values) * 0.01 ], ...
       [ errs(1,2), errs(1,2) * 0.01^2 ], 'k-');

% Make the plot pretty.
xlabel('Stepsize h');
ylabel('Error');
title('Empirical Convergence Rates');
legend('Forward Finite Difference', 'Richardson Extrapolation');
set(handle_standard, 'LineWidth', 2);
set(handle_richardson, 'LineWidth', 2);

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% Now plot the same curve on a linear scale.  Demonstrates that a linear
% scale does not display the convergence properties nearly so well as a
% loglog scale
figure;
handle_standard = plot(h_values, errs(:, i_standard), 'bo--');
hold on;
handle_richardson = plot(h_values(1:end-1), ...
                         errs(1:end-1, i_richardson), 'b+-');

% Make the plot pretty.
xlabel('Stepsize h');
ylabel('Error');
title('Empirical Convergence Rates (poorly presented)');
legend('Forward Finite Difference', 'Richardson Extrapolation');
set(handle_standard, 'LineWidth', 2);
set(handle_richardson, 'LineWidth', 2);

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% Now plot the approximations on a linear scale.  Using this plot, we can
% show Richardson's extrapolation graphically.
figure;
handle_standard = plot(h_values, approxs(:, i_standard), 'bo--');
hold on;
handle_richardson = plot(h_values(1:end-1), ...
                         approxs(1:end-1, i_richardson), 'b+-');

% Now let's show how Richardson's works for the third (and fourth) values
% of stepsize h.
fh = approxs(3, i_standard);
fh2 = approxs(4, i_standard);
c0 = 2 * fh2 - fh;
c1 = 2 / h_values(3) * (fh - fh2);
handle_interp = plot([ 0, h_values(2) ], [ c0, c0 + c1 * h_values(2) ], 'r:');

% Make the plot pretty.
xlabel('Stepsize h');
ylabel('Approximation');
title('Demonstrating Richardson Extrapolation Graphically');
legend('Forward Finite Difference', 'Richardson Extrapolation', ...
       'Interpolant v(h)');
set(handle_standard, 'LineWidth', 2);
set(handle_richardson, 'LineWidth', 2);
set(handle_interp, 'LineWidth', 2);



%---------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------------------------------------------------------------
function approx = standard_approx(f, x, h)
% standard_approx: Finite difference approx of the first dervative.
%
%   approx = standard_approx(f, x, h)
%
% A "standard" finite difference approximation (not using Richardson's
% extrapolation) for the first derivative.
%
% In this case, the first order accurate forward finite difference
% approximation.
%
% Input Parameters:
%
%   f: The function whose derivative is being approximated.
%
%   x: We will approximate f'(x).
%
%   h: The stepsize to use in the approximation.
%
% Output Parameters:
%
%   approx: The approximation.
  
  approx = (feval(f, x + h) - feval(f, x)) / h;

  
  
%---------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------------------------------------------------------------
function approx = richardson_approx(fh, fh2)
% richardson_approx: Approximation using Richardson's extrapolation.
%
%   approx = richardson_approx(fh, fh2)
% 
% A formula derived using Richardson's extrapolation to improve any first
% order accurate approximation to be at least second order accurate.
%
% Input Parameters:
%
%   fh: The approximation for some stepsize h.
%
%   fh2: The approximation for stepsize h / 2.
%
% Output Parameters:
%
%   approx: The approximation.
  
  approx = 2 * fh2 - fh;