Mercurial > hg > octave-nkf

comparison scripts/optimization/fzero.m @ 8305:368b504777a8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
implement fzero
author Jaroslav Hajek <highegg@gmail.com>
date Fri, 31 Oct 2008 08:06:45 +0100
parents
children 77b8d4aa2743
comparison
equal deleted inserted replaced
8304:eeaee297c0da 8305:368b504777a8
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)