Mercurial > hg > octave-nkf
annotate scripts/optimization/fsolve.m @ 8592:dacfd030633a
handle sparse jacobians in fsolve
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Mon, 26 Jan 2009 15:46:33 +0100 |
| parents | c136d313206a |
| children | 8833c0b18eb2 |
| rev | line source |
|---|---|
| 8590 | 1 ## Copyright (C) 2008, 2009 VZLU Prague, a.s. |
|
8306
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 |
| 8466 | 21 ## -*- texinfo -*- |
| 8513 | 22 ## @deftypefn{Function File} {} fsolve (@var{fcn}, @var{x0}, @var{options}) |
| 8466 | 23 ## @deftypefnx{Function File} {[@var{x}, @var{fvec}, @var{info}, @var{output}, @var{fjac}]} = fsolve (@var{fcn}, @dots{}) |
| 8513 | 24 ## Solve a system of nonlinear equations defined by the function @var{fcn}. |
| 8466 | 25 ## @var{fcn} should accepts a vector (array) defining the unknown variables, |
| 26 ## and return a vector of left-hand sides of the equations. Right-hand sides | |
| 27 ## are defined to be zeros. | |
|
8514
39867b4aca52
fsolve.m: additional doc fix
John W. Eaton <jwe@octave.org>
parents:
8513
diff
changeset
|
28 ## In other words, this function attempts to determine a vector @var{x} such |
|
39867b4aca52
fsolve.m: additional doc fix
John W. Eaton <jwe@octave.org>
parents:
8513
diff
changeset
|
29 ## that @code{@var{fcn} (@var{x})} gives (approximately) all zeros. |
| 8466 | 30 ## @var{x0} determines a starting guess. The shape of @var{x0} is preserved |
| 31 ## in all calls to @var{fcn}, but otherwise it is treated as a column vector. | |
|
8515
ec2715c76039
fzero.m, fsolve.m: additional doc fixes
John W. Eaton <jwe@octave.org>
parents:
8514
diff
changeset
|
32 ## @var{options} is a structure specifying additional options. |
|
ec2715c76039
fzero.m, fsolve.m: additional doc fixes
John W. Eaton <jwe@octave.org>
parents:
8514
diff
changeset
|
33 ## Currently, @code{fsolve} recognizes these options: |
|
ec2715c76039
fzero.m, fsolve.m: additional doc fixes
John W. Eaton <jwe@octave.org>
parents:
8514
diff
changeset
|
34 ## @code{"FunValCheck"}, @code{"OutputFcn"}, @code{"TolX"}, |
| 8590 | 35 ## @code{"TolFun"}, @code{"MaxIter"}, @code{"MaxFunEvals"}, |
| 36 ## @code{"Jacobian"} and @code{"Updating"}. | |
| 8466 | 37 ## |
| 8513 | 38 ## If @code{"Jacobian"} is @code{"on"}, it specifies that @var{fcn}, |
| 39 ## called with 2 output arguments, also returns the Jacobian matrix | |
| 40 ## of right-hand sides at the requested point. @code{"TolX"} specifies | |
| 41 ## the termination tolerance in the unknown variables, while | |
| 42 ## @code{"TolFun"} is a tolerance for equations. Default is @code{1e1*eps} | |
| 43 ## for @code{"TolX"} and @code{1e2*eps} for @code{"TolFun"}. | |
| 8590 | 44 ## If @code{"Updating"} is true, the function will attempt to use Broyden |
| 45 ## updates to update the Jacobian, in order to reduce the amount of jacobian | |
| 46 ## calculations. | |
| 47 ## If your user function always calculates the Jacobian (regardless of number | |
| 48 ## of output arguments), this option provides no advantage and should be set to | |
| 49 ## false. | |
| 50 ## | |
| 8466 | 51 ## For description of the other options, see @code{optimset}. |
| 52 ## | |
| 53 ## On return, @var{fval} contains the value of the function @var{fcn} | |
| 54 ## evaluated at @var{x}, and @var{info} may be one of the following values: | |
| 55 ## | |
| 56 ## @table @asis | |
| 57 ## @item 1 | |
| 58 ## Converged to a solution point. Relative residual error is less than specified | |
| 59 ## by TolFun. | |
| 60 ## @item 2 | |
| 61 ## Last relative step size was less that TolX. | |
| 62 ## @item 3 | |
| 63 ## Last relative decrease in residual was less than TolF. | |
| 64 ## @item 0 | |
| 65 ## Iteration limit exceeded. | |
| 66 ## @item -3 | |
| 67 ## The trust region radius became excessively small. | |
| 68 ## @end table | |
| 69 ## | |
| 70 ## Note: If you only have a single nonlinear equation of one variable, using | |
| 71 ## @code{fzero} is usually a much better idea. | |
| 8513 | 72 ## @seealso{fzero, optimset} |
| 8466 | 73 ## @end deftypefn |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
74 |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
75 function [x, fvec, info, output, fjac] = fsolve (fcn, x0, options) |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
76 |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
77 if (nargin < 3) |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
78 options = struct (); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
79 endif |
|
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 xsiz = size (x0); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
82 n = numel (x0); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
83 |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8466
diff
changeset
|
84 has_jac = strcmpi (optimget (options, "Jacobian", "off"), "on"); |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
85 maxiter = optimget (options, "MaxIter", Inf); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
86 maxfev = optimget (options, "MaxFunEvals", Inf); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
87 outfcn = optimget (options, "OutputFcn"); |
| 8590 | 88 updating = optimget (options, "Updating", true); |
| 89 | |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8466
diff
changeset
|
90 funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on"); |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
91 |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
92 if (funvalchk) |
| 8507 | 93 ## Replace fun with a guarded version. |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
94 fun = @(x) guarded_eval (fun, x); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
95 endif |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
96 |
| 8466 | 97 ## These defaults are rather stringent. I think that normally, user |
| 98 ## prefers accuracy to performance. | |
|
8306
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 macheps = eps (class (x0)); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
101 |
| 8590 | 102 tolx = optimget (options, "TolX", sqrt (macheps)); |
| 103 tolf = optimget (options, "TolFun", sqrt (macheps)); | |
|
8306
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 factor = 100; |
| 8466 | 106 ## FIXME: TypicalX corresponds to user scaling (???) |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
107 autodg = true; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
108 |
|
8467
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8466
diff
changeset
|
109 niter = 1; |
|
77b8d4aa2743
fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents:
8466
diff
changeset
|
110 nfev = 0; |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
111 |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
112 x = x0(:); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
113 info = 0; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
114 |
| 8507 | 115 ## Outer loop. |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
116 while (niter < maxiter && nfev < maxfev && ! info) |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
117 |
| 8507 | 118 ## Calculate function value and Jacobian (possibly via FD). |
| 119 ## Handle arbitrary shapes of x and f and remember them. | |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
120 if (has_jac) |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
121 [fvec, fjac] = fcn (reshape (x, xsiz)); |
|
8592
dacfd030633a
handle sparse jacobians in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8590
diff
changeset
|
122 ## If the jacobian is sparse, disable Broyden updating. |
|
dacfd030633a
handle sparse jacobians in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8590
diff
changeset
|
123 if (issparse (fjac)) |
|
dacfd030633a
handle sparse jacobians in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8590
diff
changeset
|
124 updating = false; |
|
dacfd030633a
handle sparse jacobians in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8590
diff
changeset
|
125 endif |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
126 nfev ++; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
127 else |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
128 [fvec, fjac] = __fdjac__ (fcn, reshape (x, xsiz)); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
129 nfev += 1 + length (x); |
|
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 fsiz = size (fvec); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
132 fvec = fvec(:); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
133 fn = norm (fvec); |
|
8383
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
134 m = length (fvec); |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
135 n = length (x); |
| 8590 | 136 |
| 137 ## For square and overdetermined systems, we update a (pivoted) QR | |
| 138 ## factorization of the jacobian to avoid solving a full system in each | |
| 139 ## step. In this case, we pass a triangular matrix to __dogleg__. | |
| 140 ## Pivoted QR is used for slightly better robustness and invariance | |
| 141 ## w.r.t. permutations of variables. | |
| 142 useqr = updating && m >= n && n > 10; | |
| 143 | |
| 144 if (useqr) | |
| 145 ## Get QR factorization of the jacobian, optionally with column pivoting. | |
| 146 ## TODO: pivoting is only done in the first step, to get invariance | |
| 147 ## w.r.t. permutations of variables. Maybe it could be beneficial to | |
| 148 ## repivot occassionally? | |
| 149 if (niter == 1) | |
| 150 [q, r, p] = qr (fjac, 0); | |
| 151 ## p is a column vector. Blame Matlab for the inconsistency. | |
| 152 ## Octave can handle permutation matrices gracefully. | |
| 153 p = eye (n)(:, p); | |
| 154 else | |
| 155 [q, r] = qr (fjac*p, 0); | |
|
8383
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
156 endif |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
157 endif |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
158 |
| 8590 | 159 ## Get column norms, use them as scaling factors. |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
160 jcn = norm (fjac, 'columns').'; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
161 if (niter == 1) |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
162 if (autodg) |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
163 dg = jcn; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
164 dg(dg == 0) = 1; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
165 endif |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
166 xn = norm (dg .* x); |
| 8590 | 167 ## FIXME: something better? |
| 168 delta = max (factor * xn, 1); | |
|
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 |
| 8507 | 171 ## Rescale if necessary. |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
172 if (autodg) |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
173 dg = max (dg, jcn); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
174 endif |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
175 |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
176 nfail = 0; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
177 nsuc = 0; |
| 8507 | 178 ## Inner loop. |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
179 while (niter <= maxiter && nfev < maxfev && ! info) |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
180 |
| 8590 | 181 if (useqr) |
| 182 tr_mat = r; | |
| 183 tr_vec = q'*fvec; | |
| 184 else | |
| 185 tr_mat = fjac; | |
| 186 tr_vec = fvec; | |
| 187 endif | |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
188 |
| 8590 | 189 ## Get trust-region model (dogleg) minimizer. |
| 190 if (useqr) | |
| 191 qtf = q'*fvec; | |
| 192 s = - __dogleg__ (r, qtf, dg, delta); | |
| 193 w = qtf + r * s; | |
| 194 s = p * s; | |
| 195 else | |
| 196 s = - __dogleg__ (fjac, fvec, dg, delta); | |
| 197 w = fvec + fjac * s; | |
|
8395
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
198 endif |
| 8590 | 199 |
|
8395
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
200 sn = norm (dg .* s); |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
201 if (niter == 1) |
|
8395
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
202 delta = min (delta, sn); |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
203 endif |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
204 |
|
8395
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
205 fvec1 = fcn (reshape (x + s, xsiz)) (:); |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
206 fn1 = norm (fvec1); |
| 8590 | 207 nfev ++; |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
208 |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
209 if (fn1 < fn) |
| 8507 | 210 ## Scaled actual reduction. |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
211 actred = 1 - (fn1/fn)^2; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
212 else |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
213 actred = -1; |
|
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 |
| 8507 | 216 ## Scaled predicted reduction, and ratio. |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
217 t = norm (w); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
218 if (t < fn) |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
219 prered = 1 - (t/fn)^2; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
220 ratio = actred / prered; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
221 else |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
222 prered = 0; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
223 ratio = 0; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
224 endif |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
225 |
| 8507 | 226 ## Update delta. |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
227 if (ratio < 0.1) |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
228 nsuc = 0; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
229 nfail ++; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
230 delta *= 0.5; |
| 8590 | 231 if (delta <= 1e1*macheps*xn) |
| 8507 | 232 ## Trust region became uselessly small. |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
233 info = -3; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
234 break; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
235 endif |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
236 else |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
237 nfail = 0; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
238 nsuc ++; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
239 if (abs (1-ratio) <= 0.1) |
|
8395
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
240 delta = 2*sn; |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
241 elseif (ratio >= 0.5 || nsuc > 1) |
|
8395
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
242 delta = max (delta, 2*sn); |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
243 endif |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
244 endif |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
245 |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
246 if (ratio >= 1e-4) |
| 8507 | 247 ## Successful iteration. |
|
8395
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
248 x += s; |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
249 xn = norm (dg .* x); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
250 fvec = fvec1; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
251 fn = fn1; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
252 niter ++; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
253 endif |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
254 |
| 8590 | 255 ## FIXME: should outputfcn be only called after a successful iteration? |
| 256 if (outfcn) | |
| 257 optimvalues.iter = niter; | |
| 258 optimvalues.funccount = nfev; | |
| 259 optimvalues.fval = fn; | |
| 260 optimvalues.searchdirection = s; | |
| 261 state = 'iter'; | |
| 262 stop = outfcn (x, optimvalues, state); | |
| 263 if (stop) | |
| 264 info = -1; | |
| 265 break; | |
| 266 endif | |
| 267 endif | |
| 268 | |
| 8466 | 269 ## Tests for termination conditions. A mysterious place, anything |
| 270 ## can happen if you change something here... | |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
271 |
| 8466 | 272 ## The rule of thumb (which I'm not sure M*b is quite following) |
| 273 ## is that for a tolerance that depends on scaling, only 0 makes | |
| 274 ## sense as a default value. But 0 usually means uselessly long | |
| 275 ## iterations, so we need scaling-independent tolerances wherever | |
| 276 ## possible. | |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
277 |
| 8466 | 278 ## FIXME -- why tolf*n*xn? If abs (e) ~ abs(x) * eps is a vector |
| 279 ## of perturbations of x, then norm (fjac*e) <= eps*n*xn, i.e. by | |
| 280 ## tolf ~ eps we demand as much accuracy as we can expect. | |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
281 if (fn <= tolf*n*xn) |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
282 info = 1; |
| 8466 | 283 ## The following tests done only after successful step. |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
284 elseif (actred > 0) |
| 8466 | 285 ## This one is classic. Note that we use scaled variables again, |
| 286 ## but compare to scaled step, so nothing bad. | |
|
8395
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
287 if (sn <= tolx*xn) |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
288 info = 2; |
| 8466 | 289 ## Again a classic one. It seems weird to use the same tolf |
| 290 ## for two different tests, but that's what M*b manual appears | |
| 291 ## to say. | |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
292 elseif (actred < tolf) |
|
8383
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
293 info = 3; |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
294 endif |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
295 endif |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
296 |
| 8507 | 297 ## Criterion for recalculating jacobian. |
| 8590 | 298 if (! updating || nfail == 2) |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
299 break; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
300 endif |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
301 |
| 8507 | 302 ## Compute the scaled Broyden update. |
| 8590 | 303 if (useqr) |
| 304 u = (fvec1 - q*w) / sn; | |
| 305 v = dg .* ((dg .* s) / sn); | |
| 306 v = p' * v; | |
| 307 | |
| 308 ## Update the QR factorization. | |
| 309 [q, r] = qrupdate (q, r, u, v); | |
| 310 else | |
| 311 u = (fvec1 - w); | |
| 312 v = dg .* ((dg .* s) / sn); | |
| 313 | |
| 314 ## update the jacobian | |
| 315 fjac += u * v'; | |
|
8395
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
316 endif |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
317 endwhile |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
318 endwhile |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
319 |
| 8507 | 320 ## Restore original shapes. |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
321 x = reshape (x, xsiz); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
322 fvec = reshape (fvec, fsiz); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
323 |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
324 output.iterations = niter; |
| 8590 | 325 output.funcCount = nfev; |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
326 |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
327 endfunction |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
328 |
| 8466 | 329 ## An assistant function that evaluates a function handle and checks for |
| 330 ## bad results. | |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
331 function fx = guarded_eval (fun, x) |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
332 fx = fun (x); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
333 if (! all (isreal (fx))) |
|
8468
866492035ecf
fsolve.m, fzero.m: undo part of previous change
John W. Eaton <jwe@octave.org>
parents:
8467
diff
changeset
|
334 error ("fsolve:notreal", "fsolve: non-real value encountered"); |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
335 elseif (any (isnan (fx))) |
|
8468
866492035ecf
fsolve.m, fzero.m: undo part of previous change
John W. Eaton <jwe@octave.org>
parents:
8467
diff
changeset
|
336 error ("fsolve:isnan", "fsolve: NaN value encountered"); |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
337 endif |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
338 endfunction |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
339 |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
340 %!function retval = f (p) |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
341 %! x = p(1); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
342 %! y = p(2); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
343 %! z = p(3); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
344 %! retval = zeros (3, 1); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
345 %! retval(1) = sin(x) + y**2 + log(z) - 7; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
346 %! retval(2) = 3*x + 2**y -z**3 + 1; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
347 %! retval(3) = x + y + z - 5; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
348 %!test |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
349 %! x_opt = [ 0.599054; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
350 %! 2.395931; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
351 %! 2.005014 ]; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
352 %! tol = 1.0e-5; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
353 %! [x, fval, info] = fsolve (@f, [ 0.5; 2.0; 2.5 ]); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
354 %! assert (info > 0); |
|
8383
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
355 %! assert (norm (x - x_opt, Inf) < tol); |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
356 %! assert (norm (fval) < tol); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
357 |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
358 %!function retval = f (p) |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
359 %! x = p(1); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
360 %! y = p(2); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
361 %! z = p(3); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
362 %! w = p(4); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
363 %! retval = zeros (4, 1); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
364 %! retval(1) = 3*x + 4*y + exp (z + w) - 1.007; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
365 %! retval(2) = 6*x - 4*y + exp (3*z + w) - 11; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
366 %! 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
|
367 %! retval(4) = x^2 + 2*y^3 + z - w - 4; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
368 %!test |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
369 %! x_opt = [ -0.767297326653401, 0.590671081117440, 1.47190018629642, -1.52719341133957 ]; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
370 %! tol = 1.0e-5; |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
371 %! [x, fval, info] = fsolve (@f, [-1, 1, 2, -1]); |
|
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
372 %! assert (info > 0); |
|
8383
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
373 %! assert (norm (x - x_opt, Inf) < tol); |
|
8306
43795cf108d0
initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff
changeset
|
374 %! assert (norm (fval) < tol); |
|
8383
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
375 |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
376 %!function retval = f (p) |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
377 %! x = p(1); |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
378 %! y = p(2); |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
379 %! z = p(3); |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
380 %! retval = zeros (3, 1); |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
381 %! retval(1) = sin(x) + y**2 + log(z) - 7; |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
382 %! retval(2) = 3*x + 2**y -z**3 + 1; |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
383 %! retval(3) = x + y + z - 5; |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
384 %! 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
|
385 %!test |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
386 %! x_opt = [ 0.599054; |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
387 %! 2.395931; |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
388 %! 2.005014 ]; |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
389 %! tol = 1.0e-5; |
| 8590 | 390 %! [x, fval, info] = fsolve (@f, [ 0.5; 2.0; 2.5 ]); |
|
8383
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
391 %! assert (info > 0); |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
392 %! assert (norm (x - x_opt, Inf) < tol); |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
393 %! assert (norm (fval) < tol); |
|
a762d9daa700
allow overdetermined systems in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8306
diff
changeset
|
394 |
|
8395
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
395 %!function retval = f (p) |
|
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
396 %! x = p(1); |
|
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
397 %! y = p(2); |
|
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
398 %! z = p(3); |
|
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
399 %! retval = zeros (3, 1); |
|
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
400 %! retval(1) = sin(x) + y**2 + log(z) - 7; |
|
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
401 %! retval(2) = 3*x + 2**y -z**3 + 1; |
|
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
402 %! retval(3) = x + y + z - 5; |
|
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
403 %!test |
|
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
404 %! x_opt = [ 0.599054; |
|
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
405 %! 2.395931; |
|
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
406 %! 2.005014 ]; |
|
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
407 %! tol = 1.0e-5; |
| 8590 | 408 %! opt = optimset ("Updating", "qrp"); |
|
8395
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
409 %! [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
|
410 %! assert (info > 0); |
|
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
411 %! assert (norm (x - x_opt, Inf) < tol); |
|
88fd356b0d95
optionally allow pivoting in fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
8383
diff
changeset
|
412 %! assert (norm (fval) < tol); |
| 8590 | 413 |
| 414 %!test | |
| 415 %! b0 = 3; | |
| 416 %! a0 = 0.2; | |
| 417 %! x = 0:.5:5; | |
| 418 %! noise = 1e-5 * sin (100*x); | |
| 419 %! y = exp (-a0*x) + b0 + noise; | |
| 420 %! c_opt = [a0, b0]; | |
| 421 %! tol = 1e-5; | |
| 422 %! | |
| 423 %! [c, fval, info, output] = fsolve (@(c) (exp(-c(1)*x) + c(2) - y), [0, 0]); | |
| 424 %! assert (info > 0); | |
| 425 %! assert (norm (c - c_opt, Inf) < tol); | |
| 426 %! assert (norm (fval) < norm (noise)); | |
| 427 |