Mercurial > hg > octave-lyh
annotate scripts/optimization/fzero.m @ 11346:30f54b3b9953
drop support of 'i' option in lookup
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Sun, 12 Dec 2010 21:58:30 +0100 |
| parents | 693e22af08ae |
| children | fd0a3ac60b0e |
| rev | line source |
|---|---|
| 8920 | 1 ## Copyright (C) 2008, 2009 VZLU Prague, a.s. |
| 8305 | 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 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
21 ## -*- texinfo -*- |
|
10687
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
22 ## @deftypefn {Function File} {} fzero (@var{fun}, @var{x0}) |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
23 ## @deftypefnx {Function File} {} fzero (@var{fun}, @var{x0}, @var{options}) |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
24 ## @deftypefnx {Function File} {[@var{x}, @var{fval}, @var{info}, @var{output}] =} fzero (@dots{}) |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
25 ## Find a zero of a univariate function. |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
26 ## |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
27 ## @var{fun} should be a function handle or name. @var{x0} should be a |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
28 ## two-element vector specifying two points which bracket a zero. In |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
29 ## other words, there must be a change in sign of the function between |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
30 ## @var{x0}(1) and @var{x0}(2). More mathematically, the following must hold |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
31 ## |
| 10470 | 32 ## @example |
| 33 ## sign (@var{fun}(@var{x0}(1))) * sign (@var{fun}(@var{x0}(2))) <= 0 | |
| 34 ## @end example | |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
35 ## |
|
10687
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
36 ## If @var{x0} is a single scalar then several nearby and distant |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
37 ## values are probed in an attempt to obtain a valid bracketing. If this |
| 10470 | 38 ## is not successful, the function fails. |
| 39 ## @var{options} is a structure specifying additional options. | |
| 40 ## Currently, @code{fzero} | |
|
8515
ec2715c76039
fzero.m, fsolve.m: additional doc fixes
John W. Eaton <jwe@octave.org>
parents:
8507
diff
changeset
|
41 ## recognizes these options: @code{"FunValCheck"}, @code{"OutputFcn"}, |
|
ec2715c76039
fzero.m, fsolve.m: additional doc fixes
John W. Eaton <jwe@octave.org>
parents:
8507
diff
changeset
|
42 ## @code{"TolX"}, @code{"MaxIter"}, @code{"MaxFunEvals"}. |
|
10687
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
43 ## For a description of these options, see @ref{doc-optimset,,optimset}. |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
44 ## |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
45 ## On exit, the function returns @var{x}, the approximate zero point |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
46 ## and @var{fval}, the function value thereof. |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
47 ## @var{info} is an exit flag that can have these values: |
|
10297
ed88ea036716
improve docs of fzero/fminbnd
Jaroslav Hajek <highegg@gmail.com>
parents:
10296
diff
changeset
|
48 ## |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
49 ## @itemize |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
50 ## @item 1 |
|
10687
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
51 ## The algorithm converged to a solution. |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
52 ## |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
53 ## @item 0 |
|
10687
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
54 ## Maximum number of iterations or function evaluations has been reached. |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
55 ## |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
56 ## @item -1 |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
57 ## The algorithm has been terminated from user output function. |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
58 ## |
|
10337
7c97da90fc8f
check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents:
10297
diff
changeset
|
59 ## @item -5 |
|
7c97da90fc8f
check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents:
10297
diff
changeset
|
60 ## The algorithm may have converged to a singular point. |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
61 ## @end itemize |
|
10687
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
62 ## |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
63 ## @var{output} is a structure containing runtime information about the |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
64 ## @code{fzero} algorithm. Fields in the structure are: |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
65 ## |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
66 ## @itemize |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
67 ## @item iterations |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
68 ## Number of iterations through loop. |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
69 ## |
|
10687
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
70 ## @item nfev |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
71 ## Number of function evaluations. |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
72 ## |
|
10687
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
73 ## @item bracketx |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
74 ## A two-element vector with the final bracketing of the zero along the x-axis. |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
75 ## |
|
10687
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
76 ## @item brackety |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
77 ## A two-element vector with the final bracketing of the zero along the y-axis. |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
78 ## @end itemize |
|
8515
ec2715c76039
fzero.m, fsolve.m: additional doc fixes
John W. Eaton <jwe@octave.org>
parents:
8507
diff
changeset
|
79 ## @seealso{optimset, fsolve} |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
80 ## @end deftypefn |
| 8305 | 81 |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
82 ## This is essentially the ACM algorithm 748: Enclosing Zeros of |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
83 ## Continuous Functions due to Alefeld, Potra and Shi, ACM Transactions |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
84 ## on Mathematical Software, Vol. 21, No. 3, September 1995. Although |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
85 ## the workflow should be the same, the structure of the algorithm has |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
86 ## been transformed non-trivially; instead of the authors' approach of |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
87 ## sequentially calling building blocks subprograms we implement here a |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
88 ## FSM version using one interior point determination and one bracketing |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
89 ## per iteration, thus reducing the number of temporary variables and |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
90 ## simplifying the algorithm structure. Further, this approach reduces |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
91 ## the need for external functions and error handling. The algorithm has |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
92 ## also been slightly modified. |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
93 |
|
8648
ff61b53eb294
optimization: use PKG_ADD: comments instead of PKG_ADD file
John W. Eaton <jwe@octave.org>
parents:
8596
diff
changeset
|
94 ## PKG_ADD: __all_opts__ ("fzero"); |
|
ff61b53eb294
optimization: use PKG_ADD: comments instead of PKG_ADD file
John W. Eaton <jwe@octave.org>
parents:
8596
diff
changeset
|
95 |
| 8305 | 96 function [x, fval, info, output] = fzero (fun, x0, options = struct ()) |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
97 |
|
8596
8833c0b18eb2
enable default settings queries in optim funcs
Jaroslav Hajek <highegg@gmail.com>
parents:
8515
diff
changeset
|
98 ## Get default options if requested. |
|
8833c0b18eb2
enable default settings queries in optim funcs
Jaroslav Hajek <highegg@gmail.com>
parents:
8515
diff
changeset
|
99 if (nargin == 1 && ischar (fun) && strcmp (fun, 'defaults')) |
| 10296 | 100 x = optimset ("MaxIter", Inf, "MaxFunEvals", Inf, "TolX", 1e-8, \ |
|
8596
8833c0b18eb2
enable default settings queries in optim funcs
Jaroslav Hajek <highegg@gmail.com>
parents:
8515
diff
changeset
|
101 "OutputFcn", [], "FunValCheck", "off"); |
|
8833c0b18eb2
enable default settings queries in optim funcs
Jaroslav Hajek <highegg@gmail.com>
parents:
8515
diff
changeset
|
102 return; |
|
8833c0b18eb2
enable default settings queries in optim funcs
Jaroslav Hajek <highegg@gmail.com>
parents:
8515
diff
changeset
|
103 endif |
|
8833c0b18eb2
enable default settings queries in optim funcs
Jaroslav Hajek <highegg@gmail.com>
parents:
8515
diff
changeset
|
104 |
| 8305 | 105 if (nargin < 2 || nargin > 3) |
| 106 print_usage (); | |
| 107 endif | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
108 |
| 8305 | 109 if (ischar (fun)) |
|
9464
e598248a060d
safer str2func use in optim functions
Jaroslav Hajek <highegg@gmail.com>
parents:
9143
diff
changeset
|
110 fun = str2func (fun, "global"); |
| 8305 | 111 endif |
| 112 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
113 ## TODO |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
114 ## displev = optimget (options, "Display", "notify"); |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
115 funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on"); |
| 8305 | 116 outfcn = optimget (options, "OutputFcn"); |
| 10296 | 117 tolx = optimget (options, "TolX", 1e-8); |
| 8305 | 118 maxiter = optimget (options, "MaxIter", Inf); |
| 119 maxfev = optimget (options, "MaxFunEvals", Inf); | |
| 120 | |
| 121 persistent mu = 0.5; | |
| 122 | |
| 123 if (funvalchk) | |
| 8507 | 124 ## Replace fun with a guarded version. |
| 8305 | 125 fun = @(x) guarded_eval (fun, x); |
| 126 endif | |
| 127 | |
| 8507 | 128 ## The default exit flag if exceeded number of iterations. |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
129 info = 0; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
130 niter = 0; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
131 nfev = 0; |
| 8305 | 132 |
| 133 x = fval = a = fa = b = fb = NaN; | |
| 10296 | 134 eps = eps (class (x0)); |
| 8305 | 135 |
| 8507 | 136 ## Prepare... |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
137 a = x0(1); |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
138 fa = fun (a); |
| 8305 | 139 nfev = 1; |
| 140 if (length (x0) > 1) | |
| 141 b = x0(2); | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
142 fb = fun (b); |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
143 nfev += 1; |
| 8305 | 144 else |
| 8507 | 145 ## Try to get b. |
| 8305 | 146 if (a == 0) |
| 147 aa = 1; | |
| 148 else | |
| 149 aa = a; | |
| 150 endif | |
| 151 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] | |
| 152 fb = fun (b); nfev += 1; | |
| 153 if (sign (fa) * sign (fb) <= 0) | |
| 154 break; | |
| 155 endif | |
| 156 endfor | |
| 157 endif | |
| 158 | |
| 159 if (b < a) | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
160 u = a; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
161 a = b; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
162 b = u; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
163 |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
164 fu = fa; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
165 fa = fb; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
166 fb = fu; |
| 8305 | 167 endif |
| 168 | |
| 169 if (! (sign (fa) * sign (fb) <= 0)) | |
|
8468
866492035ecf
fsolve.m, fzero.m: undo part of previous change
John W. Eaton <jwe@octave.org>
parents:
8467
diff
changeset
|
170 error ("fzero:bracket", "fzero: not a valid initial bracketing"); |
| 8305 | 171 endif |
| 172 | |
|
10337
7c97da90fc8f
check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents:
10297
diff
changeset
|
173 slope0 = (fb - fa) / (b - a); |
|
7c97da90fc8f
check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents:
10297
diff
changeset
|
174 |
| 10470 | 175 if (fa == 0) |
| 176 b = a; | |
| 177 fb = fa; | |
| 178 elseif (fb == 0) | |
| 179 a = b; | |
| 180 fa = fb; | |
| 181 endif | |
| 182 | |
| 8305 | 183 itype = 1; |
| 184 | |
| 185 if (abs (fa) < abs (fb)) | |
| 186 u = a; fu = fa; | |
| 187 else | |
| 188 u = b; fu = fb; | |
| 189 endif | |
| 190 | |
| 191 d = e = u; | |
| 192 fd = fe = fu; | |
| 193 mba = mu*(b - a); | |
| 194 while (niter < maxiter && nfev < maxfev) | |
| 195 switch (itype) | |
| 196 case 1 | |
| 8507 | 197 ## The initial test. |
| 8305 | 198 if (b - a <= 2*(2 * abs (u) * eps + tolx)) |
| 199 x = u; fval = fu; | |
| 200 info = 1; | |
| 201 break; | |
| 202 endif | |
| 203 if (abs (fa) <= 1e3*abs (fb) && abs (fb) <= 1e3*abs (fa)) | |
| 8507 | 204 ## Secant step. |
| 8305 | 205 c = u - (a - b) / (fa - fb) * fu; |
| 206 else | |
| 8507 | 207 ## Bisection step. |
| 8305 | 208 c = 0.5*(a + b); |
| 209 endif | |
| 210 d = u; fd = fu; | |
| 211 itype = 5; | |
| 212 case {2, 3} | |
| 213 l = length (unique ([fa, fb, fd, fe])); | |
| 214 if (l == 4) | |
| 8507 | 215 ## Inverse cubic interpolation. |
| 8305 | 216 q11 = (d - e) * fd / (fe - fd); |
| 217 q21 = (b - d) * fb / (fd - fb); | |
| 218 q31 = (a - b) * fa / (fb - fa); | |
| 219 d21 = (b - d) * fd / (fd - fb); | |
| 220 d31 = (a - b) * fb / (fb - fa); | |
| 221 q22 = (d21 - q11) * fb / (fe - fb); | |
| 222 q32 = (d31 - q21) * fa / (fd - fa); | |
| 223 d32 = (d31 - q21) * fd / (fd - fa); | |
| 224 q33 = (d32 - q22) * fa / (fe - fa); | |
| 225 c = a + q31 + q32 + q33; | |
| 226 endif | |
| 227 if (l < 4 || sign (c - a) * sign (c - b) > 0) | |
| 8507 | 228 ## Quadratic interpolation + newton. |
| 8305 | 229 a0 = fa; |
| 230 a1 = (fb - fa)/(b - a); | |
| 231 a2 = ((fd - fb)/(d - b) - a1) / (d - a); | |
| 8507 | 232 ## Modification 1: this is simpler and does not seem to be worse. |
| 8305 | 233 c = a - a0/a1; |
| 234 if (a2 != 0) | |
| 235 c = a - a0/a1; | |
| 236 for i = 1:itype | |
| 237 pc = a0 + (a1 + a2*(c - b))*(c - a); | |
| 238 pdc = a1 + a2*(2*c - a - b); | |
| 239 if (pdc == 0) | |
| 240 c = a - a0/a1; | |
| 241 break; | |
| 242 endif | |
| 243 c -= pc/pdc; | |
| 244 endfor | |
| 245 endif | |
| 246 endif | |
| 247 itype += 1; | |
| 248 case 4 | |
| 8507 | 249 ## Double secant step. |
| 8305 | 250 c = u - 2*(b - a)/(fb - fa)*fu; |
| 8507 | 251 ## Bisect if too far. |
| 8305 | 252 if (abs (c - u) > 0.5*(b - a)) |
| 253 c = 0.5 * (b + a); | |
| 254 endif | |
| 255 itype = 5; | |
| 256 case 5 | |
| 8507 | 257 ## Bisection step. |
| 8305 | 258 c = 0.5 * (b + a); |
| 259 itype = 2; | |
| 260 endswitch | |
| 261 | |
| 8507 | 262 ## Don't let c come too close to a or b. |
| 8305 | 263 delta = 2*0.7*(2 * abs (u) * eps + tolx); |
| 264 if ((b - a) <= 2*delta) | |
| 265 c = (a + b)/2; | |
| 266 else | |
| 267 c = max (a + delta, min (b - delta, c)); | |
| 268 endif | |
| 269 | |
| 8507 | 270 ## Calculate new point. |
| 8305 | 271 x = c; |
| 272 fval = fc = fun (c); | |
| 273 niter ++; nfev ++; | |
| 274 | |
| 8507 | 275 ## Modification 2: skip inverse cubic interpolation if |
| 276 ## nonmonotonicity is detected. | |
| 8305 | 277 if (sign (fc - fa) * sign (fc - fb) >= 0) |
| 8507 | 278 ## The new point broke monotonicity. |
| 279 ## Disable inverse cubic. | |
| 8305 | 280 fe = fc; |
| 281 else | |
| 282 e = d; fe = fd; | |
| 283 endif | |
| 284 | |
| 8507 | 285 ## Bracketing. |
| 8305 | 286 if (sign (fa) * sign (fc) < 0) |
| 287 d = b; fd = fb; | |
| 288 b = c; fb = fc; | |
| 289 elseif (sign (fb) * sign (fc) < 0) | |
| 290 d = a; fd = fa; | |
| 291 a = c; fa = fc; | |
| 292 elseif (fc == 0) | |
| 293 a = b = c; fa = fb = fc; | |
| 294 info = 1; | |
| 295 break; | |
| 296 else | |
| 8507 | 297 ## This should never happen. |
|
8468
866492035ecf
fsolve.m, fzero.m: undo part of previous change
John W. Eaton <jwe@octave.org>
parents:
8467
diff
changeset
|
298 error ("fzero:bracket", "fzero: zero point is not bracketed"); |
| 8305 | 299 endif |
| 300 | |
| 8507 | 301 ## If there's an output function, use it now. |
| 8305 | 302 if (outfcn) |
| 10296 | 303 optv.funccount = nfev; |
| 8305 | 304 optv.fval = fval; |
| 305 optv.iteration = niter; | |
| 306 if (outfcn (x, optv, "iter")) | |
| 307 info = -1; | |
| 308 break; | |
| 309 endif | |
| 310 endif | |
| 311 | |
| 312 if (abs (fa) < abs (fb)) | |
| 313 u = a; fu = fa; | |
| 314 else | |
| 315 u = b; fu = fb; | |
| 316 endif | |
| 317 if (b - a <= 2*(2 * abs (u) * eps + tolx)) | |
| 318 info = 1; | |
| 319 break; | |
| 320 endif | |
| 321 | |
| 8507 | 322 ## Skip bisection step if successful reduction. |
| 8305 | 323 if (itype == 5 && (b - a) <= mba) |
| 324 itype = 2; | |
| 325 endif | |
| 326 if (itype == 2) | |
| 327 mba = mu * (b - a); | |
| 328 endif | |
| 329 endwhile | |
| 330 | |
|
10337
7c97da90fc8f
check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents:
10297
diff
changeset
|
331 ## Check solution for a singularity by examining slope |
|
7c97da90fc8f
check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents:
10297
diff
changeset
|
332 if (info == 1) |
|
10394
bc475cd49147
avoid div by zero in fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
10337
diff
changeset
|
333 if ((b - a) != 0 && abs ((fb - fa)/(b - a) / slope0) > max (1e6, 0.5/(eps+tolx))) |
|
10337
7c97da90fc8f
check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents:
10297
diff
changeset
|
334 info = -5; |
|
7c97da90fc8f
check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents:
10297
diff
changeset
|
335 endif |
|
7c97da90fc8f
check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents:
10297
diff
changeset
|
336 endif |
|
7c97da90fc8f
check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents:
10297
diff
changeset
|
337 |
| 8305 | 338 output.iterations = niter; |
| 10296 | 339 output.funcCount = nfev; |
|
10687
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
340 output.bracketx = [a, b]; |
|
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
10470
diff
changeset
|
341 output.brackety = [fa, fb]; |
| 8305 | 342 |
| 343 endfunction | |
| 344 | |
| 8507 | 345 ## An assistant function that evaluates a function handle and checks for |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
346 ## bad results. |
| 8305 | 347 function fx = guarded_eval (fun, x) |
| 348 fx = fun (x); | |
| 349 fx = fx(1); | |
| 350 if (! isreal (fx)) | |
|
8468
866492035ecf
fsolve.m, fzero.m: undo part of previous change
John W. Eaton <jwe@octave.org>
parents:
8467
diff
changeset
|
351 error ("fzero:notreal", "fzero: non-real value encountered"); |
| 8305 | 352 elseif (isnan (fx)) |
|
8468
866492035ecf
fsolve.m, fzero.m: undo part of previous change
John W. Eaton <jwe@octave.org>
parents:
8467
diff
changeset
|
353 error ("fzero:isnan", "fzero: NaN value encountered"); |
| 8305 | 354 endif |
| 355 endfunction | |
| 356 | |
| 10296 | 357 %!shared opt0 |
| 358 %! opt0 = optimset ("tolx", 0); | |
| 359 %!assert(fzero(@cos, [0, 3], opt0), pi/2, 10*eps) | |
| 360 %!assert(fzero(@(x) x^(1/3) - 1e-8, [0,1], opt0), 1e-24, 1e-22*eps) |