1 ## Copyright (C) 2008-2012 VZLU Prague, a.s.
3 ## This file is part of Octave.
5 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by
7 ## the Free Software Foundation; either version 3 of the License, or (at
8 ## your option) any later version.
10 ## Octave is distributed in the hope that it will be useful, but
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 ## General Public License for more details.
15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>.
19 ## Author: Jaroslav Hajek <highegg@gmail.com>
22 ## @deftypefn {Function File} {} fzero (@var{fun}, @var{x0})
23 ## @deftypefnx {Function File} {} fzero (@var{fun}, @var{x0}, @var{options})
24 ## @deftypefnx {Function File} {[@var{x}, @var{fval}, @var{info}, @var{output}] =} fzero (@dots{})
25 ## Find a zero of a univariate function.
27 ## @var{fun} is a function handle, inline function, or string
28 ## containing the name of the function to evaluate.
29 ## @var{x0} should be a two-element vector specifying two points which
30 ## bracket a zero. In other words, there must be a change in sign of the
31 ## function between @var{x0}(1) and @var{x0}(2). More mathematically, the
32 ## following must hold
35 ## sign (@var{fun}(@var{x0}(1))) * sign (@var{fun}(@var{x0}(2))) <= 0
38 ## If @var{x0} is a single scalar then several nearby and distant
39 ## values are probed in an attempt to obtain a valid bracketing. If this
40 ## is not successful, the function fails.
41 ## @var{options} is a structure specifying additional options.
42 ## Currently, @code{fzero}
43 ## recognizes these options: @code{"FunValCheck"}, @code{"OutputFcn"},
44 ## @code{"TolX"}, @code{"MaxIter"}, @code{"MaxFunEvals"}.
45 ## For a description of these options, see @ref{doc-optimset,,optimset}.
47 ## On exit, the function returns @var{x}, the approximate zero point
48 ## and @var{fval}, the function value thereof.
49 ## @var{info} is an exit flag that can have these values:
53 ## The algorithm converged to a solution.
56 ## Maximum number of iterations or function evaluations has been reached.
59 ## The algorithm has been terminated from user output function.
62 ## The algorithm may have converged to a singular point.
65 ## @var{output} is a structure containing runtime information about the
66 ## @code{fzero} algorithm. Fields in the structure are:
70 ## Number of iterations through loop.
73 ## Number of function evaluations.
76 ## A two-element vector with the final bracketing of the zero along the x-axis.
79 ## A two-element vector with the final bracketing of the zero along the y-axis.
81 ## @seealso{optimset, fsolve}
84 ## This is essentially the ACM algorithm 748: Enclosing Zeros of
85 ## Continuous Functions due to Alefeld, Potra and Shi, ACM Transactions
86 ## on Mathematical Software, Vol. 21, No. 3, September 1995. Although
87 ## the workflow should be the same, the structure of the algorithm has
88 ## been transformed non-trivially; instead of the authors' approach of
89 ## sequentially calling building blocks subprograms we implement here a
90 ## FSM version using one interior point determination and one bracketing
91 ## per iteration, thus reducing the number of temporary variables and
92 ## simplifying the algorithm structure. Further, this approach reduces
93 ## the need for external functions and error handling. The algorithm has
94 ## also been slightly modified.
96 ## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup.
97 ## PKG_ADD: [~] = __all_opts__ ("fzero");
99 function [x, fval, info, output] = fzero (fun, x0, options = struct ())
101 ## Get default options if requested.
102 if (nargin == 1 && ischar (fun) && strcmp (fun, 'defaults'))
103 x = optimset ("MaxIter", Inf, "MaxFunEvals", Inf, "TolX", 1e-8, \
104 "OutputFcn", [], "FunValCheck", "off");
108 if (nargin < 2 || nargin > 3)
113 fun = str2func (fun, "global");
117 ## displev = optimget (options, "Display", "notify");
118 funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
119 outfcn = optimget (options, "OutputFcn");
120 tolx = optimget (options, "TolX", 1e-8);
121 maxiter = optimget (options, "MaxIter", Inf);
122 maxfev = optimget (options, "MaxFunEvals", Inf);
127 ## Replace fun with a guarded version.
128 fun = @(x) guarded_eval (fun, x);
131 ## The default exit flag if exceeded number of iterations.
136 x = fval = a = fa = b = fb = NaN;
137 eps = eps (class (x0));
154 for b = [0.9*aa, 1.1*aa, aa-1, aa+1, 0.5*aa 1.5*aa, -aa, 2*aa, -10*aa, 10*aa]
155 fb = fun (b); nfev += 1;
156 if (sign (fa) * sign (fb) <= 0)
172 if (! (sign (fa) * sign (fb) <= 0))
173 error ("fzero:bracket", "fzero: not a valid initial bracketing");
176 slope0 = (fb - fa) / (b - a);
188 if (abs (fa) < abs (fb))
197 while (niter < maxiter && nfev < maxfev)
201 if (b - a <= 2*(2 * abs (u) * eps + tolx))
206 if (abs (fa) <= 1e3*abs (fb) && abs (fb) <= 1e3*abs (fa))
208 c = u - (a - b) / (fa - fb) * fu;
216 l = length (unique ([fa, fb, fd, fe]));
218 ## Inverse cubic interpolation.
219 q11 = (d - e) * fd / (fe - fd);
220 q21 = (b - d) * fb / (fd - fb);
221 q31 = (a - b) * fa / (fb - fa);
222 d21 = (b - d) * fd / (fd - fb);
223 d31 = (a - b) * fb / (fb - fa);
224 q22 = (d21 - q11) * fb / (fe - fb);
225 q32 = (d31 - q21) * fa / (fd - fa);
226 d32 = (d31 - q21) * fd / (fd - fa);
227 q33 = (d32 - q22) * fa / (fe - fa);
228 c = a + q31 + q32 + q33;
230 if (l < 4 || sign (c - a) * sign (c - b) > 0)
231 ## Quadratic interpolation + newton.
233 a1 = (fb - fa)/(b - a);
234 a2 = ((fd - fb)/(d - b) - a1) / (d - a);
235 ## Modification 1: this is simpler and does not seem to be worse.
240 pc = a0 + (a1 + a2*(c - b))*(c - a);
241 pdc = a1 + a2*(2*c - a - b);
252 ## Double secant step.
253 c = u - 2*(b - a)/(fb - fa)*fu;
254 ## Bisect if too far.
255 if (abs (c - u) > 0.5*(b - a))
265 ## Don't let c come too close to a or b.
266 delta = 2*0.7*(2 * abs (u) * eps + tolx);
267 if ((b - a) <= 2*delta)
270 c = max (a + delta, min (b - delta, c));
273 ## Calculate new point.
278 ## Modification 2: skip inverse cubic interpolation if
279 ## nonmonotonicity is detected.
280 if (sign (fc - fa) * sign (fc - fb) >= 0)
281 ## The new point broke monotonicity.
282 ## Disable inverse cubic.
289 if (sign (fa) * sign (fc) < 0)
292 elseif (sign (fb) * sign (fc) < 0)
296 a = b = c; fa = fb = fc;
300 ## This should never happen.
301 error ("fzero:bracket", "fzero: zero point is not bracketed");
304 ## If there's an output function, use it now.
306 optv.funccount = nfev;
308 optv.iteration = niter;
309 if (outfcn (x, optv, "iter"))
315 if (abs (fa) < abs (fb))
320 if (b - a <= 2*(2 * abs (u) * eps + tolx))
325 ## Skip bisection step if successful reduction.
326 if (itype == 5 && (b - a) <= mba)
334 ## Check solution for a singularity by examining slope
336 if ((b - a) != 0 && abs ((fb - fa)/(b - a) / slope0) > max (1e6, 0.5/(eps+tolx)))
341 output.iterations = niter;
342 output.funcCount = nfev;
343 output.bracketx = [a, b];
344 output.brackety = [fa, fb];
348 ## An assistant function that evaluates a function handle and checks for
350 function fx = guarded_eval (fun, x)
354 error ("fzero:notreal", "fzero: non-real value encountered");
356 error ("fzero:isnan", "fzero: NaN value encountered");
361 %! opt0 = optimset ("tolx", 0);
362 %!assert(fzero(@cos, [0, 3], opt0), pi/2, 10*eps)
363 %!assert(fzero(@(x) x^(1/3) - 1e-8, [0,1], opt0), 1e-24, 1e-22*eps)