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 
.
, 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
3 件のコメント
Walter Roberson
2019 年 4 月 2 日
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
2019 年 4 月 2 日
@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
2019 年 4 月 4 日
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
採用された回答
Simply use the debugger to test, what's going on:
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);
2 件のコメント
Stephen23
2019 年 4 月 4 日
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
2019 年 4 月 8 日
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.
その他の回答 (1 件)
Agostino Dorano
2019 年 4 月 4 日
1 件のコメント
Walter Roberson
2019 年 4 月 4 日
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.
カテゴリ
ヘルプ センター および File Exchange で Live Scripts and Functions についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
