Mercurial > octave-nkf
view scripts/optimization/fzero.m @ 17338:1c89599167a6
maint: End m-files with 1 blank line.
Simplifies automated grammarchecking script.
* scripts/@ftp/ascii.m, scripts/@ftp/binary.m, scripts/@ftp/cd.m,
scripts/@ftp/close.m, scripts/@ftp/delete.m, scripts/@ftp/dir.m,
scripts/@ftp/display.m, scripts/@ftp/ftp.m, scripts/@ftp/loadobj.m,
scripts/@ftp/mget.m, scripts/@ftp/mkdir.m, scripts/@ftp/mput.m,
scripts/@ftp/rename.m, scripts/@ftp/rmdir.m, scripts/@ftp/saveobj.m,
scripts/audio/lin2mu.m, scripts/audio/loadaudio.m, scripts/audio/mu2lin.m,
scripts/audio/record.m, scripts/audio/saveaudio.m, scripts/audio/setaudio.m,
scripts/deprecated/__error_text__.m, scripts/deprecated/cut.m,
scripts/deprecated/error_text.m, scripts/deprecated/isstr.m,
scripts/deprecated/polyderiv.m, scripts/deprecated/studentize.m,
scripts/deprecated/sylvester_matrix.m, scripts/general/bicubic.m,
scripts/general/celldisp.m, scripts/general/colon.m,
scripts/general/cplxpair.m, scripts/general/del2.m, scripts/general/display.m,
scripts/general/isdir.m, scripts/general/isequaln.m, scripts/general/loadobj.m,
scripts/general/private/__isequal__.m, scripts/general/private/__splinen__.m,
scripts/general/profexplore.m, scripts/general/quadgk.m,
scripts/general/randi.m, scripts/general/repmat.m, scripts/general/saveobj.m,
scripts/geometry/delaunay.m, scripts/help/__unimplemented__.m,
scripts/help/doc_cache_create.m, scripts/help/get_first_help_sentence.m,
scripts/help/help.m, scripts/help/print_usage.m,
scripts/help/private/__additional_help_message__.m,
scripts/help/private/__strip_html_tags__.m, scripts/help/type.m,
scripts/image/imfinfo.m, scripts/image/imformats.m, scripts/image/imread.m,
scripts/image/imwrite.m, scripts/image/private/__imfinfo__.m,
scripts/image/private/__imread__.m, scripts/image/private/__imwrite__.m,
scripts/image/private/imageIO.m, scripts/image/private/imwrite_filename.m,
scripts/image/private/ind2x.m, scripts/io/beep.m, scripts/io/strread.m,
scripts/io/textread.m, scripts/io/textscan.m, scripts/linear-algebra/krylov.m,
scripts/linear-algebra/subspace.m, scripts/miscellaneous/bug_report.m,
scripts/miscellaneous/bunzip2.m, scripts/miscellaneous/cast.m,
scripts/miscellaneous/copyfile.m, scripts/miscellaneous/debug.m,
scripts/miscellaneous/dir.m, scripts/miscellaneous/dump_prefs.m,
scripts/miscellaneous/error_ids.m, scripts/miscellaneous/fileattrib.m,
scripts/miscellaneous/gunzip.m, scripts/miscellaneous/isdeployed.m,
scripts/miscellaneous/ismac.m, scripts/miscellaneous/mex.m,
scripts/miscellaneous/mexext.m, scripts/miscellaneous/mkoctfile.m,
scripts/miscellaneous/movefile.m, scripts/miscellaneous/namelengthmax.m,
scripts/miscellaneous/news.m, scripts/miscellaneous/pack.m,
scripts/miscellaneous/perl.m,
scripts/miscellaneous/private/display_info_file.m,
scripts/miscellaneous/python.m, scripts/miscellaneous/rmappdata.m,
scripts/miscellaneous/run.m, scripts/miscellaneous/tar.m,
scripts/miscellaneous/tempname.m, scripts/miscellaneous/untar.m,
scripts/miscellaneous/unzip.m, scripts/miscellaneous/what.m,
scripts/miscellaneous/zip.m, scripts/optimization/fminunc.m,
scripts/optimization/fsolve.m, scripts/optimization/fzero.m,
scripts/optimization/glpk.m, scripts/optimization/optimget.m,
scripts/optimization/optimset.m, scripts/optimization/qp.m,
scripts/optimization/sqp.m, scripts/path/pathdef.m, scripts/pkg/pkg.m,
scripts/pkg/private/build.m, scripts/pkg/private/describe.m,
scripts/pkg/private/dirempty.m, scripts/pkg/private/get_forge_download.m,
scripts/pkg/private/get_forge_pkg.m,
scripts/pkg/private/get_unsatisfied_deps.m, scripts/pkg/private/install.m,
scripts/pkg/private/is_architecture_dependent.m,
scripts/pkg/private/list_forge_packages.m, scripts/pkg/private/rebuild.m,
scripts/pkg/private/shell.m, scripts/pkg/private/uninstall.m,
scripts/plot/axes.m, scripts/plot/box.m, scripts/plot/closereq.m,
scripts/plot/diffuse.m, scripts/plot/ezpolar.m, scripts/plot/findfigs.m,
scripts/plot/gco.m, scripts/plot/guidata.m, scripts/plot/guihandles.m,
scripts/plot/hdl2struct.m, scripts/plot/linkprop.m, scripts/plot/peaks.m,
scripts/plot/print.m, scripts/plot/private/__add_datasource__.m,
scripts/plot/private/__axis_label__.m, scripts/plot/private/__clabel__.m,
scripts/plot/private/__color_str_rgb__.m, scripts/plot/private/__contour__.m,
scripts/plot/private/__default_plot_options__.m,
scripts/plot/private/__errcomm__.m, scripts/plot/private/__file_filter__.m,
scripts/plot/private/__fltk_file_filter__.m,
scripts/plot/private/__getlegenddata__.m,
scripts/plot/private/__gnuplot_open_stream__.m,
scripts/plot/private/__gnuplot_print__.m,
scripts/plot/private/__go_draw_axes__.m,
scripts/plot/private/__interp_cube__.m, scripts/plot/private/__is_function__.m,
scripts/plot/private/__line__.m, scripts/plot/private/__marching_cube__.m,
scripts/plot/private/__next_line_style__.m, scripts/plot/private/__patch__.m,
scripts/plot/private/__pie__.m, scripts/plot/private/__pltopt__.m,
scripts/plot/private/__quiver__.m, scripts/plot/private/__scatter__.m,
scripts/plot/private/__stem__.m, scripts/plot/private/__uigetdir_fltk__.m,
scripts/plot/private/__uigetfile_fltk__.m,
scripts/plot/private/__uiobject_split_args__.m,
scripts/plot/private/__uiputfile_fltk__.m, scripts/plot/refresh.m,
scripts/plot/saveas.m, scripts/plot/shg.m, scripts/plot/specular.m,
scripts/plot/sphere.m, scripts/plot/struct2hdl.m, scripts/plot/subplot.m,
scripts/plot/uicontextmenu.m, scripts/plot/uicontrol.m, scripts/plot/uipanel.m,
scripts/plot/uipushtool.m, scripts/plot/uiresume.m,
scripts/plot/uitoggletool.m, scripts/plot/uitoolbar.m, scripts/plot/uiwait.m,
scripts/plot/waitforbuttonpress.m, scripts/polynomial/pchip.m,
scripts/polynomial/polyeig.m, scripts/polynomial/ppval.m,
scripts/prefs/addpref.m, scripts/prefs/getpref.m, scripts/prefs/ispref.m,
scripts/prefs/private/loadprefs.m, scripts/prefs/private/prefsfile.m,
scripts/prefs/private/saveprefs.m, scripts/prefs/setpref.m,
scripts/set/private/validargs.m, scripts/set/unique.m,
scripts/signal/arch_fit.m, scripts/signal/arch_rnd.m,
scripts/signal/arch_test.m, scripts/signal/arma_rnd.m,
scripts/signal/durbinlevinson.m, scripts/signal/fractdiff.m,
scripts/signal/freqz.m, scripts/signal/freqz_plot.m, scripts/signal/hurst.m,
scripts/signal/periodogram.m, scripts/signal/private/rectangle_lw.m,
scripts/signal/private/rectangle_sw.m, scripts/signal/private/triangle_sw.m,
scripts/signal/spectral_adf.m, scripts/signal/spectral_xdf.m,
scripts/signal/stft.m, scripts/signal/synthesis.m, scripts/signal/yulewalker.m,
scripts/sparse/colperm.m, scripts/sparse/eigs.m, scripts/sparse/etreeplot.m,
scripts/sparse/gmres.m, scripts/sparse/private/__sprand_impl__.m,
scripts/sparse/spdiags.m, scripts/sparse/sprandn.m, scripts/specfun/bessel.m,
scripts/specfun/betaln.m, scripts/specfun/expint.m,
scripts/special-matrix/gallery.m, scripts/startup/__finish__.m,
scripts/statistics/base/qqplot.m, scripts/statistics/distributions/tcdf.m,
scripts/statistics/distributions/wienrnd.m,
scripts/statistics/models/logistic_regression.m,
scripts/statistics/models/private/logistic_regression_derivatives.m,
scripts/statistics/models/private/logistic_regression_likelihood.m,
scripts/statistics/tests/anova.m, scripts/statistics/tests/bartlett_test.m,
scripts/statistics/tests/chisquare_test_homogeneity.m,
scripts/statistics/tests/chisquare_test_independence.m,
scripts/statistics/tests/cor_test.m,
scripts/statistics/tests/f_test_regression.m,
scripts/statistics/tests/hotelling_test.m,
scripts/statistics/tests/hotelling_test_2.m,
scripts/statistics/tests/kolmogorov_smirnov_test_2.m,
scripts/statistics/tests/kruskal_wallis_test.m,
scripts/statistics/tests/manova.m, scripts/statistics/tests/mcnemar_test.m,
scripts/statistics/tests/prop_test_2.m, scripts/statistics/tests/run_test.m,
scripts/statistics/tests/sign_test.m, scripts/statistics/tests/t_test.m,
scripts/statistics/tests/t_test_2.m,
scripts/statistics/tests/t_test_regression.m,
scripts/statistics/tests/u_test.m, scripts/statistics/tests/var_test.m,
scripts/statistics/tests/welch_test.m,
scripts/statistics/tests/wilcoxon_test.m, scripts/statistics/tests/z_test.m,
scripts/statistics/tests/z_test_2.m, scripts/strings/strcat.m,
scripts/strings/strjoin.m, scripts/strings/strsplit.m,
scripts/testfun/__have_feature__.m, scripts/testfun/__printf_assert__.m,
scripts/testfun/__prog_output_assert__.m, scripts/testfun/__run_test_suite__.m,
scripts/time/clock.m, scripts/time/datenum.m, scripts/ui/errordlg.m,
scripts/ui/private/message_dialog.m: End m-files with 1 blank line.
author | Rik <rik@octave.org> |
---|---|
date | Wed, 28 Aug 2013 08:33:02 -0700 |
parents | bc924baa2c4e |
children | d63878346099 |
line wrap: on
line source
## Copyright (C) 2008-2012 VZLU Prague, a.s. ## ## This file is part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 3 of the License, or (at ## your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## <http://www.gnu.org/licenses/>. ## ## Author: Jaroslav Hajek <highegg@gmail.com> ## -*- texinfo -*- ## @deftypefn {Function File} {} fzero (@var{fun}, @var{x0}) ## @deftypefnx {Function File} {} fzero (@var{fun}, @var{x0}, @var{options}) ## @deftypefnx {Function File} {[@var{x}, @var{fval}, @var{info}, @var{output}] =} fzero (@dots{}) ## Find a zero of a univariate function. ## ## @var{fun} is a function handle, inline function, or string ## containing the name of the function to evaluate. ## @var{x0} should be a two-element vector specifying two points which ## bracket a zero. In other words, there must be a change in sign of the ## function between @var{x0}(1) and @var{x0}(2). More mathematically, the ## following must hold ## ## @example ## sign (@var{fun}(@var{x0}(1))) * sign (@var{fun}(@var{x0}(2))) <= 0 ## @end example ## ## If @var{x0} is a single scalar then several nearby and distant ## values are probed in an attempt to obtain a valid bracketing. If this ## is not successful, the function fails. ## @var{options} is a structure specifying additional options. ## Currently, @code{fzero} ## recognizes these options: @qcode{"FunValCheck"}, @qcode{"OutputFcn"}, ## @qcode{"TolX"}, @qcode{"MaxIter"}, @qcode{"MaxFunEvals"}. ## For a description of these options, see @ref{XREFoptimset,,optimset}. ## ## On exit, the function returns @var{x}, the approximate zero point ## and @var{fval}, the function value thereof. ## @var{info} is an exit flag that can have these values: ## ## @itemize ## @item 1 ## The algorithm converged to a solution. ## ## @item 0 ## Maximum number of iterations or function evaluations has been reached. ## ## @item -1 ## The algorithm has been terminated from user output function. ## ## @item -5 ## The algorithm may have converged to a singular point. ## @end itemize ## ## @var{output} is a structure containing runtime information about the ## @code{fzero} algorithm. Fields in the structure are: ## ## @itemize ## @item iterations ## Number of iterations through loop. ## ## @item nfev ## Number of function evaluations. ## ## @item bracketx ## A two-element vector with the final bracketing of the zero along the x-axis. ## ## @item brackety ## A two-element vector with the final bracketing of the zero along the y-axis. ## @end itemize ## @seealso{optimset, fsolve} ## @end deftypefn ## This is essentially the ACM algorithm 748: Enclosing Zeros of ## Continuous Functions due to Alefeld, Potra and Shi, ACM Transactions ## on Mathematical Software, Vol. 21, No. 3, September 1995. Although ## the workflow should be the same, the structure of the algorithm has ## been transformed non-trivially; instead of the authors' approach of ## sequentially calling building blocks subprograms we implement here a ## FSM version using one interior point determination and one bracketing ## per iteration, thus reducing the number of temporary variables and ## simplifying the algorithm structure. Further, this approach reduces ## the need for external functions and error handling. The algorithm has ## also been slightly modified. ## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup. ## PKG_ADD: [~] = __all_opts__ ("fzero"); function [x, fval, info, output] = fzero (fun, x0, options = struct ()) ## Get default options if requested. if (nargin == 1 && ischar (fun) && strcmp (fun, 'defaults')) x = optimset ("MaxIter", Inf, "MaxFunEvals", Inf, "TolX", 1e-8, "OutputFcn", [], "FunValCheck", "off"); return; endif if (nargin < 2 || nargin > 3) print_usage (); endif if (ischar (fun)) fun = str2func (fun, "global"); endif ## TODO ## displev = optimget (options, "Display", "notify"); funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on"); outfcn = optimget (options, "OutputFcn"); tolx = optimget (options, "TolX", 1e-8); maxiter = optimget (options, "MaxIter", Inf); maxfev = optimget (options, "MaxFunEvals", Inf); persistent mu = 0.5; if (funvalchk) ## Replace fun with a guarded version. fun = @(x) guarded_eval (fun, x); endif ## The default exit flag if exceeded number of iterations. info = 0; niter = 0; nfev = 0; x = fval = a = fa = b = fb = NaN; eps = eps (class (x0)); ## Prepare... a = x0(1); fa = fun (a); nfev = 1; if (length (x0) > 1) b = x0(2); fb = fun (b); nfev += 1; else ## Try to get b. if (a == 0) aa = 1; else aa = a; endif 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] fb = fun (b); nfev += 1; if (sign (fa) * sign (fb) <= 0) break; endif endfor endif if (b < a) u = a; a = b; b = u; fu = fa; fa = fb; fb = fu; endif if (! (sign (fa) * sign (fb) <= 0)) error ("fzero:bracket", "fzero: not a valid initial bracketing"); endif slope0 = (fb - fa) / (b - a); if (fa == 0) b = a; fb = fa; elseif (fb == 0) a = b; fa = fb; endif itype = 1; if (abs (fa) < abs (fb)) u = a; fu = fa; else u = b; fu = fb; endif d = e = u; fd = fe = fu; mba = mu*(b - a); while (niter < maxiter && nfev < maxfev) switch (itype) case 1 ## The initial test. if (b - a <= 2*(2 * abs (u) * eps + tolx)) x = u; fval = fu; info = 1; break; endif if (abs (fa) <= 1e3*abs (fb) && abs (fb) <= 1e3*abs (fa)) ## Secant step. c = u - (a - b) / (fa - fb) * fu; else ## Bisection step. c = 0.5*(a + b); endif d = u; fd = fu; itype = 5; case {2, 3} l = length (unique ([fa, fb, fd, fe])); if (l == 4) ## Inverse cubic interpolation. q11 = (d - e) * fd / (fe - fd); q21 = (b - d) * fb / (fd - fb); q31 = (a - b) * fa / (fb - fa); d21 = (b - d) * fd / (fd - fb); d31 = (a - b) * fb / (fb - fa); q22 = (d21 - q11) * fb / (fe - fb); q32 = (d31 - q21) * fa / (fd - fa); d32 = (d31 - q21) * fd / (fd - fa); q33 = (d32 - q22) * fa / (fe - fa); c = a + q31 + q32 + q33; endif if (l < 4 || sign (c - a) * sign (c - b) > 0) ## Quadratic interpolation + newton. a0 = fa; a1 = (fb - fa)/(b - a); a2 = ((fd - fb)/(d - b) - a1) / (d - a); ## Modification 1: this is simpler and does not seem to be worse. c = a - a0/a1; if (a2 != 0) c = a - a0/a1; for i = 1:itype pc = a0 + (a1 + a2*(c - b))*(c - a); pdc = a1 + a2*(2*c - a - b); if (pdc == 0) c = a - a0/a1; break; endif c -= pc/pdc; endfor endif endif itype += 1; case 4 ## Double secant step. c = u - 2*(b - a)/(fb - fa)*fu; ## Bisect if too far. if (abs (c - u) > 0.5*(b - a)) c = 0.5 * (b + a); endif itype = 5; case 5 ## Bisection step. c = 0.5 * (b + a); itype = 2; endswitch ## Don't let c come too close to a or b. delta = 2*0.7*(2 * abs (u) * eps + tolx); if ((b - a) <= 2*delta) c = (a + b)/2; else c = max (a + delta, min (b - delta, c)); endif ## Calculate new point. x = c; fval = fc = fun (c); niter ++; nfev ++; ## Modification 2: skip inverse cubic interpolation if ## nonmonotonicity is detected. if (sign (fc - fa) * sign (fc - fb) >= 0) ## The new point broke monotonicity. ## Disable inverse cubic. fe = fc; else e = d; fe = fd; endif ## Bracketing. if (sign (fa) * sign (fc) < 0) d = b; fd = fb; b = c; fb = fc; elseif (sign (fb) * sign (fc) < 0) d = a; fd = fa; a = c; fa = fc; elseif (fc == 0) a = b = c; fa = fb = fc; info = 1; break; else ## This should never happen. error ("fzero:bracket", "fzero: zero point is not bracketed"); endif ## If there's an output function, use it now. if (outfcn) optv.funccount = nfev; optv.fval = fval; optv.iteration = niter; if (outfcn (x, optv, "iter")) info = -1; break; endif endif if (abs (fa) < abs (fb)) u = a; fu = fa; else u = b; fu = fb; endif if (b - a <= 2*(2 * abs (u) * eps + tolx)) info = 1; break; endif ## Skip bisection step if successful reduction. if (itype == 5 && (b - a) <= mba) itype = 2; endif if (itype == 2) mba = mu * (b - a); endif endwhile ## Check solution for a singularity by examining slope if (info == 1) if ((b - a) != 0 && abs ((fb - fa)/(b - a) / slope0) > max (1e6, 0.5/(eps+tolx))) info = -5; endif endif output.iterations = niter; output.funcCount = nfev; output.bracketx = [a, b]; output.brackety = [fa, fb]; endfunction ## An assistant function that evaluates a function handle and checks for ## bad results. function fx = guarded_eval (fun, x) fx = fun (x); fx = fx(1); if (! isreal (fx)) error ("fzero:notreal", "fzero: non-real value encountered"); elseif (isnan (fx)) error ("fzero:isnan", "fzero: NaN value encountered"); endif endfunction %!shared opt0 %! opt0 = optimset ("tolx", 0); %!assert (fzero (@cos, [0, 3], opt0), pi/2, 10*eps) %!assert (fzero (@(x) x^(1/3) - 1e-8, [0,1], opt0), 1e-24, 1e-22*eps)