octave-lyh: scripts/optimization/fsolve.m annotate

annotate scripts/optimization/fsolve.m @ 8507:cadc73247d65

style fixes

author	John W. Eaton <jwe@octave.org>
date	Tue, 13 Jan 2009 14:08:36 -0500
parents	866492035ecf
children	352d3245d4c1

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
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	21 ## -- texinfo --
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	22 ## @deftypefn{Function File} {} fsolve(@var{fcn}, @var{x0}, @var{options})
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	23 ## @deftypefnx{Function File} {[@var{x}, @var{fvec}, @var{info}, @var{output}, @var{fjac}]} = fsolve (@var{fcn}, @dots{})
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	24 ## Solves a system of nonlinear equations defined by the function @var{fcn}.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	25 ## @var{fcn} should accepts a vector (array) defining the unknown variables,
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	26 ## and return a vector of left-hand sides of the equations. Right-hand sides
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	27 ## are defined to be zeros.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	28 ## In other words, this function attempts to determine a vector @var{X} such
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	29 ## that @code{@var{fcn}(@var{X})} gives (approximately) all zeros.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	30 ## @var{x0} determines a starting guess. The shape of @var{x0} is preserved
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	31 ## in all calls to @var{fcn}, but otherwise it is treated as a column vector.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	32 ## @var{options} is a structure specifying additional options. Currently, fsolve
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	33 ## recognizes these options: FunValCheck, OutputFcn, TolX, TolFun, MaxIter,
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	34 ## MaxFunEvals and Jacobian.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	35 ##
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	36 ## If Jacobian is 'on', it specifies that @var{fcn}, called with 2 output arguments,
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	37 ## also returns the Jacobian matrix of right-hand sides at the requested point.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	38 ## TolX specifies the termination tolerance in the unknown variables, while
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	39 ## TolFun is a tolerance for equations. Default is @code{1e1*eps}
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	40 ## for TolX and @code{1e2*eps} for TolFun.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	41 ## For description of the other options, see @code{optimset}.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	42 ##
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	43 ## On return, @var{fval} contains the value of the function @var{fcn}
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	44 ## evaluated at @var{x}, and @var{info} may be one of the following values:
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	45 ##
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	46 ## @table @asis
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	47 ## @item 1
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	48 ## Converged to a solution point. Relative residual error is less than specified
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	49 ## by TolFun.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	50 ## @item 2
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	51 ## Last relative step size was less that TolX.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	52 ## @item 3
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	53 ## Last relative decrease in residual was less than TolF.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	54 ## @item 0
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	55 ## Iteration limit exceeded.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	56 ## @item -3
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	57 ## The trust region radius became excessively small.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	58 ## @end table
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	59 ##
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	60 ## Note: If you only have a single nonlinear equation of one variable, using
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	61 ## @code{fzero} is usually a much better idea.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	62 ## @seealso{fzero,optimset}
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	63 ## @end deftypefn
8306 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
8467 77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options John W. Eaton <jwe@octave.org> parents: 8466 diff changeset	74 has_jac = strcmpi (optimget (options, "Jacobian", "off"), "on");
8306 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);
8467 77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options John W. Eaton <jwe@octave.org> parents: 8466 diff changeset	79 funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
8306 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)
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	82 ## Replace fun with a guarded version.
8306 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
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	86 ## These defaults are rather stringent. I think that normally, user
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	87 ## prefers accuracy to performance.
8306 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);
8467 77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options John W. Eaton <jwe@octave.org> parents: 8466 diff changeset	92 tolf = optimget (options, "TolFun", 1e2*macheps);
8306 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;
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	95 ## FIXME: TypicalX corresponds to user scaling (???)
8306 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
8467 77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options John W. Eaton <jwe@octave.org> parents: 8466 diff changeset	98 niter = 1;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options John W. Eaton <jwe@octave.org> parents: 8466 diff changeset	99 nfev = 0;
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	100
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	101 x = x0(:);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	102 info = 0;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	103
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	104 ## Outer loop.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	105 while (niter < maxiter && nfev < maxfev && ! info)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	106
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	107 ## Calculate function value and Jacobian (possibly via FD).
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	108 ## Handle arbitrary shapes of x and f and remember them.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	109 if (has_jac)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	110 [fvec, fjac] = fcn (reshape (x, xsiz));
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	111 nfev ++;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	112 else
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	113 [fvec, fjac] = __fdjac__ (fcn, reshape (x, xsiz));
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	114 nfev += 1 + length (x);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	115 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	116 fsiz = size (fvec);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	117 fvec = fvec(:);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	118 fn = norm (fvec);
8383 a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	119 m = length (fvec);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	120 n = length (x);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	121 if (m < n)
8468 866492035ecf fsolve.m, fzero.m: undo part of previous change John W. Eaton <jwe@octave.org> parents: 8467 diff changeset	122 error ("fsolve:under", "cannot solve underdetermined systems");
8383 a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	123 elseif (m > n && niter == 1)
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	124 if (isempty (optimget (options, "TolFun")))
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	125 warning ("an overdetermined system cannot usually be solved exactly; consider specifying the TolFun option");
8383 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
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 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
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	136 ## Get column norms, use them as scaling factor.
8306 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
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	147 ## Rescale if necessary.
8306 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;
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	154 ## Inner loop.
8306 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
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	159 ## Get TR model (dogleg) minimizer
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 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)
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	175 ## Scaled actual reduction.
8306 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
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 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
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	196 ## Update delta.
8306 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)
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	202 ## Trust region became uselessly small.
8306 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)
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 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
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	225 ## Tests for termination conditions. A mysterious place, anything
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	226 ## can happen if you change something here...
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	227
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	228 ## The rule of thumb (which I'm not sure M*b is quite following)
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	229 ## is that for a tolerance that depends on scaling, only 0 makes
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	230 ## sense as a default value. But 0 usually means uselessly long
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	231 ## iterations, so we need scaling-independent tolerances wherever
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	232 ## possible.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	233
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	234 ## FIXME -- why tolfnxn? If abs (e) ~ abs(x) * eps is a vector
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	235 ## of perturbations of x, then norm (fjace) <= epsn*xn, i.e. by
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	236 ## 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	237 if (fn <= tolfnxn)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	238 info = 1;
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	239 ## The following tests done only after successful step.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	240 elseif (actred > 0)
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	241 ## This one is classic. Note that we use scaled variables again,
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	242 ## but compare to scaled step, so nothing bad.
8395 88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	243 if (sn <= tolx*xn)
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	244 info = 2;
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	245 ## Again a classic one. It seems weird to use the same tolf
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	246 ## for two different tests, but that's what M*b manual appears
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	247 ## to say.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	248 elseif (actred < tolf)
8383 a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	249 info = 3;
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	250 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	251 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	252
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	253 ## Criterion for recalculating jacobian.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	254 if (nfail == 2)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	255 break;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	256 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	257
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	258 ## Compute the scaled Broyden update.
8395 88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	259 u = (fvec1 - q*w) / sn;
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	260 v = dg .* ((dg .* s) / sn);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	261 if (pivoting)
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	262 v = p'*v;
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	263 endif
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	264
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	265 ## Update the QR factorization.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	266 [q, r] = qrupdate (q, r, u, v);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	267
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	268 endwhile
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	269 endwhile
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	270
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	271 ## Restore original shapes.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	272 x = reshape (x, xsiz);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	273 fvec = reshape (fvec, fsiz);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	274
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	275 output.iterations = niter;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	276 output.funcCount = niter + 2;
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 endfunction
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	279
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	280 ## An assistant function that evaluates a function handle and checks for
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	281 ## bad results.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	282 function fx = guarded_eval (fun, x)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	283 fx = fun (x);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	284 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	285 error ("fsolve:notreal", "fsolve: non-real value encountered");
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	286 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	287 error ("fsolve:isnan", "fsolve: NaN value encountered");
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	288 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	289 endfunction
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	290
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	291 %!function retval = f (p)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	292 %! x = p(1);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	293 %! y = p(2);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	294 %! z = p(3);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	295 %! retval = zeros (3, 1);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	296 %! retval(1) = sin(x) + y**2 + log(z) - 7;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	297 %! retval(2) = 3x + 2y -z*3 + 1;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	298 %! retval(3) = x + y + z - 5;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	299 %!test
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	300 %! x_opt = [ 0.599054;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	301 %! 2.395931;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	302 %! 2.005014 ];
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	303 %! tol = 1.0e-5;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	304 %! [x, fval, info] = fsolve (@f, [ 0.5; 2.0; 2.5 ]);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	305 %! assert (info > 0);
8383 a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	306 %! assert (norm (x - x_opt, Inf) < tol);
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	307 %! assert (norm (fval) < tol);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	308
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	309 %!function retval = f (p)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	310 %! x = p(1);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	311 %! y = p(2);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	312 %! z = p(3);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	313 %! w = p(4);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	314 %! retval = zeros (4, 1);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	315 %! retval(1) = 3x + 4y + exp (z + w) - 1.007;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	316 %! retval(2) = 6x - 4y + exp (3*z + w) - 11;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	317 %! retval(3) = x^4 - 4y^2 + 6z - 8*w - 20;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	318 %! retval(4) = x^2 + 2*y^3 + z - w - 4;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	319 %!test
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	320 %! x_opt = [ -0.767297326653401, 0.590671081117440, 1.47190018629642, -1.52719341133957 ];
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	321 %! tol = 1.0e-5;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	322 %! [x, fval, info] = fsolve (@f, [-1, 1, 2, -1]);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	323 %! assert (info > 0);
8383 a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	324 %! assert (norm (x - x_opt, Inf) < tol);
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	325 %! assert (norm (fval) < tol);
8383 a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	326
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	327 %!function retval = f (p)
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	328 %! x = p(1);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	329 %! y = p(2);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	330 %! z = p(3);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	331 %! retval = zeros (3, 1);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	332 %! retval(1) = sin(x) + y**2 + log(z) - 7;
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	333 %! retval(2) = 3x + 2y -z*3 + 1;
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	334 %! retval(3) = x + y + z - 5;
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	335 %! retval(4) = xx + y - zlog(z) - 1.36;
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	336 %!test
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	337 %! x_opt = [ 0.599054;
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	338 %! 2.395931;
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	339 %! 2.005014 ];
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	340 %! tol = 1.0e-5;
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	341 %! [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	342 %! assert (info > 0);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	343 %! assert (norm (x - x_opt, Inf) < tol);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	344 %! assert (norm (fval) < tol);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	345
8395 88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	346 %!function retval = f (p)
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	347 %! x = p(1);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	348 %! y = p(2);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	349 %! z = p(3);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	350 %! retval = zeros (3, 1);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	351 %! retval(1) = sin(x) + y**2 + log(z) - 7;
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	352 %! retval(2) = 3x + 2y -z*3 + 1;
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	353 %! retval(3) = x + y + z - 5;
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	354 %!test
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	355 %! x_opt = [ 0.599054;
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	356 %! 2.395931;
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	357 %! 2.005014 ];
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	358 %! tol = 1.0e-5;
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	359 %! opt = optimset ("Pivoting", true);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	360 %! [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	361 %! assert (info > 0);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	362 %! assert (norm (x - x_opt, Inf) < tol);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	363 %! assert (norm (fval) < tol);

Mercurial > hg > octave-lyh

annotate scripts/optimization/fsolve.m @ 8507:cadc73247d65