% Steffensen's method for 1D root finding.
% K. Ming Leung, 04/24/05

clear all; format long g;
fcn = inline('(x-1)^2'); 
fcn = inline('cos(x)'); %
fcn = inline('x*x-1');
%fcn = inline('x*x-4*sin(x)');   % define the function using inline.
x = 3;                          % initial guess for the root.
tol = 100000000*eps;                    % set the tolerance

dx = 1;                         % the change in x (the step)
loop = 1;                       % loop counter
itMax = 50;                     % set the maximum number of iterations
fprintf('\n  loop          x                      f(x)                     dx\n');

while abs(dx) > tol && loop < itMax
    f = feval(fcn,x);           % function value at x
    ff = feval(fcn,x+f);        % function value at x+f
    dx = (f*f)/(f-ff);          % the change in x
    x = x + dx;                 % update x 
    fprintf('%4d %20.14f %25.14e %25.14e\n',loop,x,f,dx);
    %disp([loop x f dx]);        % display intermediate results
    loop = loop+1;              % increment counter
end

% Results for the function x*x-4*sin(x):
% loop          x                      f(x)                     dx
%    1     2.43505086117835     8.43551996776053e+000    -5.64949138821653e-001
%    2     2.08694437317843     3.33264347351268e+000    -3.48106487999916e-001
%    3     1.98020118336022     8.76430032570448e-001    -1.06743189818208e-001
%    4     1.93990029096082     2.51765166153126e-001    -4.03008923994033e-002
%    5     1.93387949970100     3.26091721726476e-002    -6.02079125981839e-003
%    6     1.93375381674873     6.64900427585824e-004    -1.25682952271405e-004
%    7     1.93375376282703     2.85120195453459e-007    -5.39216991425897e-008
%    8     1.93375376282702     5.24025267623074e-014    -9.90945217620556e-015
%    9     1.93375376282702    -4.44089209850063e-016     8.88178419700125e-017
%
% dx is less than eps and so the result for the root is expected to be good to eps.
% From the function values for loop 5-8, it appears that the convergence is quadratic.
% Similar result if we consider dx.
%
% Testing the method using other functions:
% Convergence is quadratic at simple roots.
% Sometimes it is difficult to have very highly accurate roots due to 
% the near cancellation of f(x+f) and f(x) in computing dx (line 19).
% Just like Newton's method, convergence depends whether the initial
% guess is near the root.
% 
% Results for the function (x-1)^2, double root at x=1:
%   loop          x                      f(x)                     dx
%    1     2.50000000000000     4.00000000000000e+000    -5.00000000000000e-001
%    2     2.07142857142857     2.25000000000000e+000    -4.28571428571429e-001
%    3     1.72259136212625     1.14795918367347e+000    -3.48837209302325e-001
%    4     1.45718562692185     5.22138276619464e-001    -2.65405735204393e-001
%    5     1.27112494761755     2.09018697463928e-001    -1.86060679304304e-001
%    6     1.15174580561923     7.35087372206185e-002    -1.19379141998314e-001
%    7     1.08122362533988     2.30267895230306e-002    -7.05221802793594e-002
%    8     1.04219676424194     6.59727731335243e-003    -3.90268610979325e-002
%    9     1.02153432613569     1.78056691249010e-003    -2.06624381062537e-002
%   10     1.01088185990878     4.63727202118221e-004    -1.06524662269057e-002
%   11     1.00547037347304     1.18414875074388e-004    -5.41148643574709e-003
%   12     1.00274264757622     2.99249859344999e-005    -2.72772589681785e-003
%   13     1.00137320174177     7.52211572733693e-006    -1.36944583444666e-003
%   14     1.00068707196815     1.88568302360511e-006    -6.86129773621954e-004
%   15     1.00034365396054     4.72067889417193e-007    -3.43418007611140e-004
%   16     1.00017185649983     1.18098044593860e-007    -1.71797460709576e-004
%   17     1.00008593563284     2.95346565334688e-008    -8.59208669872514e-005
%   18     1.00004296966197     7.38493299191385e-009    -4.29659708748701e-005
%   19     1.00002148529349     1.84639184954265e-009    -2.14843684796301e-005
%   20     1.00001074276258     4.61617836232076e-010    -1.07425309024190e-005
%   21     1.00000537141337     1.15406947954166e-010    -5.37134921438877e-006
%   22     1.00000268570927     2.88520815971757e-011    -2.68570410071623e-006
%   23     1.00000134287221     7.21303428225131e-012    -1.34283706159182e-006
%   24     1.00000067140620     1.80330576788053e-012    -6.71466011757772e-007
%   25     1.00000033567660     4.50786280885950e-013    -3.35729598211824e-007
%   26     1.00000016768601     1.12678778792218e-013    -1.67990585574852e-007
%   27     1.00000008408404     2.81185989065354e-014    -8.36019748806373e-008
%   28     1.00000004225090     7.07012543874527e-015    -4.18331397041824e-008
%   29     1.00000002102101     1.78513841016001e-015    -2.12298841719436e-008
%   
% above reults obtaind using a tolerance of
% tol = 100000000*eps; 
% We have convergence if the tolerance is set to a value
% large than 100000000*eps. which is about 1e-8 (the square root of eps).
% Setting a smaller tolerance causes a run-time error due to
% near cancellation of f(x+f) and f(x).
% Consequently, the resulting root is accurate only to about 1e-8.

