octave_packages/m/optimization/fzero.m

   1 ## Copyright (C) 2008-2012 VZLU Prague, a.s.
   2 ##
   3 ## This file is part of Octave.
   4 ##
   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.
   9 ##
  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.
  14 ##
  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/>.
  18 ##
  19 ## Author: Jaroslav Hajek <highegg@gmail.com>
  20
  21 ## -*- texinfo -*-
  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.
  26 ##
  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
  33 ##
  34 ## @example
  35 ## sign (@var{fun}(@var{x0}(1))) * sign (@var{fun}(@var{x0}(2))) <= 0
  36 ## @end example
  37 ##
  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}.
  46 ##
  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:
  50 ##
  51 ## @itemize
  52 ## @item 1
  53 ##  The algorithm converged to a solution.
  54 ##
  55 ## @item 0
  56 ##  Maximum number of iterations or function evaluations has been reached.
  57 ##
  58 ## @item -1
  59 ## The algorithm has been terminated from user output function.
  60 ##
  61 ## @item -5
  62 ## The algorithm may have converged to a singular point.
  63 ## @end itemize
  64 ##
  65 ## @var{output} is a structure containing runtime information about the
  66 ## @code{fzero} algorithm.  Fields in the structure are:
  67 ##
  68 ## @itemize
  69 ## @item iterations
  70 ##  Number of iterations through loop.
  71 ##
  72 ## @item nfev
  73 ##  Number of function evaluations.
  74 ##
  75 ## @item bracketx
  76 ##  A two-element vector with the final bracketing of the zero along the x-axis.
  77 ##
  78 ## @item brackety
  79 ##  A two-element vector with the final bracketing of the zero along the y-axis.
  80 ## @end itemize
  81 ## @seealso{optimset, fsolve}
  82 ## @end deftypefn
  83
  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.
  95
  96 ## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup.
  97 ## PKG_ADD: [~] = __all_opts__ ("fzero");
  98
  99 function [x, fval, info, output] = fzero (fun, x0, options = struct ())
 100
 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");
 105     return;
 106   endif
 107
 108   if (nargin < 2 || nargin > 3)
 109     print_usage ();
 110   endif
 111
 112   if (ischar (fun))
 113     fun = str2func (fun, "global");
 114   endif
 115
 116   ## TODO
 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);
 123
 124   persistent mu = 0.5;
 125
 126   if (funvalchk)
 127     ## Replace fun with a guarded version.
 128     fun = @(x) guarded_eval (fun, x);
 129   endif
 130
 131   ## The default exit flag if exceeded number of iterations.
 132   info = 0;
 133   niter = 0;
 134   nfev = 0;
 135
 136   x = fval = a = fa = b = fb = NaN;
 137   eps = eps (class (x0));
 138
 139   ## Prepare...
 140   a = x0(1);
 141   fa = fun (a);
 142   nfev = 1;
 143   if (length (x0) > 1)
 144     b = x0(2);
 145     fb = fun (b);
 146     nfev += 1;
 147   else
 148     ## Try to get b.
 149     if (a == 0)
 150       aa = 1;
 151     else
 152       aa = a;
 153     endif
 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)
 157         break;
 158       endif
 159     endfor
 160   endif
 161
 162   if (b < a)
 163     u = a;
 164     a = b;
 165     b = u;
 166
 167     fu = fa;
 168     fa = fb;
 169     fb = fu;
 170   endif
 171
 172   if (! (sign (fa) * sign (fb) <= 0))
 173     error ("fzero:bracket", "fzero: not a valid initial bracketing");
 174   endif
 175
 176   slope0 = (fb - fa) / (b - a);
 177
 178   if (fa == 0)
 179     b = a;
 180     fb = fa;
 181   elseif (fb == 0)
 182     a = b;
 183     fa = fb;
 184   endif
 185
 186   itype = 1;
 187
 188   if (abs (fa) < abs (fb))
 189     u = a; fu = fa;
 190   else
 191     u = b; fu = fb;
 192   endif
 193
 194   d = e = u;
 195   fd = fe = fu;
 196   mba = mu*(b - a);
 197   while (niter < maxiter && nfev < maxfev)
 198     switch (itype)
 199     case 1
 200       ## The initial test.
 201       if (b - a <= 2*(2 * abs (u) * eps + tolx))
 202         x = u; fval = fu;
 203         info = 1;
 204         break;
 205       endif
 206       if (abs (fa) <= 1e3*abs (fb) && abs (fb) <= 1e3*abs (fa))
 207         ## Secant step.
 208         c = u - (a - b) / (fa - fb) * fu;
 209       else
 210         ## Bisection step.
 211         c = 0.5*(a + b);
 212       endif
 213       d = u; fd = fu;
 214       itype = 5;
 215     case {2, 3}
 216       l = length (unique ([fa, fb, fd, fe]));
 217       if (l == 4)
 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;
 229       endif
 230       if (l < 4 || sign (c - a) * sign (c - b) > 0)
 231         ## Quadratic interpolation + newton.
 232         a0 = fa;
 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.
 236         c = a - a0/a1;
 237         if (a2 != 0)
 238           c = a - a0/a1;
 239           for i = 1:itype
 240             pc = a0 + (a1 + a2*(c - b))*(c - a);
 241             pdc = a1 + a2*(2*c - a - b);
 242             if (pdc == 0)
 243               c = a - a0/a1;
 244               break;
 245             endif
 246             c -= pc/pdc;
 247           endfor
 248         endif
 249       endif
 250       itype += 1;
 251     case 4
 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))
 256         c = 0.5 * (b + a);
 257       endif
 258       itype = 5;
 259     case 5
 260       ## Bisection step.
 261       c = 0.5 * (b + a);
 262       itype = 2;
 263     endswitch
 264
 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)
 268       c = (a + b)/2;
 269     else
 270       c = max (a + delta, min (b - delta, c));
 271     endif
 272
 273     ## Calculate new point.
 274     x = c;
 275     fval = fc = fun (c);
 276     niter ++; nfev ++;
 277
 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.
 283       fe = fc;
 284     else
 285       e = d; fe = fd;
 286     endif
 287
 288     ## Bracketing.
 289     if (sign (fa) * sign (fc) < 0)
 290       d = b; fd = fb;
 291       b = c; fb = fc;
 292     elseif (sign (fb) * sign (fc) < 0)
 293       d = a; fd = fa;
 294       a = c; fa = fc;
 295     elseif (fc == 0)
 296       a = b = c; fa = fb = fc;
 297       info = 1;
 298       break;
 299     else
 300       ## This should never happen.
 301       error ("fzero:bracket", "fzero: zero point is not bracketed");
 302     endif
 303
 304     ## If there's an output function, use it now.
 305     if (outfcn)
 306       optv.funccount = nfev;
 307       optv.fval = fval;
 308       optv.iteration = niter;
 309       if (outfcn (x, optv, "iter"))
 310         info = -1;
 311         break;
 312       endif
 313     endif
 314
 315     if (abs (fa) < abs (fb))
 316       u = a; fu = fa;
 317     else
 318       u = b; fu = fb;
 319     endif
 320     if (b - a <= 2*(2 * abs (u) * eps + tolx))
 321       info = 1;
 322       break;
 323     endif
 324
 325     ## Skip bisection step if successful reduction.
 326     if (itype == 5 && (b - a) <= mba)
 327       itype = 2;
 328     endif
 329     if (itype == 2)
 330       mba = mu * (b - a);
 331     endif
 332   endwhile
 333
 334   ## Check solution for a singularity by examining slope
 335   if (info == 1)
 336     if ((b - a) != 0 && abs ((fb - fa)/(b - a) / slope0) > max (1e6, 0.5/(eps+tolx)))
 337       info = -5;
 338     endif
 339   endif
 340
 341   output.iterations = niter;
 342   output.funcCount = nfev;
 343   output.bracketx = [a, b];
 344   output.brackety = [fa, fb];
 345
 346 endfunction
 347
 348 ## An assistant function that evaluates a function handle and checks for
 349 ## bad results.
 350 function fx = guarded_eval (fun, x)
 351   fx = fun (x);
 352   fx = fx(1);
 353   if (! isreal (fx))
 354     error ("fzero:notreal", "fzero: non-real value encountered");
 355   elseif (isnan (fx))
 356     error ("fzero:isnan", "fzero: NaN value encountered");
 357   endif
 358 endfunction
 359
 360 %!shared opt0
 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)