Mercurial > hg > octave-nkf

annotate scripts/optimization/fzero.m @ 20830:b65888ec820e draft default tip gccjit

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
dmalcom gcc jit import
author Stefan Mahr <dac922@gmx.de>
date Fri, 27 Feb 2015 16:59:36 +0100
parents f1d0f506ee78
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
19898
4197fc428c7d maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents: 18186
diff changeset
1 ## Copyright (C) 2008-2015 VZLU Prague, a.s.
8305
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 -*-
10687
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
22 ## @deftypefn {Function File} {} fzero (@var{fun}, @var{x0})
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
23 ## @deftypefnx {Function File} {} fzero (@var{fun}, @var{x0}, @var{options})
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
24 ## @deftypefnx {Function File} {[@var{x}, @var{fval}, @var{info}, @var{output}] =} fzero (@dots{})
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
25 ## Find a zero of a univariate function.
10687
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
26 ##
20375
f1d0f506ee78 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20038
diff changeset
27 ## @var{fun} is a function handle, inline function, or string containing the
f1d0f506ee78 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20038
diff changeset
28 ## name of the function to evaluate.
f1d0f506ee78 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20038
diff changeset
29 ##
12632
2dbac27e0e40 doc: miscellaneous touch-ups to documentation strings
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
30 ## @var{x0} should be a two-element vector specifying two points which
2dbac27e0e40 doc: miscellaneous touch-ups to documentation strings
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
31 ## bracket a zero. In other words, there must be a change in sign of the
2dbac27e0e40 doc: miscellaneous touch-ups to documentation strings
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
32 ## function between @var{x0}(1) and @var{x0}(2). More mathematically, the
2dbac27e0e40 doc: miscellaneous touch-ups to documentation strings
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
33 ## following must hold
10821
693e22af08ae Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents: 10687
diff changeset
34 ##
10470
9500a66118dc improve fzero
Jaroslav Hajek <highegg@gmail.com>
parents: 10394
diff changeset
35 ## @example
9500a66118dc improve fzero
Jaroslav Hajek <highegg@gmail.com>
parents: 10394
diff changeset
36 ## sign (@var{fun}(@var{x0}(1))) * sign (@var{fun}(@var{x0}(2))) <= 0
9500a66118dc improve fzero
Jaroslav Hajek <highegg@gmail.com>
parents: 10394
diff changeset
37 ## @end example
10821
693e22af08ae Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents: 10687
diff changeset
38 ##
20375
f1d0f506ee78 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20038
diff changeset
39 ## If @var{x0} is a single scalar then several nearby and distant values are
f1d0f506ee78 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20038
diff changeset
40 ## probed in an attempt to obtain a valid bracketing. If this is not
f1d0f506ee78 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20038
diff changeset
41 ## successful, the function fails.
f1d0f506ee78 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20038
diff changeset
42 ##
f1d0f506ee78 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20038
diff changeset
43 ## @var{options} is a structure specifying additional options. Currently,
f1d0f506ee78 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20038
diff changeset
44 ## @code{fzero} recognizes these options:
f1d0f506ee78 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20038
diff changeset
45 ## @qcode{"FunValCheck"}, @qcode{"OutputFcn"}, @qcode{"TolX"},
f1d0f506ee78 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20038
diff changeset
46 ## @qcode{"MaxIter"}, @qcode{"MaxFunEvals"}.
17097
e7a059a9a644 doc: Use XREF as anchor prefix in documentation for clearer results in Info viewer.
Rik <rik@octave.org>
parents: 16772
diff changeset
47 ## For a description of these options, see @ref{XREFoptimset,,optimset}.
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
48 ##
20375
f1d0f506ee78 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20038
diff changeset
49 ## On exit, the function returns @var{x}, the approximate zero point and
f1d0f506ee78 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20038
diff changeset
50 ## @var{fval}, the function value thereof.
f1d0f506ee78 doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20038
diff changeset
51 ##
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
52 ## @var{info} is an exit flag that can have these values:
10297
ed88ea036716 improve docs of fzero/fminbnd
Jaroslav Hajek <highegg@gmail.com>
parents: 10296
diff changeset
53 ##
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
54 ## @itemize
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
55 ## @item 1
10687
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
56 ## The algorithm converged to a solution.
10821
693e22af08ae Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents: 10687
diff changeset
57 ##
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
58 ## @item 0
10687
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
59 ## Maximum number of iterations or function evaluations has been reached.
10821
693e22af08ae Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents: 10687
diff changeset
60 ##
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
61 ## @item -1
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
62 ## The algorithm has been terminated from user output function.
10821
693e22af08ae Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents: 10687
diff changeset
63 ##
10337
7c97da90fc8f check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents: 10297
diff changeset
64 ## @item -5
7c97da90fc8f check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents: 10297
diff changeset
65 ## The algorithm may have converged to a singular point.
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
66 ## @end itemize
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
67 ##
10687
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
68 ## @var{output} is a structure containing runtime information about the
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
69 ## @code{fzero} algorithm. Fields in the structure are:
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
70 ##
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
71 ## @itemize
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
72 ## @item iterations
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
73 ## Number of iterations through loop.
10821
693e22af08ae Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents: 10687
diff changeset
74 ##
10687
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
75 ## @item nfev
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
76 ## Number of function evaluations.
10821
693e22af08ae Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents: 10687
diff changeset
77 ##
10687
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
78 ## @item bracketx
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
79 ## A two-element vector with the final bracketing of the zero along the x-axis.
10821
693e22af08ae Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents: 10687
diff changeset
80 ##
10687
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
81 ## @item brackety
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
82 ## A two-element vector with the final bracketing of the zero along the y-axis.
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
83 ## @end itemize
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
84 ## @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
85 ## @end deftypefn
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
86
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
87 ## 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
88 ## 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
89 ## 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
90 ## 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
91 ## 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
92 ## 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
93 ## 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
94 ## 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
95 ## 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
96 ## 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
97 ## 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
98
13027
b9a89ca0fb75 prevent optimization functions from setting ans in workspace at startup
John W. Eaton <jwe@octave.org>
parents: 12632
diff changeset
99 ## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup.
b9a89ca0fb75 prevent optimization functions from setting ans in workspace at startup
John W. Eaton <jwe@octave.org>
parents: 12632
diff changeset
100 ## PKG_ADD: [~] = __all_opts__ ("fzero");
8648
ff61b53eb294 optimization: use PKG_ADD: comments instead of PKG_ADD file
John W. Eaton <jwe@octave.org>
parents: 8596
diff changeset
101
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
102 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
103
8596
8833c0b18eb2 enable default settings queries in optim funcs
Jaroslav Hajek <highegg@gmail.com>
parents: 8515
diff changeset
104 ## Get default options if requested.
8833c0b18eb2 enable default settings queries in optim funcs
Jaroslav Hajek <highegg@gmail.com>
parents: 8515
diff changeset
105 if (nargin == 1 && ischar (fun) && strcmp (fun, 'defaults'))
18186
6e0686858b6a Make fzero default tolerance Matlab compatible (bug #41000)
Felipe G. Nievinski <fgnievinski@gmail.com>
parents: 18132
diff changeset
106 x = optimset ("MaxIter", Inf, "MaxFunEvals", Inf, "TolX", eps,
14552
86854d032a37 maint: miscellaneous style fixes for .m files
John W. Eaton <jwe@octave.org>
parents: 14363
diff changeset
107 "OutputFcn", [], "FunValCheck", "off");
8596
8833c0b18eb2 enable default settings queries in optim funcs
Jaroslav Hajek <highegg@gmail.com>
parents: 8515
diff changeset
108 return;
8833c0b18eb2 enable default settings queries in optim funcs
Jaroslav Hajek <highegg@gmail.com>
parents: 8515
diff changeset
109 endif
8833c0b18eb2 enable default settings queries in optim funcs
Jaroslav Hajek <highegg@gmail.com>
parents: 8515
diff changeset
110
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
111 if (nargin < 2 || nargin > 3)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
112 print_usage ();
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
113 endif
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
114
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
115 if (ischar (fun))
9464
e598248a060d safer str2func use in optim functions
Jaroslav Hajek <highegg@gmail.com>
parents: 9143
diff changeset
116 fun = str2func (fun, "global");
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
117 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
118
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
119 ## TODO
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
120 ## 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
121 funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
122 outfcn = optimget (options, "OutputFcn");
18186
6e0686858b6a Make fzero default tolerance Matlab compatible (bug #41000)
Felipe G. Nievinski <fgnievinski@gmail.com>
parents: 18132
diff changeset
123 tolx = optimget (options, "TolX", eps);
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
124 maxiter = optimget (options, "MaxIter", Inf);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
125 maxfev = optimget (options, "MaxFunEvals", Inf);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
126
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
127 persistent mu = 0.5;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
128
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
129 if (funvalchk)
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
130 ## Replace fun with a guarded version.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
131 fun = @(x) guarded_eval (fun, x);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
132 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
133
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
134 ## 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
135 info = 0;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
136 niter = 0;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
137 nfev = 0;
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
138
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
139 x = fval = a = fa = b = fb = NaN;
10296
035ac548a67e implement fminbnd
Jaroslav Hajek <highegg@gmail.com>
parents: 9464
diff changeset
140 eps = eps (class (x0));
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
141
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
142 ## Prepare...
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
143 a = x0(1);
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
144 fa = fun (a);
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
145 nfev = 1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
146 if (length (x0) > 1)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
147 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
148 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
149 nfev += 1;
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
150 else
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
151 ## Try to get b.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
152 if (a == 0)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
153 aa = 1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
154 else
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
155 aa = a;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
156 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
157 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
158 fb = fun (b); nfev += 1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
159 if (sign (fa) * sign (fb) <= 0)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
160 break;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
161 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
162 endfor
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
163 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
164
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
165 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
166 u = a;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
167 a = b;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
168 b = u;
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
169
8467
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
170 fu = fa;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
171 fa = fb;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options
John W. Eaton <jwe@octave.org>
parents: 8305
diff changeset
172 fb = fu;
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
173 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
174
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
175 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
176 error ("fzero:bracket", "fzero: not a valid initial bracketing");
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
177 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
178
10337
7c97da90fc8f check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents: 10297
diff changeset
179 slope0 = (fb - fa) / (b - a);
7c97da90fc8f check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents: 10297
diff changeset
180
10470
9500a66118dc improve fzero
Jaroslav Hajek <highegg@gmail.com>
parents: 10394
diff changeset
181 if (fa == 0)
9500a66118dc improve fzero
Jaroslav Hajek <highegg@gmail.com>
parents: 10394
diff changeset
182 b = a;
9500a66118dc improve fzero
Jaroslav Hajek <highegg@gmail.com>
parents: 10394
diff changeset
183 fb = fa;
9500a66118dc improve fzero
Jaroslav Hajek <highegg@gmail.com>
parents: 10394
diff changeset
184 elseif (fb == 0)
9500a66118dc improve fzero
Jaroslav Hajek <highegg@gmail.com>
parents: 10394
diff changeset
185 a = b;
9500a66118dc improve fzero
Jaroslav Hajek <highegg@gmail.com>
parents: 10394
diff changeset
186 fa = fb;
9500a66118dc improve fzero
Jaroslav Hajek <highegg@gmail.com>
parents: 10394
diff changeset
187 endif
9500a66118dc improve fzero
Jaroslav Hajek <highegg@gmail.com>
parents: 10394
diff changeset
188
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
189 itype = 1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
190
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
191 if (abs (fa) < abs (fb))
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
192 u = a; fu = fa;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
193 else
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
194 u = b; fu = fb;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
195 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
196
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
197 d = e = u;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
198 fd = fe = fu;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
199 mba = mu*(b - a);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
200 while (niter < maxiter && nfev < maxfev)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
201 switch (itype)
17174
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
202 case 1
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
203 ## The initial test.
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
204 if (b - a <= 2*(2 * abs (u) * eps + tolx))
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
205 x = u; fval = fu;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
206 info = 1;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
207 break;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
208 endif
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
209 if (abs (fa) <= 1e3*abs (fb) && abs (fb) <= 1e3*abs (fa))
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
210 ## Secant step.
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
211 c = u - (a - b) / (fa - fb) * fu;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
212 else
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
213 ## Bisection step.
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
214 c = 0.5*(a + b);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
215 endif
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
216 d = u; fd = fu;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
217 itype = 5;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
218 case {2, 3}
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
219 l = length (unique ([fa, fb, fd, fe]));
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
220 if (l == 4)
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
221 ## Inverse cubic interpolation.
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
222 q11 = (d - e) * fd / (fe - fd);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
223 q21 = (b - d) * fb / (fd - fb);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
224 q31 = (a - b) * fa / (fb - fa);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
225 d21 = (b - d) * fd / (fd - fb);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
226 d31 = (a - b) * fb / (fb - fa);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
227 q22 = (d21 - q11) * fb / (fe - fb);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
228 q32 = (d31 - q21) * fa / (fd - fa);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
229 d32 = (d31 - q21) * fd / (fd - fa);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
230 q33 = (d32 - q22) * fa / (fe - fa);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
231 c = a + q31 + q32 + q33;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
232 endif
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
233 if (l < 4 || sign (c - a) * sign (c - b) > 0)
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
234 ## Quadratic interpolation + newton.
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
235 a0 = fa;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
236 a1 = (fb - fa)/(b - a);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
237 a2 = ((fd - fb)/(d - b) - a1) / (d - a);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
238 ## Modification 1: this is simpler and does not seem to be worse.
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
239 c = a - a0/a1;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
240 if (a2 != 0)
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
241 c = a - a0/a1;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
242 for i = 1:itype
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
243 pc = a0 + (a1 + a2*(c - b))*(c - a);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
244 pdc = a1 + a2*(2*c - a - b);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
245 if (pdc == 0)
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
246 c = a - a0/a1;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
247 break;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
248 endif
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
249 c -= pc/pdc;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
250 endfor
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
251 endif
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
252 endif
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
253 itype += 1;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
254 case 4
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
255 ## Double secant step.
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
256 c = u - 2*(b - a)/(fb - fa)*fu;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
257 ## Bisect if too far.
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
258 if (abs (c - u) > 0.5*(b - a))
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
259 c = 0.5 * (b + a);
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
260 endif
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
261 itype = 5;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
262 case 5
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
263 ## Bisection step.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
264 c = 0.5 * (b + a);
17174
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 17097
diff changeset
265 itype = 2;
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
266 endswitch
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
267
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
268 ## Don't let c come too close to a or b.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
269 delta = 2*0.7*(2 * abs (u) * eps + tolx);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
270 if ((b - a) <= 2*delta)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
271 c = (a + b)/2;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
272 else
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
273 c = max (a + delta, min (b - delta, c));
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
274 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
275
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
276 ## Calculate new point.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
277 x = c;
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
278 fval = fc = fun (c);
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
279 niter ++; nfev ++;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
280
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
281 ## Modification 2: skip inverse cubic interpolation if
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
282 ## nonmonotonicity is detected.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
283 if (sign (fc - fa) * sign (fc - fb) >= 0)
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
284 ## The new point broke monotonicity.
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
285 ## Disable inverse cubic.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
286 fe = fc;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
287 else
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
288 e = d; fe = fd;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
289 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
290
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
291 ## Bracketing.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
292 if (sign (fa) * sign (fc) < 0)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
293 d = b; fd = fb;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
294 b = c; fb = fc;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
295 elseif (sign (fb) * sign (fc) < 0)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
296 d = a; fd = fa;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
297 a = c; fa = fc;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
298 elseif (fc == 0)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
299 a = b = c; fa = fb = fc;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
300 info = 1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
301 break;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
302 else
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
303 ## 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
304 error ("fzero:bracket", "fzero: zero point is not bracketed");
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
305 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
306
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
307 ## If there's an output function, use it now.
18132
8ca7b1906a41 fzero.m: Correctly detct if OutputFcn is present (bug #40889).
Rik <rik@octave.org>
parents: 17744
diff changeset
308 if (! isempty (outfcn))
10296
035ac548a67e implement fminbnd
Jaroslav Hajek <highegg@gmail.com>
parents: 9464
diff changeset
309 optv.funccount = nfev;
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
310 optv.fval = fval;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
311 optv.iteration = niter;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
312 if (outfcn (x, optv, "iter"))
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
313 info = -1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
314 break;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
315 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
316 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
317
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
318 if (abs (fa) < abs (fb))
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
319 u = a; fu = fa;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
320 else
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
321 u = b; fu = fb;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
322 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
323 if (b - a <= 2*(2 * abs (u) * eps + tolx))
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
324 info = 1;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
325 break;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
326 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
327
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
328 ## Skip bisection step if successful reduction.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
329 if (itype == 5 && (b - a) <= mba)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
330 itype = 2;
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
331 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
332 if (itype == 2)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
333 mba = mu * (b - a);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
334 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
335 endwhile
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
336
10337
7c97da90fc8f check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents: 10297
diff changeset
337 ## Check solution for a singularity by examining slope
7c97da90fc8f check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents: 10297
diff changeset
338 if (info == 1)
20038
9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19898
diff changeset
339 if ((b - a) != 0
9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19898
diff changeset
340 && abs ((fb - fa)/(b - a) / slope0) > max (1e6, 0.5/(eps+tolx)))
10337
7c97da90fc8f check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents: 10297
diff changeset
341 info = -5;
7c97da90fc8f check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents: 10297
diff changeset
342 endif
7c97da90fc8f check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents: 10297
diff changeset
343 endif
7c97da90fc8f check for discontinuities in fzero
Rik <rdrider0-list@yahoo.com>
parents: 10297
diff changeset
344
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
345 output.iterations = niter;
10296
035ac548a67e implement fminbnd
Jaroslav Hajek <highegg@gmail.com>
parents: 9464
diff changeset
346 output.funcCount = nfev;
10687
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
347 output.bracketx = [a, b];
a8ce6bdecce5 Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents: 10470
diff changeset
348 output.brackety = [fa, fb];
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
349
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
350 endfunction
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
351
8507
cadc73247d65 style fixes
John W. Eaton <jwe@octave.org>
parents: 8468
diff changeset
352 ## 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
353 ## bad results.
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
354 function fx = guarded_eval (fun, x)
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
355 fx = fun (x);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
356 fx = fx(1);
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
357 if (! isreal (fx))
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
358 error ("fzero:notreal", "fzero: non-real value encountered");
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
359 elseif (isnan (fx))
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
360 error ("fzero:isnan", "fzero: NaN value encountered");
8305
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
361 endif
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
362 endfunction
368b504777a8 implement fzero
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
363
17338
1c89599167a6 maint: End m-files with 1 blank line.
Rik <rik@octave.org>
parents: 17281
diff changeset
364
10296
035ac548a67e implement fminbnd
Jaroslav Hajek <highegg@gmail.com>
parents: 9464
diff changeset
365 %!shared opt0
035ac548a67e implement fminbnd
Jaroslav Hajek <highegg@gmail.com>
parents: 9464
diff changeset
366 %! opt0 = optimset ("tolx", 0);
14363
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14138
diff changeset
367 %!assert (fzero (@cos, [0, 3], opt0), pi/2, 10*eps)
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14138
diff changeset
368 %!assert (fzero (@(x) x^(1/3) - 1e-8, [0,1], opt0), 1e-24, 1e-22*eps)
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14138
diff changeset
369