HELP! Bisection method fails to found zero

I tried to calculate the zero for the function , but the zero isn't found in . I have the same problem for other function like .

I paste here my code (part of my function) :

TOLF = eps;
nit= 0; % number of iterations
m = xo(1)+0.5*diff(xo); % midpoint
while(abs(diff(xo))> p.Results.TOL*max(abs(xo)) && ...
        abs(fun(m))>TOLF && nit<p.Results.NMAX) % stopping condition
    if(fun(m)*fun(xo(1))<0)
        xo(2)= m;
    else
        xo(1)= m;
    end
    nit=nit+1;
    m = xo(1)+0.5*diff(xo);
end

[Walter Roberson]

Your second case, x^4, has the problem that it has no zero crossing. Bisection methods can only be used if you have the endpoints of an interval in which there is at least one zero crossing.

[Jan]

@Agostino Dorano: Please post code, which we can run. What are the initial values of xo? What is p.Results.TOL and p.Results.NMAX? How do you define fun? What exactly do you observe? "zero isn't found" does not reveal the details.

[Stephen23]

Agostino Dorano's "Answer" moved here:

% Function that calculates the zero of functions.
% This function implements the bisection algorithm to found zero of a function.
function [x, output] = fzer0(fun,xo,varargin)

% CONTROLLO NUMERI DI PARAMETRI
narginchk(2, 5);

% CREO L'INPUT PARSER
p = inputParser;

% CONTROLLO SUI TIPI DEI VALORI DI INGRESSO
addRequired (p, 'fun', @(x) assert(isa(x,'function_handle'),'First argoument must be a function handle of a in R function'));
addRequired (p, 'xo' , @(x) assert( isreal(x) && isnumeric(x) && isvector(x) && (length(x)==2) && all(isfinite(x))  ,'xo must be a vector of 2 real numbers')); % minore di zero perche deve avere proprio uno zero crossing
% SOSTITUIRE ADDPARAMETER CON ADDOPTIONAL PER RENDERE POSIZIONALE
addParameter(p,'TOL',eps,@(x) assert(isreal(x) && isnumeric(x) && all(isfinite(x)) && x>0 && x>=eps, 'TOL value must be a real value greather than 0 and eps and lower than inf'));
addParameter(p, 'NMAX', 500,  @(x) assert(isreal(x) && isnumeric(x) && isfinite(x) && x>0 && floor(x)==x,'NMAX value must be a integer number greather than 0'));
addParameter(p, 'g','N', @(x) assert(ismember(x, {'Y','N'}),'g charater must be Y or N. See doc fzer0') );
parse (p, fun,xo, varargin{:});

xo = p.Results.xo; % array of two elemens that rapresents the interval extreme

% CONTROLLO DI APPLICABILITA'
if (fun(xo(1))*fun(xo(2))>=0)
    error("zero or more zero in xo");
end

TOLF = eps; 
nit= 0;
m = xo(1)+0.5*diff(xo);

%TOLX = p.Results.TOL*max(abs(xo)); sostituisco la sera del due aprile col
% codice riportato sotto===> dubbio: possiamo evitarlo? anzi, per il
% momento lascio tolx definito come : TOLX = p.Results.TOL*max(abs(xo));

% SET TOLX IN MANIERA SICURA
if(max(abs(xo))>realmin/p.Results.TOL)
    TOLX = p.Results.TOL * max(abs(xo));
else
    TOLX = realmin;
end

%%%%%%%% BISOGNA CONTROLLARE SE UNO DEI DUE ESTREMI è UNO ZERO => SE TALE
%%%%%%%% CONDIZIONE è VERA DEVO SETTARE LA VARIABILE DI RITORNO E TERMINARE
%%%%%%%% => ma tutto cio non è corretto

while(abs(diff(xo))>= TOLX && abs(fun(m))>=TOLF && nit<p.Results.NMAX) 
    if(fun(m)*fun(xo(1))<0)
        xo(2)= m;
    else
        xo(1)= m;
    end
    nit=nit+1;
    m = xo(1)+0.5*diff(xo);
    % CONTROLLO DI ROBUSTEZZA PER IL SET DI TOLX
    if(max(abs(xo))>realmin/p.Results.TOL)
        TOLX = p.Results.TOL * max(abs(xo));
    else
        TOLX = realmin;
    end
    
end

if(abs(diff(xo))> p.Results.TOL*max(abs(xo)) && abs(fun(m))>TOLF && nit==p.Results.NMAX)
    warning("limit reached... calculation stopped!");
end
% return var
nargoutchk(0,2);
x = m;
if (nargout ==2 )
    output.fx = fun(x);
    output.niter = nit;
end

end

sorry for the late reply. I paste my complete code here. I want ask you also if to utilize InputParser is a good method to validate inputs for my function (better in performance(execution time) than nargin + if/ switch or exist) .

For "zero not found " i paste you the result of mine program comparated with fzero of matlab. them looks like pretty different . Thanks a lot.

>> fzer0(@(x) x^3, [-2 1]) % this is my function

ans =

     3.814697265625000e-06

>> fzero(@(x) x^3, [-2 1]) % default matlab fzero

ans =

     1.001465536366818e-16 
     
     % both functions were calculted with tolerance equal to eps 

Answer 1

[Jan]

Simply use the debugger to test, what's going on:

https://www.mathworks.com/help/matlab/matlab_prog/debugging-process-and-features.html

You can check programmatically also, why the loop is stoped:

ready = false;
while ~ready
    ...
    readyTOL = (abs(diff(xo)) >= TOLX;
    readyABS = abs(fun(m)) >= TOLF;
    readyNIT = nit<p.Results.NMAX;
    ready    = (readyTOL || readyABS || readINT);
end
output.stopped = [readyTOL, readyABS, readINT];
...

Now you can see, which condition has stopped the loop.

By the way, this test is not useful:

if (abs(diff(xo))> p.Results.TOL*max(abs(xo)) && ...
        abs(fun(m))>TOLF && ...
        nit==p.Results.NMAX)
    warning("limit reached... calculation stopped!");
end

It is very unlikly, that the maximum iteration number was reached and the absolute tolerance and the relative tolerance for the interval length at the same time.

I assume you want this:

if nit == p.Results.NMAX
    warning("Maximum number of iterations reached!");
end

This warns, if the found solution is not accurate, because the process stopped by the iteration counter.

What about replacing:

    if(max(abs(xo))>realmin/p.Results.TOL)
        TOLX = p.Results.TOL * max(abs(xo));
    else
        TOLX = realmin;
    end

by

TOLX = max(p.Results.TOL * mmax(abs(xo)), realmin);

[Stephen23]

Agostino Dorano's "Answer" moved here:

@Jan Thanks for your reply, i've updated my code following yours advices. I post here my updated code.

function [x, output] = fzer0v3(fun,xo,varargin)

% CONTROLLO NUMERI DI PARAMETRI -- check inputs number
narginchk(2, 5);

% INPUT PARSER
p = inputParser;

% ERROR INPUTS CHECK => strongly wrong
addRequired (p, 'fun', @(x) assert(isa(x,'function_handle'),'First argoument must be a function handle of a in R function'));
addRequired (p, 'xo' , @CheckXoErr); 
addOptional(p,'TOL',eps,@CheckTolErr);
addOptional(p, 'NMAX', 500,@CheckNmaxErr);
addOptional(p, 'g','N', @(x) assert(ismember(x, {'Y','N'}),'g charater must be Y or N. See doc fzer0') );
parse (p, fun,xo, varargin{:});

% tolleranza aggiungere controllo <=10^-1 => controllo solo che sia minore di 1 % settare tolleranza di default nel caso sia minore di 0 ed eps => fatto % controllo nmax e settaggio a 500 nel caso sia minore uguale di zero => % fatto

% applicability check
if (fun(xo(1))*fun(xo(2))>0)
    error("zero or more zero in xo");
end

xo = p.Results.xo;
TOL = p.Results.TOL;
NMAX = p.Results.NMAX;

% WARNING AND CORRECTION OF TOLERANCE AND NMAX (if inputs are "little wrong")
switch nargin
    case 5; [TOL, NMAX]= CheckNMAXWarn(TOL,NMAX);
    case 4; [TOL, NMAX]= CheckNMAXWarn(TOL,NMAX);
    case 3; TOL = CheckTOLWarn(TOL);
    otherwise
end

% SET OF SOME VARIABLES
TOLF = eps; 
nit= 0;
m = xo(1)+0.5*diff(xo);
TOLX = max(TOL * max(abs(xo)), realmin);
%%%%%%%% BISOGNA CONTROLLARE SE UNO DEI DUE ESTREMI è UNO ZERO => SE TALE CONDIZIONE è VERA DEVO SETTARE LA VARIABILE DI RITORNO E TERMINARE=> ma tutto cio non è corretto
while(abs(diff(xo))>= TOLX && abs(fun(m))>=TOLF && nit<NMAX) 
    if(fun(m)*fun(xo(1))<0)
        xo(2)= m;
    else
        xo(1)= m;
    end
    nit=nit+1;
    m = xo(1)+0.5*diff(xo);
    TOLX = max(TOL*max(abs(xo)), realmin);
end

% if(abs(diff(xo))> p.Results.TOL*max(abs(xo)) && abs(fun(m))>TOLF && nit==p.Results.NMAX)
if(nit==p.Results.NMAX)
    warning("limit reached... calculation stopped!");
end

if p.Results.g == 'Y'
    plotFunction(fun, p.Results.xo, m);
end
nargoutchk(0,2);
x = m;
if (nargout ==2 )
    output.fx = fun(x);
    output.niter = nit;
end
end

function [TOL, NMAX]=CheckNMAXWarn(tol,nmax) % FUNCTION THAT CONTROL AND CORRECT NMAX VALUE. IT CALL CHECKTOLWARN FUNCTION
NMAX = nmax;
if nmax<=0 
    warning("Specified NMAX is too low. Set to default value.");
    NMAX = 500;
end
TOL = CheckTOLWarn(tol);
end

function TOL = CheckTOLWarn(tol) % FUNCTION THAT CONTROL AND CORRECT TOLERANCE VALUE.
TOL = tol;
if tol < eps || tol >=1
    warning("tolerance value is wrong, see doc fzer0. Set to default value.");
    TOL = eps; 
end
end

function TF = CheckXoErr(xo)
TF = false;
if (~isreal(xo)||~isnumeric(xo))
    error('xo must be composed by only real number.');
elseif ~all(isfinite(xo))
    error("xo elements can't be + - Inf or NaN.");
elseif (~isvector(xo)||~(length(xo)==2))
    error('xo must be a vector of 2 elements.');
else
    TF = true;
end
end

%%%%% support function

function TF = CheckTolErr(tol)
TF = false;
if (~isreal(tol)||~isnumeric(tol)||~(length(tol)==1))
    error('TOL must be a real number.');
elseif ~isfinite(tol)
    error("TOL can't be + - Inf or NaN.");
else
    TF = true;
end
end

function TF = CheckNmaxErr(nmax)
TF = false;
if (~isreal(nmax)||~isnumeric(nmax)||~(length(nmax)==1)||~isfinite(nmax))
    error('NMAX must be a real number.');
elseif ~(floor(nmax)==nmax)
    error("NMAX must be an integer value");
else
    TF = true;
end
end

Now my question is: if I have constraints on the execution time , is this the best way to check function's inputs (errors and warnings)?

I searched in the matlab doc and on the web, but i didn't found nothing of interesting.

I'm curious and I'd like to know if there are better way to implements these check.

I'll wait for your reply.

[Jan]

It depends, on what you call "better". The inputParser is not highly efficient, but exhaustively tested. If you write this set of tests by your own, you can trust it only with some unit-tests which prove the correct working and rejecting of bad inputs. So maybe te code runs some microseconds faster, therefore you need 1 hour for programming and testing.

Answer 2

[Agostino Dorano]

@Jan Thanks for your reply, i've updated my code following yours advices. I post here my updated code.

function [x, output] = fzer0v3(fun,xo,varargin)

% CONTROLLO NUMERI DI PARAMETRI -- check inputs number
narginchk(2, 5);

% INPUT PARSER
p = inputParser;

% ERROR INPUTS CHECK => strongly wrong
addRequired (p, 'fun', @(x) assert(isa(x,'function_handle'),'First argoument must be a function handle of a in R function'));
addRequired (p, 'xo' , @CheckXoErr); 
addOptional(p,'TOL',eps,@CheckTolErr);
addOptional(p, 'NMAX', 500,@CheckNmaxErr);
addOptional(p, 'g','N', @(x) assert(ismember(x, {'Y','N'}),'g charater must be Y or N. See doc fzer0') );
parse (p, fun,xo, varargin{:});

% tolleranza aggiungere controllo <=10^-1 => controllo solo che sia minore di 1 % settare tolleranza di default nel caso sia minore di 0 ed eps => fatto % controllo nmax e settaggio a 500 nel caso sia minore uguale di zero => % fatto

% applicability check
if (fun(xo(1))*fun(xo(2))>0)
    error("zero or more zero in xo");
end

xo = p.Results.xo;
TOL = p.Results.TOL;
NMAX = p.Results.NMAX;

% WARNING AND CORRECTION OF TOLERANCE AND NMAX (if inputs are "little wrong")
switch nargin
    case 5; [TOL, NMAX]= CheckNMAXWarn(TOL,NMAX);
    case 4; [TOL, NMAX]= CheckNMAXWarn(TOL,NMAX);
    case 3; TOL = CheckTOLWarn(TOL);
    otherwise
end

% SET OF SOME VARIABLES
TOLF = eps; 
nit= 0;
m = xo(1)+0.5*diff(xo);
TOLX = max(TOL * max(abs(xo)), realmin);
%%%%%%%% BISOGNA CONTROLLARE SE UNO DEI DUE ESTREMI è UNO ZERO => SE TALE CONDIZIONE è VERA DEVO SETTARE LA VARIABILE DI RITORNO E TERMINARE=> ma tutto cio non è corretto
while(abs(diff(xo))>= TOLX && abs(fun(m))>=TOLF && nit<NMAX) 
    if(fun(m)*fun(xo(1))<0)
        xo(2)= m;
    else
        xo(1)= m;
    end
    nit=nit+1;
    m = xo(1)+0.5*diff(xo);
    TOLX = max(TOL*max(abs(xo)), realmin);
end

% if(abs(diff(xo))> p.Results.TOL*max(abs(xo)) && abs(fun(m))>TOLF && nit==p.Results.NMAX)
if(nit==p.Results.NMAX)
    warning("limit reached... calculation stopped!");
end

if p.Results.g == 'Y'
    plotFunction(fun, p.Results.xo, m);
end
nargoutchk(0,2);
x = m;
if (nargout ==2 )
    output.fx = fun(x);
    output.niter = nit;
end
end

function [TOL, NMAX]=CheckNMAXWarn(tol,nmax) % FUNCTION THAT CONTROL AND CORRECT NMAX VALUE. IT CALL CHECKTOLWARN FUNCTION
NMAX = nmax;
if nmax<=0 
    warning("Specified NMAX is too low. Set to default value.");
    NMAX = 500;
end
TOL = CheckTOLWarn(tol);
end

function TOL = CheckTOLWarn(tol) % FUNCTION THAT CONTROL AND CORRECT TOLERANCE VALUE.
TOL = tol;
if tol < eps || tol >=1
    warning("tolerance value is wrong, see doc fzer0. Set to default value.");
    TOL = eps; 
end
end

function TF = CheckXoErr(xo)
TF = false;
if (~isreal(xo)||~isnumeric(xo))
    error('xo must be composed by only real number.');
elseif ~all(isfinite(xo))
    error("xo elements can't be + - Inf or NaN.");
elseif (~isvector(xo)||~(length(xo)==2))
    error('xo must be a vector of 2 elements.');
else
    TF = true;
end
end

%%%%% support function

function TF = CheckTolErr(tol)
TF = false;
if (~isreal(tol)||~isnumeric(tol)||~(length(tol)==1))
    error('TOL must be a real number.');
elseif ~isfinite(tol)
    error("TOL can't be + - Inf or NaN.");
else
    TF = true;
end
end

function TF = CheckNmaxErr(nmax)
TF = false;
if (~isreal(nmax)||~isnumeric(nmax)||~(length(nmax)==1)||~isfinite(nmax))
    error('NMAX must be a real number.');
elseif ~(floor(nmax)==nmax)
    error("NMAX must be an integer value");
else
    TF = true;
end
end

Now my question is: if I have constraints on the execution time , is this the best way to check function's inputs (errors and warnings)?

I searched in the matlab doc and on the web, but i didn't found nothing of interesting.

I'm curious and I'd like to know if there are better way to implements these check.

I'll wait for your reply.

[Walter Roberson]

If your code must be as efficient as possible, then No, using addRequired and related routines is not as efficient as you could possibly get. They involve multiple function calls. The calls do not take long, but it is more efficient at execution time to "hard code" the work to be done rather than to call a function to do the work. You can probably save a good third of a microsecond of execution time by doing the checks yourself.