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