Menu
▾
▴
[OctDev] fzero.m doesn't pass varargin to feval
|
From: Ben B. <bar...@al...> - 2005-05-10 18:45:16
|
fzero doesn't seem to pass varargin to feval. There are 6 places=20
including the function statement itself where varargin (function=20
statement) or varargin{:} (when using feval) is needed. I have tested=20
this modified fzero.m for varargin functionality and include it below.
## Copyright (C) 2004 =A3ukasz Bodzon, <lll...@o2...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 US=
A
##
## REVISION HISTORY
##
## 2004-07-20, Piotr Krzyzanowski, <pio...@mi...>:
## Options parameter and fall back to fsolve if only scalar APPROX argume=
nt
## supplied
##
## 2004-07-01, Lukasz Bodzon:
## Replaced f(a)*f(b) < 0 criterion by a more robust
## sign(f(a)) ~=3D sign(f(b))
##
## 2004-06-18, Lukasz Bodzon:
## Original implementation of Brent's method of finding a zero of a scala=
r
## function
## -*- texinfo -*-
## @deftypefn {Function File} {} [X, FX, INFO] =3D fzero (FCN, APPROX,=20
OPTIONS)
##
## Given FCN, the name of a function of the form `F (X)', and an initial
## approximation APPROX, `fzero' solves the scalar nonlinear equation=20
such that
## `F(X) =3D=3D 0'. Depending on APPROX, `fzero' uses different algorithm=
s=20
to solve
## the problem: either the Brent's method or the Powell's method of=20
`fsolve'.
##
## @table @asis
## @item INPUT ARGUMENTS
## @end table
##
## @table @asis
## @item APPROX can be a vector with two components,
## @example
## A =3D APPROX(1) and B =3D APPROX(2),
## @end example
## which localizes the zero of F, that is, it is assumed that X lies=20
between A and
## B. If APPROX is a scalar, it is treated as an initial guess for X.
##
## If APPROX is a vector of length 2 and F takes different signs at A and=
B,
## F(A)*F(B) < 0, then the Brent's zero finding algorithm [1] is used=20
with error
## tolerance criterion
## @example
## reltol*|X|+abstol (see OPTIONS).
## @end example
## This algorithm combines
## superlinear convergence (for sufficiently regular functions) with the
## robustness of bisection.
##
## Whether F has identical signs at A and B, or APPROX is a single=20
scalar value,
## then `fzero' falls back to another method and `fsolve(FCN, X0)' is=20
called, with
## the starting value X0 equal to (A+B)/2 or APPROX, respectively. Only=20
absolute
## residual tolerance, abstol, is used then, due to the limitations of=20
the `fsolve_options'
## function. See OPTIONS and `help fsolve' for details.
##
## @item OPTIONS is a structure, with the following fields:
##
## @table @asis
## @item 'abstol' - absolute (error for Brent's or residual for fsolve)
## tolerance. Default =3D 1e-6.
##
## @item 'reltol' - relative error tolerance (only Brent's method).=20
Default =3D 1e-6.
##
## @item 'prl' - print level, how much diagnostics to print. Default =3D =
0, no
## diagnostics output.
## @end table
##
## If OPTIONS argument is omitted, or a specific field is not present in =
the
## OPTIONS structure, default values will be used.
## @end table
##
## @table @asis
## @item OUTPUT ARGUMENTS
## @end table
##
## @table @asis
## @item The computed approximation to the zero of FCN is returned in X.=20
FX is then equal
## to FCN(X). If the iteration converged, INFO =3D=3D 1. If Brent's metho=
d=20
is used,
## and the function seems discontinuous, INFO is set to -5. If fsolve is=20
used,
## INFO is determined by its convergence.
## @end table
##
## @table @asis
## @item EXAMPLES
## @end table
##
## @example
## fzero('sin',[-2 1]) will use Brent's method to find the solution to
## sin(x) =3D 0 in the interval [-2, 1]
## @end example
##
## @example
## [x, fx, info] =3D fzero('sin',-2) will use fsolve to find a solution t=
o
## sin(x)=3D0 near -2.
## @end example
##
## @example
## options.abstol =3D 1e-2; fzero('sin',-2, options) will use fsolve to
## find a solution to sin(x)=3D0 near -2 with the absolute tolerance 1e-2.
## @end example
##
## @table @asis
## @item REFERENCES
## [1] Brent, R. P. "Algorithms for minimization without derivatives"=20
(1971).
## @end table
## @end deftypefn
## @seealso{fsolve}
function [Z, FZ, INFO] =3Dfzero(Func,bracket,options,varargin)
if (nargin < 2)
usage("[x, fx, info] =3D fzero(@fcn, [lo,hi]|start, options)");
endif
if !isstr(Func) && !isa(Func,"function handle") && !isa(Func,"inline=20
function")
error("fzero expects a function as the first argument");
endif
bracket =3D bracket(:);
if all(length(bracket)!=3D[1,2])
error("fzero expects an initial value or a range");
endif
set_default_options =3D false;
if (nargin >=3D 2) % check for the options
if (nargin =3D=3D 2)
set_default_options =3D true;
options =3D [];
else % nargin > 2
if ~isstruct(options)
warning('Options incorrect. Setting default values.');
set_default_options =3D true;
end
end
end
if ~isfield(options,'abstol')
options.abstol =3D 1e-6;
end
if ~isfield(options,'reltol')
options.reltol =3D 1e-6;
end
% if ~isfield(options,'maxit')
% options.maxit =3D 100;
% end
if ~isfield(options,'prl')
options.prl =3D 0; % no diagnostics output
end
fcount =3D 0; % counts function evaluations
if (length(bracket) > 1)
a =3D bracket(1); b =3D bracket(2);
use_brent =3D true;
else
b =3D bracket;
use_brent =3D false;
end
if (use_brent)
fa=3Dfeval(Func,a,varargin{:}); fcount=3Dfcount+1;
fb=3Dfeval(Func,b,varargin{:}); fcount=3Dfcount+1;
BOO=3Dtrue;
tol=3Doptions.reltol*abs(b)+options.abstol;
% check if one of the endpoints is the solution
if (fa =3D=3D 0.0)
BOO =3D false;
c =3D b =3D a;
fc =3D fb =3D fa;
end
if (fb =3D=3D 0.0)
BOO =3D false;
c =3D a =3D b;
fc =3D fa =3D fb;
end
if ((sign(fa) =3D=3D sign(fb)) & BOO)
warning ("fzero: equal signs at both ends of the interval.\n\
Using fsolve('%s',%g) instead", Func, 0.5*(a+b));
use_brent =3D false;
b =3D 0.5*(a+b);
endif
end
if (use_brent) % it is reasonable to call Brent's met=
hod
if options.prl > 0
fprintf(stderr,"=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D\n");
fprintf(stderr,"fzero: using Brent's method\n");
fprintf(stderr,"=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D\n");
end
c=3Da;
fc=3Dfa;
d=3Db-a;
e=3Dd;
while (BOO =3D=3D true) % convergence check
if (sign(fb) =3D=3D sign(fc)) % rename a, b, c and adjust=20
bounding interval
c=3Da;
fc=3Dfa;
d=3Db-a;
e=3Dd;
endif,
## We are preventing overflow and division by zero
## while computing the new approximation by
## linear interpolation.
## After this step, we lose the chance for using
## inverse quadratic interpolation (a=3D=3Dc).
if (abs(fc) < abs(fb))
a=3Db;
b=3Dc;
c=3Da;
fa=3Dfb;
fb=3Dfc;
fc=3Dfa;
endif,
tol=3Doptions.reltol*abs(b)+options.abstol;
m=3D0.5*(c-b);
if options.prl > 0
fprintf(stderr,'fzero: [%d feval] X =3D %8.4e\n', fcount,=
b );
if options.prl > 1
fprintf(stderr,'fzero: m =3D %8.4e e =3D %8.4e [tol =3D=
=20
%8.4e]\n', m, e, tol);
end
end
if (abs(m) > tol & fb !=3D 0)
## The second condition in following if-instruction
## prevents overflow and division by zero
## while computing the new approximation by
## inverse quadratic interpolation.
if (abs(e) < tol | abs(fa) <=3D abs(fb))
d=3Dm; % bisection
e=3Dm;
else
s=3Dfb/fa;
if (a =3D=3D c) % attempt linear interpolatio=
n
p=3D2*m*s; % (the secant method)
q=3D1-s;
else % attempt inverse quadratic=20
interpolation
q=3Dfa/fc;
r=3Dfb/fc;
p=3Ds*(2*m*q*(q-r)-(b-a)*(r-1));
q=3D(q-1)*(r-1)*(s-1);
endif,
if (p > 0) % fit signs
q=3D-q; % to the sign of (c-b)
else
p=3D-p;
endif,
s=3De;
e=3Dd;
if (2*p < 3*m*q-abs(tol*q) & p < abs(0.5*s*q))
d=3Dp/q; % accept interpolation
else % interpolation failed;
d=3Dm; % take the bisection step
e=3Dm;
endif,
endif,
a=3Db;
fa=3Dfb;
if (abs(d) > tol) % the step we take is never=20
shorter
b=3Db+d; % than tol
else
if (m > 0) % fit signs
b=3Db+tol; % to the sign of (c-b)
else
b=3Db-tol;
endif,
endif,
fb=3Dfeval(Func,b,varargin{:}); fcount=3Dfcount+1;
else
BOO=3Dfalse;
endif,
endwhile,
Z=3Db;
FZ =3D fb;
if abs(FZ) > 100*tol % large value of the residual may=20
indicate a discontinuity point
INFO =3D -5;
else
INFO =3D 1;
end
%
% TODO: test if Z may be a singular point of F (ie F is=20
discontinuous at Z
% Then return INFO =3D -5
%
if (options.prl > 0 )
fprintf(stderr,"\nfzero: summary\n");
switch(INFO)
case 1
MSG =3D "Solution converged within specified tolerance";
case -5
MSG =3D strcat("Probably a discontinuity/singularity poin=
t=20
of F()\n encountered close to X =3D ", sprintf('%8.4e',Z),...
".\n Value of the residual at X, |F(X)| =3D ",...
sprintf('%8.4e',abs(FZ)), ...
".\n Another possibility is that you use too large=20
tolerance parameters",...
".\n Currently TOL =3D ", sprintf('%8.4e', tol), ...
".\n Try fzero with smaller tolerance values");
otherwise
MSG =3D "Something strange happened"
endswitch
fprintf(stderr,' %s.\n', MSG);
fprintf(stderr,' %d function evaluations.\n', fcount);
end
else % fall back to fsolve
if options.prl > 0
fprintf(stderr,"=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D\n");
fprintf(stderr,"fzero: using fsolve\n");
fprintf(stderr,"=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D\n");
end
% check for zeros in APPROX
fb=3Dfeval(Func,b,varargin{:});
fcount=3Dfcount+1;
tol_save =3D fsolve_options('tolerance');
fsolve_options("tolerance",options.abstol);
[Z, INFO, MSG] =3D fsolve(Func, b);
fsolve_options('tolerance',tol_save);
FZ =3D feval(Func,Z,varargin{:});
if options.prl > 0
fprintf(stderr,"\nfzero: summary\n");
fprintf(stderr,' %s.\n', MSG);
end
end
endfunction;
%!## usage and error testing
%!## the Brent's method
%!test
%! options.abstol=3D0;
%! assert (fzero('sin',[-1,2],options), 0)
%!test
%! options.abstol=3D0.01;
%! options.reltol=3D1e-3;
%! assert (abs(fzero('tan',[-0.5,1.41],options)), 0, 0.01)
%!test
%! options.abstol=3D1e-3;
%! assert (abs(fzero('atan',[-(10^300),10^290],options)), 0, 1e-3)
%!test
%! testfun=3Dinline('(x-1)^3','x');
%! options.abstol=3D0;
%! options.reltol=3Deps;
%! assert (abs(fzero(testfun,[0,3],options)), 1, -eps)
%!test
%! testfun=3Dinline('x.^2-100','x');
%! options.abstol=3D1e-4;
%! assert (abs(fzero(testfun,[-9,300],options)),10,1e-4)
%!## `fsolve'
%!test
%! options.abstol=3D0.01;
%! assert (abs(fzero('tan',-0.5,options)), 0, 0.01)
%!test
%! options.abstol=3D0;
%! assert (fzero('sin',[0.5,1],options), 0)
%!
%!demo
%! bracket=3D[-1,1.2];
%! [X,FX,MSG]=3Dfzero('tan',bracket)
%!demo
%! bracket=3D1; # `fsolve' will be used
%! [X,FX,MSG]=3Dfzero('sin',bracket)
%!demo
%! bracket=3D[-1,2];
%! options.abstol=3D0; options.prl=3D1;
%! X=3Dfzero('sin',bracket,options)
%!demo
%! bracket=3D[0.5,1];
%! options.abstol=3D0; options.reltol=3Deps; options.prl=3D1;
%! fzero('sin',bracket,options)
%!demo
%! demofun=3Dinline('2*x.*exp(-4)+1 - 2*exp(-4*x)','x');
%! bracket=3D[0, 1];
%! options.abstol=3D1e-14; options.reltol=3Deps; options.prl=3D2;
%! [X,FX]=3Dfzero(demofun,bracket,options)
%!demo
%! demofun=3Dinline('x^51','x');
%! bracket=3D[-12,10];
%! # too large tolerance parameters
%! options.abstol=3D1; options.reltol=3D1; options.prl=3D1;
%! [X,FX]=3Dfzero(demofun,bracket,options)
%!demo
%! # points of discontinuity inside the bracket
%! demofun=3Dinline('0.5*(sign(x-1e-7)+sign(x+1e-7))','x');
%! bracket=3D[-5,7];
%! options.prl=3D1;
%! [X,FX]=3Dfzero(demofun,bracket,options)
%!demo
%! demofun=3Dinline('2*x*exp(-x^2)','x');
%! bracket=3D1;
%! options.abstol=3D1e-14; options.prl=3D2;
%! [X,FX]=3Dfzero(demofun,bracket,options)
%!demo
%! demofun=3Dinline('2*x.*exp(-x.^2)','x');
%! bracket=3D[-10,1];
%! options.abstol=3D1e-14; options.prl=3D2;
%! [X,FX]=3Dfzero(demofun,bracket,options)
--=20
-------------------- Benjamin E. Barrowes -------------------
Los Alamos National Laboratory bar...@al...
Biophysics Group P-21, MS-D454 Phone:(505)606-0105
Los Alamos, NM 87544 FAX:(270)294-1268
-------------------------------------------------------------
|
