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