%{
Roots (both single and multiple) of an equation f(x)=0 
using five methods
Methods: Bisectional, False position, Secant, Brute force
also using Matlab function fzero
February 2022 by AG.
%}

clear       % Clears command window
clc         % Clears command history
clearvars   % Clear variables from memory
clf('reset')% Removes anything in the figure window before simulation. 

tic         % set timer for the program

% Initial data ==============================
a = -10.0;       % left-end point
b =   2.0;       % right-end point
x0=   1.0;       % guess point for Secant and Matlab (root)
eps = 1.0e-8;   % desired uncertainity of the root
Nint= 100;      % number of subintervals for brute force method

% Name of the integrand @f
f0=@f;
%============================================

% Estimate number of steps for given tolerance
Nsteps =  ceil((1.0/log(2))*log((b-a)/eps));

format1 = ' %3i  %10.6f  %10.6f ';
format2 = '      %10.6f  %10.6f ';
format3 = ' %3i  %10.6f  %10.6f \n';

fprintf(' Root(s) of f(x) \n')
fprintf(' Interval [%4.1f,%4.1f] \n',a,b)
fprintf(' Guess Secant and Matlab x0 =%4.1f \n',x0)
fprintf(' Tolerance = %6.2e \n',eps)
% fprintf(' Steps for bisectional method %3i \n',Nsteps)
fprintf(' iter    root       f(root)\n')

% Step 0: plot the function
xf = a: 0.01: b;
yf = f(xf);
yzero = 0.0*xf;
plot(xf,yf,xf,yzero,'LineWidth',1.5)
xlabel('x');
ylabel('f(x)');
grid;

% 1: Bisectional method
[root1,flag] = root_bs(f0, a, b, eps);
if flag ~= 0
    fprintf(format1,flag, root1,f0(root1))
    fprintf('  Bisectional  \n')
else
    fprintf(' Bisectional:     No root found \n')
end

% 2: False position method
[root2,flag] = root_fp(f0, a, b, eps);
if flag ~= 0
    fprintf(format1,flag, root2,f0(root2))
    fprintf('  False position \n')
else
    fprintf(' False position:  No root found \n')
end

% 3: method of secants
x1 = x0;
x2 = x1+0.1;
[root3,flag] = root_sm(f0, a, b, x0, eps);
if flag ~= 0
    fprintf(format1,flag, root3,f0(root3))
    fprintf('  Secant \n')
else
    fprintf('Secant method - No root found \n')
end

% 4: Matlab function fzero
root4 = fzero(f0,x0);
    fprintf(format2, root4,f0(root4))
    fprintf('  MatLab \n')

% 5: Broot force method
 [Broots, nroots] = Bforce(f0, a, b, eps, Nint);
 fprintf('  \n  Brute force roots\n ')
%print solutions
if(nroots == 0)
    fprintf (' No roots found \n')
    return
else
    fprintf('  n     root       f(root)\n')     
for i=1:nroots
  fprintf(format3,i, Broots(i), f0(Broots(i)))
end
end
fprintf('\n')

toc
timeE = ceil(toc);
%fprintf('\n Elapsed time = %4i sec \n',timeE)
%==============================

function f=f(x)
% use [-10,+10]
%f = x + cos(x);

% use [-2,2], x0=1
%f = x-cos(x);                      % easy root at 0.739           

% use [-2,2] x0=1
%f = x.^3 - 2*x.^2 + 1.5*x - 1/3;  % false position fails (too slow)

% use [-2,2] x0=1
%f = 4*x.^4 -6*x.^2 -11/4;

% use [-2,2] x0=1
%f = exp(x) - 2.*sin(x) - 2; 

% very fancy
% use [-10,2]
f = exp(x).*log(x.^2) - x.*cos(x);

%f = exp(x) - x.*sin(x*pi/2.0);
end
%==============================

function [zero,flag] = root_bs(f, a, b, eps)
%{
!============================================================
! Solutions of equation f(x)=0 on [x1,x2] interval
! Method: Bisectional (closed domain) (a single root)
!------------------------------------------------------------
! input ...
! f   - function - evaluates f(x) for any x in [x1,x2]
! a   - left endpoint of initial interval
! b   - right endpoint of initial interval
! eps - desired uncertainity of the root as |b-a|<eps
! output ...
! zero  - root of the equation f(x)=0
! flag  - indicator of success
!         >0 - a single root found, flag=number of iterations
!          0 - no solutions for the bisectional method
!         <0 - not a root but singularity, flag=number of iterations
!
! Comments: Function f(x) has to change sign between x1 and x2
!           Max number of iterations - 999 (accuracy (b-a)/2*999
%}

iter = 999;     % Max number of iterations
flag = 0;       % set success to 0;
zero = 99.99;   % if there is no root
% check the bisection condition
    if f(a)*f(b) > 0.0
        return
    end
% finding a single root
    xl = a;
    xr = b;   
       
while abs(xr - xl) > eps
     flag = flag + 1;
     x0 = (xr + xl)/2.0;
     if (f(xl) * f(x0)) <= 0.0
        xr = x0;
        else
        xl = x0;
     end
     if(flag >= iter) 
         break
     end
end
    zero = x0;
end
% end of root_bs



function [zero,flag] = root_fp(f, a, b, eps)
%{
!============================================================
! Solutions of equation f(x)=0 on [x1,x2] interval
! Method: False position (closed domain) (a single root)
!------------------------------------------------------------
The same description as root_bs
%}

iter = 999;     % Max number of iterations
flag = 0;       % initial values for iterations
zero = 99.99;   % if there is no root
% check the bisection condition
    if f(a)*f(b) > 0.0
        return
    end
% finding a single root
    xl = a;
    xr = b;
    x0 = (xr + xl)/2.0;
       
while abs(f(x0)) > eps
     flag = flag + 1;
     x0 = xr - (f(xr)*(xr - xl))/(f(xr)-f(xl));     
     if (f(xl) * f(x0)) <= 0.0
        xr = x0;
        else
        xl = x0;
     end
     if(flag >= iter) 
         break
     end
end
zero = x0;
end


function [zero,flag] = root_sm(f, x1, x2, x0, eps)
%{
!============================================================
! Solutions of equation f(x)=0 on [x1,x2] interval
! Method: Secant method (a single root)
!------------------------------------------------------------
! input ...
! f   - function - evaluates f(x) for any x in [x1,x2]
! x1 - initial point
! x2 - second initial point
! eps - desired uncertainity of the root as |b-a|<eps
! output ...
! zero  - root of the equation f(x)=0
! flag  - indicator of success
!         >0 - a single root found, flag=number of iterations
!          0 - no solutions for the bisectional method
!         <0 - not a root but singularity, flag=number of iterations
!
! Comments: Function f(x) has to change sign between x1 and x2
!           Max number of iterations - 999 (accuracy (b-a)/(2*999)
%}

iter = 999;
flag   = 0;
while (abs(x2 - x1)) > eps
    flag = flag + 1;
    if flag == 1
        x3 = x0;
    else
        x3 = x2 - (f(x2)*(x2-x1))/(f(x2)-f(x1));
    end
    x1 = x2;
    x2 = x3;
    if(flag >= iter) 
        break
    end %end if
    flag = flag;
    if (flag == iter) 
        flag = 0;
    end %end if
end % end while
zero = x3;
end

% Brute force method
function [Broots, nroots] = Bforce(f, x1, x2, eps, n)

Broots(1)=0.0;
iter = 200;
% initialization
dx = (x2-x1)/n;
nroots = 0;

% loop over subintervals
for j=1:n
  a = x1  + (j-1)*dx;
  b = a + dx;
% check the closed domain condition f(a)*f(b)<0
  if(f(a)*f(b)>0) 
      continue
  end

% Iterative refining the solution 
  for i=1:iter
      c=(b+a)/2.0;
    if(f(a)*f(c) <= 0.0)
      b = c;
    else
      a=c;
    end
%   condition(s) to stop iterations)
    if(abs(b-a)<= eps)
        break
    end
  end
% check if it is a root or singularity
  root = (b+a)/2.0;
  if (abs(f(root)) < 1.0) 
    nroots = nroots+1;
    Broots(nroots)=root;
  end
end
end