Mercurial > hg > octave-lyh
annotate scripts/optimization/fzero.m @ 8467:77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Mon, 12 Jan 2009 13:11:35 -0500 |
| parents | 368b504777a8 |
| children | 866492035ecf |
| rev | line source |
|---|---|
| 8305 | 1 ## Copyright (C) 2008 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 | |
|
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 -*- |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
22 ## @deftypefn{Function File}{[@var{x}, @var{fval}, @var{info}, @var{output}] =} fzero (@var{fun}, @var{x0}, @var{options}) |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
23 ## Finds a zero point of a univariate function. @var{fun} should be a function |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
24 ## handle or name. @var{x0} specifies a starting point. @var{options} is a |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
25 ## structure specifying additional options. Currently, fzero recognizes these |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
26 ## options: FunValCheck, OutputFcn, TolX, MaxIter, MaxFunEvals. |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
27 ## For description of these options, see @code{optimset}. |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
28 ## |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
29 ## 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
|
30 ## 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
|
31 ## @var{info} is an exit flag that can have these values: |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
32 ## @itemize |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
33 ## @item 1 |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
34 ## The algorithm converged to a solution. |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
35 ## @item 0 |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
36 ## Maximum number of iterations or function evaluations has been exhausted. |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
37 ## @item -1 |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
38 ## The algorithm has been terminated from user output function. |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
39 ## @item -2 |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
40 ## A general unexpected error. |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
41 ## @item -3 |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
42 ## A non-real value encountered. |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
43 ## @item -4 |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
44 ## A NaN value encountered. |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
45 ## @end itemize |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
46 ## @seealso{optimset, fminbnd, fsolve} |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
47 ## @end deftypefn |
| 8305 | 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 ## 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
|
50 ## 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
|
51 ## 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
|
52 ## 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
|
53 ## 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
|
54 ## 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
|
55 ## 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
|
56 ## 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
|
57 ## 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
|
58 ## 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
|
59 ## 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
|
60 |
| 8305 | 61 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
|
62 |
| 8305 | 63 if (nargin < 2 || nargin > 3) |
| 64 print_usage (); | |
| 65 endif | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
66 |
| 8305 | 67 if (ischar (fun)) |
| 68 fun = str2func (fun); | |
| 69 endif | |
| 70 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
71 ## TODO |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
72 ## 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
|
73 funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on"); |
| 8305 | 74 outfcn = optimget (options, "OutputFcn"); |
| 75 tolx = optimget (options, "TolX", 0); | |
| 76 maxiter = optimget (options, "MaxIter", Inf); | |
| 77 maxfev = optimget (options, "MaxFunEvals", Inf); | |
| 78 | |
| 79 persistent mu = 0.5; | |
| 80 | |
| 81 if (funvalchk) | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
82 ## replace fun with a guarded version |
| 8305 | 83 fun = @(x) guarded_eval (fun, x); |
| 84 endif | |
| 85 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
86 ## the default exit flag if exceeded number of iterations |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
87 info = 0; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
88 niter = 0; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
89 nfev = 0; |
| 8305 | 90 |
| 91 x = fval = a = fa = b = fb = NaN; | |
| 92 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
93 ## prepare... |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
94 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
|
95 fa = fun (a); |
| 8305 | 96 nfev = 1; |
| 97 if (length (x0) > 1) | |
| 98 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
|
99 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
|
100 nfev += 1; |
| 8305 | 101 else |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
102 ## try to get b |
| 8305 | 103 if (a == 0) |
| 104 aa = 1; | |
| 105 else | |
| 106 aa = a; | |
| 107 endif | |
| 108 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] | |
| 109 fb = fun (b); nfev += 1; | |
| 110 if (sign (fa) * sign (fb) <= 0) | |
| 111 break; | |
| 112 endif | |
| 113 endfor | |
| 114 endif | |
| 115 | |
| 116 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
|
117 u = a; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
118 a = b; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
119 b = u; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
120 |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
121 fu = fa; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
122 fa = fb; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
123 fb = fu; |
| 8305 | 124 endif |
| 125 | |
| 126 if (! (sign (fa) * sign (fb) <= 0)) | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
127 error ("fzero: not a valid initial bracketing"); |
| 8305 | 128 endif |
| 129 | |
| 130 itype = 1; | |
| 131 | |
| 132 if (abs (fa) < abs (fb)) | |
| 133 u = a; fu = fa; | |
| 134 else | |
| 135 u = b; fu = fb; | |
| 136 endif | |
| 137 | |
| 138 d = e = u; | |
| 139 fd = fe = fu; | |
| 140 mba = mu*(b - a); | |
| 141 while (niter < maxiter && nfev < maxfev) | |
| 142 switch (itype) | |
| 143 case 1 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
144 ## the initial test |
| 8305 | 145 if (b - a <= 2*(2 * abs (u) * eps + tolx)) |
| 146 x = u; fval = fu; | |
| 147 info = 1; | |
| 148 break; | |
| 149 endif | |
| 150 if (abs (fa) <= 1e3*abs (fb) && abs (fb) <= 1e3*abs (fa)) | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
151 ## secant step |
| 8305 | 152 c = u - (a - b) / (fa - fb) * fu; |
| 153 else | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
154 ## bisection step |
| 8305 | 155 c = 0.5*(a + b); |
| 156 endif | |
| 157 d = u; fd = fu; | |
| 158 itype = 5; | |
| 159 case {2, 3} | |
| 160 l = length (unique ([fa, fb, fd, fe])); | |
| 161 if (l == 4) | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
162 ## inverse cubic interpolation |
| 8305 | 163 q11 = (d - e) * fd / (fe - fd); |
| 164 q21 = (b - d) * fb / (fd - fb); | |
| 165 q31 = (a - b) * fa / (fb - fa); | |
| 166 d21 = (b - d) * fd / (fd - fb); | |
| 167 d31 = (a - b) * fb / (fb - fa); | |
| 168 q22 = (d21 - q11) * fb / (fe - fb); | |
| 169 q32 = (d31 - q21) * fa / (fd - fa); | |
| 170 d32 = (d31 - q21) * fd / (fd - fa); | |
| 171 q33 = (d32 - q22) * fa / (fe - fa); | |
| 172 c = a + q31 + q32 + q33; | |
| 173 endif | |
| 174 if (l < 4 || sign (c - a) * sign (c - b) > 0) | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
175 ## quadratic interpolation + newton |
| 8305 | 176 a0 = fa; |
| 177 a1 = (fb - fa)/(b - a); | |
| 178 a2 = ((fd - fb)/(d - b) - a1) / (d - a); | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
179 ## modification 1: this is simpler and does not seem to be worse |
| 8305 | 180 c = a - a0/a1; |
| 181 if (a2 != 0) | |
| 182 c = a - a0/a1; | |
| 183 for i = 1:itype | |
| 184 pc = a0 + (a1 + a2*(c - b))*(c - a); | |
| 185 pdc = a1 + a2*(2*c - a - b); | |
| 186 if (pdc == 0) | |
| 187 c = a - a0/a1; | |
| 188 break; | |
| 189 endif | |
| 190 c -= pc/pdc; | |
| 191 endfor | |
| 192 endif | |
| 193 endif | |
| 194 itype += 1; | |
| 195 case 4 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
196 ## double secant step |
| 8305 | 197 c = u - 2*(b - a)/(fb - fa)*fu; |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
198 ## bisect if too far |
| 8305 | 199 if (abs (c - u) > 0.5*(b - a)) |
| 200 c = 0.5 * (b + a); | |
| 201 endif | |
| 202 itype = 5; | |
| 203 case 5 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
204 ## bisection step |
| 8305 | 205 c = 0.5 * (b + a); |
| 206 itype = 2; | |
| 207 endswitch | |
| 208 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
209 ## don't let c come too close to a or b |
| 8305 | 210 delta = 2*0.7*(2 * abs (u) * eps + tolx); |
| 211 if ((b - a) <= 2*delta) | |
| 212 c = (a + b)/2; | |
| 213 else | |
| 214 c = max (a + delta, min (b - delta, c)); | |
| 215 endif | |
| 216 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
217 ## calculate new point |
| 8305 | 218 x = c; |
| 219 fval = fc = fun (c); | |
| 220 niter ++; nfev ++; | |
| 221 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
222 ## modification 2: skip inverse cubic interpolation if |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
223 ## nonmonotonicity is detected |
| 8305 | 224 if (sign (fc - fa) * sign (fc - fb) >= 0) |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
225 ## the new point broke monotonicity. |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
226 ## disable inverse cubic |
| 8305 | 227 fe = fc; |
| 228 else | |
| 229 e = d; fe = fd; | |
| 230 endif | |
| 231 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
232 ## bracketing |
| 8305 | 233 if (sign (fa) * sign (fc) < 0) |
| 234 d = b; fd = fb; | |
| 235 b = c; fb = fc; | |
| 236 elseif (sign (fb) * sign (fc) < 0) | |
| 237 d = a; fd = fa; | |
| 238 a = c; fa = fc; | |
| 239 elseif (fc == 0) | |
| 240 a = b = c; fa = fb = fc; | |
| 241 info = 1; | |
| 242 break; | |
| 243 else | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
244 ## this should never happen. |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
245 error ("fzero: zero point is not bracketed"); |
| 8305 | 246 endif |
| 247 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
248 ## if there's an output function, use it now |
| 8305 | 249 if (outfcn) |
| 250 optv.funccount = niter + 2; | |
| 251 optv.fval = fval; | |
| 252 optv.iteration = niter; | |
| 253 if (outfcn (x, optv, "iter")) | |
| 254 info = -1; | |
| 255 break; | |
| 256 endif | |
| 257 endif | |
| 258 | |
| 259 if (abs (fa) < abs (fb)) | |
| 260 u = a; fu = fa; | |
| 261 else | |
| 262 u = b; fu = fb; | |
| 263 endif | |
| 264 if (b - a <= 2*(2 * abs (u) * eps + tolx)) | |
| 265 info = 1; | |
| 266 break; | |
| 267 endif | |
| 268 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
269 ## skip bisection step if successful reduction |
| 8305 | 270 if (itype == 5 && (b - a) <= mba) |
| 271 itype = 2; | |
| 272 endif | |
| 273 if (itype == 2) | |
| 274 mba = mu * (b - a); | |
| 275 endif | |
| 276 endwhile | |
| 277 | |
| 278 output.iterations = niter; | |
| 279 output.funcCount = niter + 2; | |
| 280 output.bracket = [a, b]; | |
| 281 output.bracketf = [fa, fb]; | |
| 282 | |
| 283 endfunction | |
| 284 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
285 ## an assistant function that evaluates a function handle and checks for |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
286 ## bad results. |
| 8305 | 287 function fx = guarded_eval (fun, x) |
| 288 fx = fun (x); | |
| 289 fx = fx(1); | |
| 290 if (! isreal (fx)) | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
291 error ("fzero: non-real value encountered"); |
| 8305 | 292 elseif (isnan (fx)) |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8305
diff
changeset
|
293 error ("fzero: NaN value encountered"); |
| 8305 | 294 endif |
| 295 endfunction | |
| 296 | |
| 297 %!assert(fzero(@cos, [0, 3]), pi/2, 10*eps) | |
| 298 %!assert(fzero(@(x) x^(1/3) - 1e-8, [0,1]), 1e-24, 1e-22*eps) |