Mercurial > hg > octave-nkf

annotate scripts/optimization/fsolve.m @ 8395:88fd356b0d95

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
optionally allow pivoting in fsolve
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 11 Dec 2008 07:51:40 +0100
parents a762d9daa700
children 4d008d9f0ccf
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");
8395
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
78 pivoting = optimget (options, "Pivoting", false);
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
79 funvalchk = strcmp (optimget (options, "FunValCheck", "off"), "on");
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
80
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
81 if (funvalchk)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
82 # replace fun with a guarded version
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
83 fun = @(x) guarded_eval (fun, x);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
84 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
85
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
86 # These defaults are rather stringent. I think that normally, user prefers
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
87 # accuracy to performance.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
88
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
89 macheps = eps (class (x0));
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
90
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
91 tolx = optimget (options, "TolX", 1e1*macheps);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
92 tolf = optimget (options, "TolFun",1e2*macheps);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
93
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
94 factor = 100;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
95 # FIXME: TypicalX corresponds to user scaling (???)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
96 autodg = true;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
97
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
98 niter = 1; nfev = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
99
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
100 x = x0(:);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
101 info = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
102
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
103 # outer loop
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
104 while (niter < maxiter && nfev < maxfev && ! info)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
105
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
106 # calc func value and jacobian (possibly via FD)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
107 # handle arbitrary shapes of x and f and remember them
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
108 if (has_jac)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
109 [fvec, fjac] = fcn (reshape (x, xsiz));
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
110 nfev ++;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
111 else
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
112 [fvec, fjac] = __fdjac__ (fcn, reshape (x, xsiz));
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
113 nfev += 1 + length (x);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
114 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
115 fsiz = size (fvec);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
116 fvec = fvec(:);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
117 fn = norm (fvec);
8383
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
118 m = length (fvec);
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
119 n = length (x);
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
120 if (m < n)
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
121 error ("fsolve:under", "cannot solve underdetermined systems");
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
122 elseif (m > n && niter == 1)
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
123 if (isempty (optimget (options, "TolFun")))
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
124 warning (["An overdetermined system cannot usually be solved exactly. ", ...
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
125 "Consider specifying the TolFun option."]);
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
126 endif
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
127 endif
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
128
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
129 # get QR factorization of the jacobian
8395
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
130 if (pivoting)
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
131 [q, r, p] = qr (fjac);
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
132 else
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
133 [q, r] = qr (fjac);
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
134 endif
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
135
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
136 # get column norms, use them as scaling factor
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
137 jcn = norm (fjac, 'columns').';
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
138 if (niter == 1)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
139 if (autodg)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
140 dg = jcn;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
141 dg(dg == 0) = 1;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
142 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
143 xn = norm (dg .* x);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
144 delta = factor * xn;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
145 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
146
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
147 # rescale if necessary
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
148 if (autodg)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
149 dg = max (dg, jcn);
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 nfail = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
153 nsuc = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
154 # inner loop
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
155 while (niter <= maxiter && nfev < maxfev && ! info)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
156
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
157 qtf = q'*fvec;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
158
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
159 # get TR model (dogleg) minimizer
8383
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
160 # in case of an overdetermined system, get lsq solution
8395
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
161 s = - __dogleg__ (r(1:n,:), qtf(1:n), dg, delta);
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
162 if (pivoting)
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
163 s = p * s;
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
164 endif
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
165 sn = norm (dg .* s);
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
166
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
167 if (niter == 1)
8395
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
168 delta = min (delta, sn);
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
169 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
170
8395
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
171 fvec1 = fcn (reshape (x + s, xsiz)) (:);
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
172 fn1 = norm (fvec1);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
173
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
174 if (fn1 < fn)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
175 # scaled actual reduction
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
176 actred = 1 - (fn1/fn)^2;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
177 else
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
178 actred = -1;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
179 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
180
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
181 # scaled predicted reduction, and ratio
8395
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
182 if (pivoting)
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
183 w = qtf + r*(p'*s);
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
184 else
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
185 w = qtf + r*s;
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
186 endif
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
187 t = norm (w);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
188 if (t < fn)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
189 prered = 1 - (t/fn)^2;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
190 ratio = actred / prered;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
191 else
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
192 prered = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
193 ratio = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
194 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
195
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
196 # update delta
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
197 if (ratio < 0.1)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
198 nsuc = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
199 nfail ++;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
200 delta *= 0.5;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
201 if (delta <= sqrt (macheps)*xn)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
202 # trust region became uselessly small
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
203 info = -3;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
204 break;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
205 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
206 else
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
207 nfail = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
208 nsuc ++;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
209 if (abs (1-ratio) <= 0.1)
8395
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
210 delta = 2*sn;
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
211 elseif (ratio >= 0.5 || nsuc > 1)
8395
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
212 delta = max (delta, 2*sn);
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
213 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
214 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
215
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
216 if (ratio >= 1e-4)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
217 # successful iteration
8395
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
218 x += s;
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
219 xn = norm (dg .* x);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
220 fvec = fvec1;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
221 fn = fn1;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
222 niter ++;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
223 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
224
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
225 # Tests for termination conditions. A mysterious place, anything can
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
226 # happen if you change something here...
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 # 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
229 # 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
230 # default value. But 0 usually means uselessly long iterations,
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
231 # so we need scaling-independent tolerances wherever possible.
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 # 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
234 # 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
235 # eps we demand as much accuracy as we can expect.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
236 if (fn <= tolf*n*xn)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
237 info = 1;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
238 # The following tests done only after successful step.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
239 elseif (actred > 0)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
240 # 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
241 # compare to scaled step, so nothing bad.
8395
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
242 if (sn <= tolx*xn)
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
243 info = 2;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
244 # 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
245 # 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
246 elseif (actred < tolf)
8383
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
247 info = 3;
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
248 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
249 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
250
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
251 # criterion for recalculating jacobian
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
252 if (nfail == 2)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
253 break;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
254 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
255
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
256 # compute the scaled Broyden update
8395
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
257 u = (fvec1 - q*w) / sn;
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
258 v = dg .* ((dg .* s) / sn);
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
259 if (pivoting)
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
260 v = p'*v;
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
261 endif
8306
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 # update the QR factorization
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
264 [q, r] = qrupdate (q, r, u, v);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
265
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
266 endwhile
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
267 endwhile
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
268
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
269 # restore original shapes
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
270 x = reshape (x, xsiz);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
271 fvec = reshape (fvec, fsiz);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
272
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
273 output.iterations = niter;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
274 output.funcCount = niter + 2;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
275
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
276 endfunction
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
277
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
278 # 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
279 # results.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
280 function fx = guarded_eval (fun, x)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
281 fx = fun (x);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
282 if (! all (isreal (fx)))
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
283 error ("fsolve:notreal", "fsolve: non-real value encountered");
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
284 elseif (any (isnan (fx)))
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
285 error ("fsolve:isnan", "fsolve: NaN value encountered");
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
286 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
287 endfunction
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
288
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
289 %!function retval = f (p)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
290 %! x = p(1);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
291 %! y = p(2);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
292 %! z = p(3);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
293 %! retval = zeros (3, 1);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
294 %! retval(1) = sin(x) + y**2 + log(z) - 7;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
295 %! retval(2) = 3*x + 2**y -z**3 + 1;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
296 %! retval(3) = x + y + z - 5;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
297 %!test
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
298 %! x_opt = [ 0.599054;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
299 %! 2.395931;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
300 %! 2.005014 ];
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
301 %! tol = 1.0e-5;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
302 %! [x, fval, info] = fsolve (@f, [ 0.5; 2.0; 2.5 ]);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
303 %! assert (info > 0);
8383
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
304 %! assert (norm (x - x_opt, Inf) < tol);
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
305 %! assert (norm (fval) < tol);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
307 %!function retval = f (p)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
308 %! x = p(1);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
309 %! y = p(2);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
310 %! z = p(3);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
311 %! w = p(4);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
312 %! retval = zeros (4, 1);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
313 %! retval(1) = 3*x + 4*y + exp (z + w) - 1.007;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
314 %! retval(2) = 6*x - 4*y + exp (3*z + w) - 11;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
315 %! 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
316 %! retval(4) = x^2 + 2*y^3 + z - w - 4;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
317 %!test
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
318 %! x_opt = [ -0.767297326653401, 0.590671081117440, 1.47190018629642, -1.52719341133957 ];
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
319 %! tol = 1.0e-5;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
320 %! [x, fval, info] = fsolve (@f, [-1, 1, 2, -1]);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
321 %! assert (info > 0);
8383
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
322 %! assert (norm (x - x_opt, Inf) < tol);
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
323 %! assert (norm (fval) < tol);
8383
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
324
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
325 %!function retval = f (p)
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
326 %! x = p(1);
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
327 %! y = p(2);
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
328 %! z = p(3);
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
329 %! retval = zeros (3, 1);
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
330 %! retval(1) = sin(x) + y**2 + log(z) - 7;
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
331 %! retval(2) = 3*x + 2**y -z**3 + 1;
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
332 %! retval(3) = x + y + z - 5;
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
333 %! retval(4) = x*x + y - z*log(z) - 1.36;
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
334 %!test
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
335 %! x_opt = [ 0.599054;
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
336 %! 2.395931;
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
337 %! 2.005014 ];
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
338 %! tol = 1.0e-5;
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
339 %! [x, fval, info] = fsolve (@f, [ 0.5; 2.0; 2.5 ], optimset ("TolFun", 1e-6));
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
340 %! assert (info > 0);
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
341 %! assert (norm (x - x_opt, Inf) < tol);
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
342 %! assert (norm (fval) < tol);
a762d9daa700 allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8306
diff changeset
343
8395
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
344 %!function retval = f (p)
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
345 %! x = p(1);
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
346 %! y = p(2);
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
347 %! z = p(3);
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
348 %! retval = zeros (3, 1);
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
349 %! retval(1) = sin(x) + y**2 + log(z) - 7;
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
350 %! retval(2) = 3*x + 2**y -z**3 + 1;
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
351 %! retval(3) = x + y + z - 5;
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
352 %!test
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
353 %! x_opt = [ 0.599054;
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
354 %! 2.395931;
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
355 %! 2.005014 ];
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
356 %! tol = 1.0e-5;
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
357 %! opt = optimset ("Pivoting", true);
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
358 %! [x, fval, info] = fsolve (@f, [ 0.5; 2.0; 2.5 ], opt);
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
359 %! assert (info > 0);
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
360 %! assert (norm (x - x_opt, Inf) < tol);
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
361 %! assert (norm (fval) < tol);
88fd356b0d95 optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents: 8383
diff changeset
362