#!/usr/bin/octave -q

% ----------------------------------------------------------------------
% J. R. Senning <jonathan.senning@gordon.edu>
% Gordon College
%
% This program uses five different root-finding algorithms to solve
% f(x) = 0.
% ----------------------------------------------------------------------

% ---------- Set up the parameters for the testing ---------------------

format long

% Display output if verbose = true
verbose = true;

% Maximum number of iterations
maxIter = 100;

% Tolerance for all algorithms.  This value determines how accurate
% the final estimate of the root will be.
tol = 1e-8;

% ---------- Define f(x), df(x) and the interval to search -------------

function y = f(x);  y = 1 - x * exp( x ); end
function y = df(x); y = -( 1 + x ) * exp( x ); end

a = 0.0;
b = 2.0;
first_guess = ( a + b ) / 2.0;

fprintf( '    ************************************************************\n' );
fprintf( '    Finding root of f(x) = 1 - x * exp(x) on [%f,%f]\n', a, b );
fprintf( '    ************************************************************\n' );

% ======================================================================
% ---------- Bisection Method ------------------------------------------
% ======================================================================

fprintf( '\n    ---- Bisection Method -------------------------------\n\n' );
[x, rate] = bisect( @f, a, b, tol, maxIter, verbose );
fprintf( '    root = ' );
disp( x );
fprintf( '    Estimated Convergence Rate = %5.2f\n', rate );

% ======================================================================
% ---------- Newton's method -------------------------------------------
% ======================================================================

fprintf( '\n    ---- Newton''s Method --------------------------------\n\n' );
[x, rate] = newton( @f, @df, first_guess, tol, maxIter, verbose );
fprintf( '    root = ' );
disp( x );
fprintf( '    Estimated Convergence Rate = %5.2f\n', rate );

% ======================================================================
% ---------- Secant method ---------------------------------------------
% ======================================================================

fprintf( '\n    ---- Secant Method ----------------------------------\n\n' );
[x, rate] = secant( @f, first_guess, tol, maxIter, verbose );
fprintf( '    root = ' );
disp( x );
fprintf( '    Estimated Convergence Rate = %5.2f\n', rate );

% ======================================================================
% ---------- Inverse Polynomial Interpolation method -------------------
% ======================================================================

fprintf( '\n    ---- Inverse Quartic Polynomial Interpolation -------\n\n' );
x0 = linspace( a, b, 4 );
[x, rate] = invpoly( @f, x0, tol, maxIter, verbose );
fprintf( '    root = ' );
disp( x );
fprintf( '    Estimated Convergence Rate = %5.2f\n', rate );

% ======================================================================
% ---------- Fixed point method ----------------------------------------
% ======================================================================

%
% The third parameter to this function should ideally be -1/df(r) where 
% r is the exact root and df(r) is the value of the first derivitive of
% f evaluated at r.  We can approximate df(f) with the centered
% difference formula
%               ( f( x + tol ) - f( x - tol ) )
%       df(r) = ------------------------------- + O(tol)
%                           2 * tol
%
% so that -1/df(r) = 2 * tol / ( f(x-tol) - f(x+tol) )
%
% For the "exact" root, we use the root returned by the last method...
% This is, of course, a little bogus since if we already have a good
% estimate of the root we don't need another one.  It illustrates,
% however, that a good choice of this parameter can really speed up
% convergence to the root.

dfr = 2 * tol / ( f( x - tol ) - f( x + tol ) );

fprintf( '\n    ---- Fixed-Point Method -----------------------------\n\n' );
[x, rate] = fixedpoint( @f, first_guess, dfr, tol, maxIter, verbose );
fprintf( '    root = ' );
disp( x );
fprintf( '    Estimated Convergence Rate = %5.2f\n', rate );

% ======================================================================
% ---------- Steffensen's method ---------------------------------------
% ======================================================================

%
% The value of dfr is computed above for the fixedpoint iteration.
%

fprintf( '\n    ---- Steffensen''s Method ----------------------------\n\n' );
[x, rate] = steffensen( @f, first_guess, dfr, tol, maxIter, verbose );
fprintf( '    root = ' );
disp( x );
fprintf( '    Estimated Convergence Rate = %5.2f\n', rate );

% End of roots_demo