Mercurial > hg > octave-nkf

annotate scripts/optimization/fsolve.m @ 8306:43795cf108d0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
initial implementation of fsolve remove old fsolve code
author Jaroslav Hajek <highegg@gmail.com>
date Sun, 28 Sep 2008 21:09:35 +0200
parents
children a762d9daa700
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
1 ## Copyright (C) 2008 VZLU Prague, a.s.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
2 ##
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
3 ## This file is part of Octave.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
4 ##
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
5 ## Octave is free software; you can redistribute it and/or modify it
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
6 ## under the terms of the GNU General Public License as published by
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
7 ## the Free Software Foundation; either version 3 of the License, or (at
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
8 ## your option) any later version.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
9 ##
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
10 ## Octave is distributed in the hope that it will be useful, but
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
13 ## General Public License for more details.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
14 ##
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
15 ## You should have received a copy of the GNU General Public License
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
16 ## along with Octave; see the file COPYING. If not, see
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
17 ## <http://www.gnu.org/licenses/>.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
18 ##
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
19 ## Author: Jaroslav Hajek <highegg@gmail.com>
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
20
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
21 # -*- texinfo -*-
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
22 # @deftypefn{Function File} {} fsolve(@var{fcn}, @var{x0}, @var{options})
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
23 # @deftypefnx{Function File} {[@var{x}, @var{fvec}, @var{info}, @var{output}, @var{fjac}]} = fsolve (@var{fcn}, @dots{})
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
24 # Solves a system of nonlinear equations defined by the function @var{fcn}.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
25 # @var{fcn} should accepts a vector (array) defining the unknown variables,
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
26 # and return a vector of left-hand sides of the equations. Right-hand sides
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
27 # are defined to be zeros.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
28 # In other words, this function attempts to determine a vector @var{X} such
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
29 # that @code{@var{fcn}(@var{X})} gives (approximately) all zeros.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
30 # @var{x0} determines a starting guess. The shape of @var{x0} is preserved
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
31 # in all calls to @var{fcn}, but otherwise it is treated as a column vector.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
32 # @var{options} is a structure specifying additional options. Currently, fsolve
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
33 # recognizes these options: FunValCheck, OutputFcn, TolX, TolFun, MaxIter,
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
34 # MaxFunEvals and Jacobian.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
35 #
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
36 # If Jacobian is 'on', it specifies that @var{fcn}, called with 2 output arguments,
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
37 # also returns the Jacobian matrix of right-hand sides at the requested point.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
38 # TolX specifies the termination tolerance in the unknown variables, while
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
39 # TolFun is a tolerance for equations. Default is @code{1e1*eps}
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
40 # for TolX and @code{1e2*eps} for TolFun.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
41 # For description of the other options, see @code{optimset}.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
42 #
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
43 # On return, @var{fval} contains the value of the function @var{fcn}
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
44 # evaluated at @var{x}, and @var{info} may be one of the following values:
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
45 #
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
46 # @table @asis
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
47 # @item 1
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
48 # Converged to a solution point. Relative residual error is less than specified
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
49 # by TolFun.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
50 # @item 2
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
51 # Last relative step size was less that TolX.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
52 # @item 3
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
53 # Last relative decrease in residual was less than TolF.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
54 # @item 0
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
55 # Iteration limit exceeded.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
56 # @item -3
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
57 # The trust region radius became excessively small.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
58 # @end table
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
59 #
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
60 # Note: If you only have a single nonlinear equation of one variable, using
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
61 # @code{fzero} is usually a much better idea.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
62 # @seealso{fzero,optimset}
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
63 # @end deftypefn
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
64
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
65 function [x, fvec, info, output, fjac] = fsolve (fcn, x0, options)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
66
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
67 if (nargin < 3)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
68 options = struct ();
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
69 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
70
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
71 xsiz = size (x0);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
72 n = numel (x0);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
73
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
74 has_jac = strcmp (optimget (options, "Jacobian", "off"), "on");
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
75 maxiter = optimget (options, "MaxIter", Inf);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
76 maxfev = optimget (options, "MaxFunEvals", Inf);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
77 outfcn = optimget (options, "OutputFcn");
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
78 funvalchk = strcmp (optimget (options, "FunValCheck", "off"), "on");
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
79
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
80 if (funvalchk)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
81 # replace fun with a guarded version
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
82 fun = @(x) guarded_eval (fun, x);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
83 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
84
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
85 # These defaults are rather stringent. I think that normally, user prefers
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
86 # accuracy to performance.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
87
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
88 macheps = eps (class (x0));
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
89
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
90 tolx = optimget (options, "TolX", 1e1*macheps);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
91 tolf = optimget (options, "TolFun",1e2*macheps);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
92
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
93 factor = 100;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
94 # FIXME: TypicalX corresponds to user scaling (???)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
95 autodg = true;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
96
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
97 niter = 1; nfev = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
98
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
99 x = x0(:);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
100 info = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
101
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
102 # outer loop
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
103 while (niter < maxiter && nfev < maxfev && ! info)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
104
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
105 # calc func value and jacobian (possibly via FD)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
106 # handle arbitrary shapes of x and f and remember them
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
107 if (has_jac)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
108 [fvec, fjac] = fcn (reshape (x, xsiz));
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
109 nfev ++;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
110 else
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
111 [fvec, fjac] = __fdjac__ (fcn, reshape (x, xsiz));
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
112 nfev += 1 + length (x);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
113 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
114 fsiz = size (fvec);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
115 fvec = fvec(:);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
116 fn = norm (fvec);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
117
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
118 # get QR factorization of the jacobian
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
119 [q, r] = qr (fjac);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
120
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
121 # get column norms, use them as scaling factor
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
122 jcn = norm (fjac, 'columns').';
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
123 if (niter == 1)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
124 if (autodg)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
125 dg = jcn;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
126 dg(dg == 0) = 1;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
127 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
128 xn = norm (dg .* x);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
129 delta = factor * xn;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
130 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
131
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
132 # rescale if necessary
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
133 if (autodg)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
134 dg = max (dg, jcn);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
135 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
136
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
137 nfail = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
138 nsuc = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
139 # inner loop
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
140 while (niter <= maxiter && nfev < maxfev && ! info)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
141
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
142 qtf = q'*fvec;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
143
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
144 # get TR model (dogleg) minimizer
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
145 p = - __dogleg__ (r, qtf, dg, delta);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
146 pn = norm (dg .* p);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
147
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
148 if (niter == 1)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
149 delta = min (delta, pn);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
150 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
151
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
152 fvec1 = fcn (reshape (x + p, xsiz)) (:);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
153 fn1 = norm (fvec1);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
154
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
155 if (fn1 < fn)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
156 # scaled actual reduction
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
157 actred = 1 - (fn1/fn)^2;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
158 else
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
159 actred = -1;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
160 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
161
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
162 # scaled predicted reduction, and ratio
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
163 w = qtf + r*p;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
164 t = norm (w);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
165 if (t < fn)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
166 prered = 1 - (t/fn)^2;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
167 ratio = actred / prered;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
168 else
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
169 prered = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
170 ratio = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
171 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
172
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
173 # update delta
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
174 if (ratio < 0.1)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
175 nsuc = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
176 nfail ++;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
177 delta *= 0.5;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
178 if (delta <= sqrt (macheps)*xn)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
179 # trust region became uselessly small
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
180 info = -3;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
181 break;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
182 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
183 else
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
184 nfail = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
185 nsuc ++;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
186 if (abs (1-ratio) <= 0.1)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
187 delta = 2*pn;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
188 elseif (ratio >= 0.5 || nsuc > 1)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
189 delta = max (delta, 2*pn);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
190 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
191 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
192
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
193 if (ratio >= 1e-4)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
194 # successful iteration
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
195 x += p;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
196 xn = norm (dg .* x);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
197 fvec = fvec1;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
198 fn = fn1;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
199 niter ++;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
200 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
201
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
202 # Tests for termination conditions. A mysterious place, anything can
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
203 # happen if you change something here...
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
204
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
205 # The rule of thumb (which I'm not sure M*b is quite following) is that
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
206 # for a tolerance that depends on scaling, only 0 makes sense as a
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
207 # default value. But 0 usually means uselessly long iterations,
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
208 # so we need scaling-independent tolerances wherever possible.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
209
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
210 # XXX: why tolf*n*xn? If abs (e) ~ abs(x) * eps is a vector of
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
211 # perturbations of x, then norm (fjac*e) <= eps*n*xn, i.e. by tolf ~
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
212 # eps we demand as much accuracy as we can expect.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
213 if (fn <= tolf*n*xn)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
214 info = 1;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
215 # The following tests done only after successful step.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
216 elseif (actred > 0)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
217 # This one is classic. Note that we use scaled variables again, but
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
218 # compare to scaled step, so nothing bad.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
219 if (pn <= tolx*xn)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
220 info = 2;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
221 # Again a classic one. It seems weird to use the same tolf for two
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
222 # different tests, but that's what M*b manual appears to say.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
223 elseif (actred < tolf)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
224 info = 3
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
225 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
226 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
227
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
228 # criterion for recalculating jacobian
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
229 if (nfail == 2)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
230 break;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
231 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
232
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
233 # compute the scaled Broyden update
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
234 u = (fvec1 - q*w) / pn;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
235 v = dg .* ((dg .* p) / pn);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
236
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
237 # update the QR factorization
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
238 [q, r] = qrupdate (q, r, u, v);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
239
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
240 endwhile
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
241 endwhile
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
242
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
243 # restore original shapes
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
244 x = reshape (x, xsiz);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
245 fvec = reshape (fvec, fsiz);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
246
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
247 output.iterations = niter;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
248 output.funcCount = niter + 2;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
249
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
250 endfunction
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
251
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
252 # an assistant function that evaluates a function handle and checks for bad
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
253 # results.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
254 function fx = guarded_eval (fun, x)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
255 fx = fun (x);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
256 if (! all (isreal (fx)))
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
257 error ("fsolve:notreal", "fsolve: non-real value encountered");
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
258 elseif (any (isnan (fx)))
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
259 error ("fsolve:isnan", "fsolve: NaN value encountered");
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
260 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
261 endfunction
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
262
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
263 %!function retval = f (p)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
264 %! x = p(1);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
265 %! y = p(2);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
266 %! z = p(3);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
267 %! retval = zeros (3, 1);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
268 %! retval(1) = sin(x) + y**2 + log(z) - 7;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
269 %! retval(2) = 3*x + 2**y -z**3 + 1;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
270 %! retval(3) = x + y + z - 5;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
271 %!test
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
272 %! x_opt = [ 0.599054;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
273 %! 2.395931;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
274 %! 2.005014 ];
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
275 %! tol = 1.0e-5;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
276 %! [x, fval, info] = fsolve (@f, [ 0.5; 2.0; 2.5 ]);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
277 %! assert (info > 0);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
278 %! assert (norm (x - x_opt, 1) < tol);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
279 %! assert (norm (fval) < tol);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
280
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
281 %!function retval = f (p)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
282 %! x = p(1);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
283 %! y = p(2);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
284 %! z = p(3);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
285 %! w = p(4);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
286 %! retval = zeros (4, 1);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
287 %! retval(1) = 3*x + 4*y + exp (z + w) - 1.007;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
288 %! retval(2) = 6*x - 4*y + exp (3*z + w) - 11;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
289 %! retval(3) = x^4 - 4*y^2 + 6*z - 8*w - 20;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
290 %! retval(4) = x^2 + 2*y^3 + z - w - 4;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
291 %!test
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
292 %! x_opt = [ -0.767297326653401, 0.590671081117440, 1.47190018629642, -1.52719341133957 ];
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
293 %! tol = 1.0e-5;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
294 %! [x, fval, info] = fsolve (@f, [-1, 1, 2, -1]);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
295 %! assert (info > 0);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
296 %! assert (norm (x - x_opt, 1) < tol);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
297 %! assert (norm (fval) < tol);