annotate scripts/optimization/fzero.m @ 8592:dacfd030633a

handle sparse jacobians in fsolve
author Jaroslav Hajek <highegg@gmail.com>
date Mon, 26 Jan 2009 15:46:33 +0100
parents ec2715c76039
children 8833c0b18eb2
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
1 ## Copyright (C) 2008 VZLU Prague, a.s.
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
2 ##
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
3 ## This file is part of Octave.
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
4 ##
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
5 ## Octave is free software; you can redistribute it and/or modify it
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
6 ## under the terms of the GNU General Public License as published by
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
7 ## the Free Software Foundation; either version 3 of the License, or (at
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
8 ## your option) any later version.
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
9 ##
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
10 ## Octave is distributed in the hope that it will be useful, but
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
13 ## General Public License for more details.
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
14 ##
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
15 ## You should have received a copy of the GNU General Public License
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
16 ## along with Octave; see the file COPYING. If not, see
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
17 ## <http://www.gnu.org/licenses/>.
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
18 ##
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
19 ## Author: Jaroslav Hajek <highegg@gmail.com>
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
20
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
21 ## -*- texinfo -*-
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
22 ## @deftypefn{Function File}{[@var{x}, @var{fval}, @var{info}, @var{output}] =} fzero (@var{fun}, @var{x0}, @var{options})
8515
ec2715c76039 fzero.m, fsolve.m: additional doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
diff changeset
23 ## Find a zero point of a univariate function. @var{fun} should be a function
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
24 ## handle or name. @var{x0} specifies a starting point. @var{options} is a
8515
ec2715c76039 fzero.m, fsolve.m: additional doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
diff changeset
25 ## structure specifying additional options. Currently, @code{fzero}
ec2715c76039 fzero.m, fsolve.m: additional doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
diff changeset
26 ## recognizes these options: @code{"FunValCheck"}, @code{"OutputFcn"},
ec2715c76039 fzero.m, fsolve.m: additional doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
diff changeset
27 ## @code{"TolX"}, @code{"MaxIter"}, @code{"MaxFunEvals"}.
ec2715c76039 fzero.m, fsolve.m: additional doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
diff changeset
28 ## For description of these options, see @ref{doc-optimset,,optimset}.
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
29 ##
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
30 ## On exit, the function returns @var{x}, the approximate zero point
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
31 ## and @var{fval}, the function value thereof.
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
32 ## @var{info} is an exit flag that can have these values:
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
33 ## @itemize
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
34 ## @item 1
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
35 ## The algorithm converged to a solution.
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
36 ## @item 0
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
37 ## Maximum number of iterations or function evaluations has been exhausted.
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
38 ## @item -1
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
39 ## The algorithm has been terminated from user output function.
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
40 ## @item -2
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
41 ## A general unexpected error.
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
42 ## @item -3
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
43 ## A non-real value encountered.
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
44 ## @item -4
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
45 ## A NaN value encountered.
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
46 ## @end itemize
8515
ec2715c76039 fzero.m, fsolve.m: additional doc fixes
John W. Eaton <jwe@octave.org>
parents: 8507
diff changeset
47 ## @seealso{optimset, fsolve}
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
48 ## @end deftypefn
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
49
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
50 ## This is essentially the ACM algorithm 748: Enclosing Zeros of
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
51 ## Continuous Functions due to Alefeld, Potra and Shi, ACM Transactions
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
52 ## on Mathematical Software, Vol. 21, No. 3, September 1995. Although
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
53 ## the workflow should be the same, the structure of the algorithm has
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
54 ## been transformed non-trivially; instead of the authors' approach of
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
55 ## sequentially calling building blocks subprograms we implement here a
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
56 ## FSM version using one interior point determination and one bracketing
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
57 ## per iteration, thus reducing the number of temporary variables and
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
58 ## simplifying the algorithm structure. Further, this approach reduces
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
59 ## the need for external functions and error handling. The algorithm has
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
60 ## also been slightly modified.
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
61
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
62 function [x, fval, info, output] = fzero (fun, x0, options = struct ())
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
63
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
64 if (nargin < 2 || nargin > 3)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
65 print_usage ();
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
66 endif
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
67
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
68 if (ischar (fun))
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
69 fun = str2func (fun);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
70 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
71
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
72 ## TODO
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
73 ## displev = optimget (options, "Display", "notify");
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
74 funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
75 outfcn = optimget (options, "OutputFcn");
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
76 tolx = optimget (options, "TolX", 0);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
77 maxiter = optimget (options, "MaxIter", Inf);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
78 maxfev = optimget (options, "MaxFunEvals", Inf);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
79
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
80 persistent mu = 0.5;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
81
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
82 if (funvalchk)
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
83 ## Replace fun with a guarded version.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
84 fun = @(x) guarded_eval (fun, x);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
85 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
86
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
87 ## The default exit flag if exceeded number of iterations.
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
88 info = 0;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
89 niter = 0;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
90 nfev = 0;
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
91
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
92 x = fval = a = fa = b = fb = NaN;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
93
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
94 ## Prepare...
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
95 a = x0(1);
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
96 fa = fun (a);
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
97 nfev = 1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
98 if (length (x0) > 1)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
99 b = x0(2);
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
100 fb = fun (b);
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
101 nfev += 1;
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
102 else
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
103 ## Try to get b.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
104 if (a == 0)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
105 aa = 1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
106 else
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
107 aa = a;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
108 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
109 for b = [0.9*aa, 1.1*aa, aa-1, aa+1, 0.5*aa 1.5*aa, -aa, 2*aa, -10*aa, 10*aa]
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
110 fb = fun (b); nfev += 1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
111 if (sign (fa) * sign (fb) <= 0)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
112 break;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
113 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
114 endfor
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
115 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
116
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
117 if (b < a)
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
118 u = a;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
119 a = b;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
120 b = u;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
121
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
122 fu = fa;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
123 fa = fb;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
124 fb = fu;
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
125 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
126
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
127 if (! (sign (fa) * sign (fb) <= 0))
8468
866492035ecf fsolve.m, fzero.m: undo part of previous change
John W. Eaton <jwe@octave.org>
parents: 8467
diff changeset
128 error ("fzero:bracket", "fzero: not a valid initial bracketing");
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
129 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
130
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
131 itype = 1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
132
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
133 if (abs (fa) < abs (fb))
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
134 u = a; fu = fa;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
135 else
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
136 u = b; fu = fb;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
137 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
138
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
139 d = e = u;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
140 fd = fe = fu;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
141 mba = mu*(b - a);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
142 while (niter < maxiter && nfev < maxfev)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
143 switch (itype)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
144 case 1
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
145 ## The initial test.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
146 if (b - a <= 2*(2 * abs (u) * eps + tolx))
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
147 x = u; fval = fu;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
148 info = 1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
149 break;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
150 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
151 if (abs (fa) <= 1e3*abs (fb) && abs (fb) <= 1e3*abs (fa))
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
152 ## Secant step.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
153 c = u - (a - b) / (fa - fb) * fu;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
154 else
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
155 ## Bisection step.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
156 c = 0.5*(a + b);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
157 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
158 d = u; fd = fu;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
159 itype = 5;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
160 case {2, 3}
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
161 l = length (unique ([fa, fb, fd, fe]));
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
162 if (l == 4)
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
163 ## Inverse cubic interpolation.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
164 q11 = (d - e) * fd / (fe - fd);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
165 q21 = (b - d) * fb / (fd - fb);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
166 q31 = (a - b) * fa / (fb - fa);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
167 d21 = (b - d) * fd / (fd - fb);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
168 d31 = (a - b) * fb / (fb - fa);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
169 q22 = (d21 - q11) * fb / (fe - fb);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
170 q32 = (d31 - q21) * fa / (fd - fa);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
171 d32 = (d31 - q21) * fd / (fd - fa);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
172 q33 = (d32 - q22) * fa / (fe - fa);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
173 c = a + q31 + q32 + q33;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
174 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
175 if (l < 4 || sign (c - a) * sign (c - b) > 0)
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
176 ## Quadratic interpolation + newton.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
177 a0 = fa;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
178 a1 = (fb - fa)/(b - a);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
179 a2 = ((fd - fb)/(d - b) - a1) / (d - a);
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
180 ## Modification 1: this is simpler and does not seem to be worse.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
181 c = a - a0/a1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
182 if (a2 != 0)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
183 c = a - a0/a1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
184 for i = 1:itype
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
185 pc = a0 + (a1 + a2*(c - b))*(c - a);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
186 pdc = a1 + a2*(2*c - a - b);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
187 if (pdc == 0)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
188 c = a - a0/a1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
189 break;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
190 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
191 c -= pc/pdc;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
192 endfor
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
193 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
194 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
195 itype += 1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
196 case 4
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
197 ## Double secant step.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
198 c = u - 2*(b - a)/(fb - fa)*fu;
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
199 ## Bisect if too far.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
200 if (abs (c - u) > 0.5*(b - a))
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
201 c = 0.5 * (b + a);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
202 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
203 itype = 5;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
204 case 5
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
205 ## Bisection step.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
206 c = 0.5 * (b + a);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
207 itype = 2;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
208 endswitch
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
209
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
210 ## Don't let c come too close to a or b.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
211 delta = 2*0.7*(2 * abs (u) * eps + tolx);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
212 if ((b - a) <= 2*delta)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
213 c = (a + b)/2;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
214 else
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
215 c = max (a + delta, min (b - delta, c));
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
216 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
217
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
218 ## Calculate new point.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
219 x = c;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
220 fval = fc = fun (c);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
221 niter ++; nfev ++;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
222
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
223 ## Modification 2: skip inverse cubic interpolation if
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
224 ## nonmonotonicity is detected.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
225 if (sign (fc - fa) * sign (fc - fb) >= 0)
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
226 ## The new point broke monotonicity.
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
227 ## Disable inverse cubic.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
228 fe = fc;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
229 else
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
230 e = d; fe = fd;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
231 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
232
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
233 ## Bracketing.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
234 if (sign (fa) * sign (fc) < 0)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
235 d = b; fd = fb;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
236 b = c; fb = fc;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
237 elseif (sign (fb) * sign (fc) < 0)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
238 d = a; fd = fa;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
239 a = c; fa = fc;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
240 elseif (fc == 0)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
241 a = b = c; fa = fb = fc;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
242 info = 1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
243 break;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
244 else
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
245 ## This should never happen.
8468
866492035ecf fsolve.m, fzero.m: undo part of previous change
John W. Eaton <jwe@octave.org>
parents: 8467
diff changeset
246 error ("fzero:bracket", "fzero: zero point is not bracketed");
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
247 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
248
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
249 ## If there's an output function, use it now.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
250 if (outfcn)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
251 optv.funccount = niter + 2;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
252 optv.fval = fval;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
253 optv.iteration = niter;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
254 if (outfcn (x, optv, "iter"))
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
255 info = -1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
256 break;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
257 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
258 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
259
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
260 if (abs (fa) < abs (fb))
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
261 u = a; fu = fa;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
262 else
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
263 u = b; fu = fb;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
264 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
265 if (b - a <= 2*(2 * abs (u) * eps + tolx))
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
266 info = 1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
267 break;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
268 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
269
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
270 ## Skip bisection step if successful reduction.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
271 if (itype == 5 && (b - a) <= mba)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
272 itype = 2;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
273 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
274 if (itype == 2)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
275 mba = mu * (b - a);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
276 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
277 endwhile
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
278
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
279 output.iterations = niter;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
280 output.funcCount = niter + 2;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
281 output.bracket = [a, b];
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
282 output.bracketf = [fa, fb];
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
283
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
284 endfunction
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
285
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
286 ## An assistant function that evaluates a function handle and checks for
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
287 ## bad results.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
288 function fx = guarded_eval (fun, x)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
289 fx = fun (x);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
290 fx = fx(1);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
291 if (! isreal (fx))
8468
866492035ecf fsolve.m, fzero.m: undo part of previous change
John W. Eaton <jwe@octave.org>
parents: 8467
diff changeset
292 error ("fzero:notreal", "fzero: non-real value encountered");
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
293 elseif (isnan (fx))
8468
866492035ecf fsolve.m, fzero.m: undo part of previous change
John W. Eaton <jwe@octave.org>
parents: 8467
diff changeset
294 error ("fzero:isnan", "fzero: NaN value encountered");
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
295 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
296 endfunction
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
297
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
298 %!assert(fzero(@cos, [0, 3]), pi/2, 10*eps)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
299 %!assert(fzero(@(x) x^(1/3) - 1e-8, [0,1]), 1e-24, 1e-22*eps)