rootfinding newtzero

[john]

Hello,

on 25 june 2012 I posted newtzero for finding real roots and imag. roots of
non-linear equations.

The problem was the syntax "fcnchk".Now i found on the internet the function
"fcnchk".It's matlab

but it works perfect in octave 3.2.4.

Yoy have to put the file "fcnchk" in the directory  \bin of octave.See below
for the programs:

function root = newtzero(f,xr,mx,tol)
%NEWTZERO finds roots of function using unique approach to Newton's Method.
% May find more than one root, even if guess is way off.  The function f 
% should be an inline object or function handle AND capable of taking 
% vector arguments.  The user can also pass in a string function. 
% NEWTZERO takes four arguments. The first argument is the function to be 
% evaluated, the second is an optional initial guess, the third is an 
% optional number of iterations, and the fourth is an absolute tolerance.  
% If the optional arguments are not supplied, the initial guess will be set 
% to 1, the number of iterations will be set to 30, and the tolerance will
% be set to 1e-13.  If the initial guess is a root, no iterations will take
% place.
%
% EXAMPLES:  
%
%         % Use any one of these three equivalent function definitions.
%         f = inline('cos(x+3).^3+(x-1).^2'); % Inline object.
%         f = 'cos(x+3).^3+(x-1).^2'; % String function.
%         f = @(x) cos(x+3).^3+(x-1).^2; % Anonymous function.
%         newtzero(f,900) % A very bad initial guess!        
%
%         fb = @(x) besselj(2,x)
%         rt = newtzero(fb); % Finds roots about the origin.
%         rt = rt(rt>=-100 & rt<=100); % Plot about origin.
%         x = -100:.1:100;
%         plot(x,real(fb(x)),rt,abs(fb(rt)),'*r')
%
%         f = @(x) x.^3 +1; 
%         newtzero(f,1) % Finds the real root;
%         newtzero(f,i) % Finds the real root and two imaginary roots.
%
%         f = @(x) x.^3 + 2*x.^2 + 3*x -exp(x)
%         newtzero(f)  % Finds two roots.
%
%         % Try it with Wilkinson's famous polynomial.
%         g = @(x)prod(bsxfun(@(x,y)(x-y),[1:20]',x.')).'; % For brevity
%         newtzero(g,0)
%         
%
% May work when the initial guess is outside the range of fzero, 
% for example compare:
%
%            fzero(f,900)   % for f as in the above example.
%
% This function may also find complex roots when fzero fails.  For example,
% try to find the roots of the following using NEWTZERO and FZERO:
%             
%            f1 = @(x) x.^(2*cos(x)) -x.^3 - sin(x)+1;            
%            ntz1 = newtzero(f1,[],[],1e-15) % NEWTZERO finds 5 roots
%            fzrt1 = fzero(f1,0) % FZERO aborts.
%
%            f2 = @(x) x.^2 + (2 + .1*i);  % could use 'roots' here.
%            ntz2 = newtzero(f2,1)  % NEWTZERO finds 2 roots
%            fzrt2 = fzero(f2,1) % FZERO fails.
% 
% See also fzero, roots
%
% Author:  Matt Fig
% Contact:  address@hidden

defaults = {1,30,1e-13};% Initial guess, max iterations, tolerance.

switch nargin  % Error checking and default assignments.
    case 1
        [xr,mx,tol] = defaults{:};
    case 2
        [mx,tol] = defaults{2:3};
    case 3
        tol = defaults{3};
end

if isempty(xr)   % User called newtzero(f,[],50,1e-3) for example.
    xr = defaults{1};
end

if isempty(mx)
    mx = 30;
end

if ~isa(xr,'double') 
    error('Only double values allowed.  See help examples.')
end

if tol < 0 || ~isreal(tol)
    error('Tolerance must be greater than zero.')
end

if mx < 1 || ~isreal(mx)
    error('Maximum number of iterations must be real and >0.')
end

[f,err] = fcnchk(f,'vectorized'); % If user passed in a string.

if ~isempty(err)
    error(['Error using NEWTZERO:',err.message])
end

if abs(f(xr))< tol 
    root = xr; % The user supplied a root as the initial guess.
    return  % Initial guess correct.
end  

LGS = logspace(0,3,220); % Can be altered for wider range or denser search.
LNS = 0:1/19:18/19;  % Search very close to initial guess too.
xr = [xr-LGS xr+LGS xr-LNS(2:end) xr+LNS].';  % Make vector of guesses.
iter = 1;  % Initialize the counter for the while loop.
mn1 = .1; % These will store the norms of the converging roots.
mn2 = 1; % See last comment.
sqrteps = sqrt(eps); % Used to make h.  See loop.
warning off MATLAB:divideByZero % WILL BE RESET AT THE END OF WHILE LOOP. 

while iter <= mx && abs(mn2-mn1) > 5*eps
    h = sqrteps*xr; % From numerical recipes, make h = h(xr)
    xr = xr-f(xr)./((f(xr+h)-f(xr-h))./(2*h)); % Newton's method.
    xr(isnan(xr) | isinf(xr)) = []; % No need for these anymore.
    mn1 = mn2; % Store the old value first.
    mn2 = norm(xr,'fro'); % This could make the loop terminate early!
    iter = iter+1;  % Increment the counter.
end

if abs(f(0)) < tol % The above method will tend to send zero root to Inf.
    xr = [xr;0];  % So explicitly check.
end

warning on MATLAB:divideByZero % Reset this guy, as promised.

% Filtering.  We want to filter out certain common results.
idxi = abs(imag(xr)) < 5e-15;  % A very small imag term is zero.
xr(idxi) = real(xr(idxi));  % Discard small imaginary terms.
idxr = abs(real(xr)) < 5e-15;  % A very small real term is zero.
xr(idxr) = complex(0,imag(xr(idxr))); % Discard small real terms.
root = xr(abs(f(xr)) < tol); % Apply the tolerance.

% Next we are going to delete repeat roots.  unique does not work in
% this case because many repeats are very close to each other but not
% equal.  For loops are fast enough here, most root vectors are short(ish).

if ~isempty(root)
    cnt = 1;  % Counter for while loop.
    
    while ~isempty(root)
        vct = abs(root - root(1))<5e-6; % Minimum spacing between roots.
        C = root(vct);  %C has roots grouped close together.
        [idx,idx] = min(abs(f(C)));  % Pick the best root per group.
        rt(cnt) = C(idx);  %#ok<AGROW>  Most root vectors are small.
        root(vct) = []; % Deplete the pool of roots.
        cnt = cnt + 1;  % Increment the counter.
    end
    
    root = sort(rt).';  % return a nice, sorted column vector
end

and here the function "fcnchk":    

function [f,msg] = fcnchk(fun,varargin) 
%FCNCHK Check FUNFUN function argument. 
%   FCNCHK(FUN,...) returns an inline object based on FUN if FUN 
%   is a string containing parentheses, variables, and math 
%   operators.  FCNCHK simply returns FUN if FUN is a function handle,  
%   or a matlab object with an feval method (such as an inline object).  
%   If FUN is a string name of a function (e.g. 'sin'), FCNCHK returns a 
%   function handle to that function. 
% 
%   FCNCHK is a helper function for FMINBND, FMINSEARCH, FZERO, etc. so they 
%   can compute with string expressions in addition to m-file functions. 
% 
%   FCNCHK(FUN,...,'vectorized') processes the string (e.g., replacing 
%   '*' with '.*') to produce a vectorized function. 
% 
%   When FUN contains an expression then FCNCHK(FUN,...) is the same as 
%   INLINE(FUN,...) except that the optional trailing argument 'vectorized' 
%   can be used to produce a vectorized function. 
% 
%   [F,MSG] = FCNCHK(...) returns an empty structure in MSG if successful 
%   or an error message structure if not. 
% 
%   See also ERROR, INLINE, @. 
 
%   Copyright 1984-2004 The MathWorks, Inc. 
%   $Revision: 1.29.4.7 $  $Date: 2004/04/16 22:05:23 $ 
 
message = ''; 
msgident = ''; 
 
nin = nargin; 
if (nin>1) && strcmp(varargin{end},'vectorized') 
    vectorizing = 1; 
    nin = nin-1; 
else 
    vectorizing = 0; 
end 
 
if ischar(fun) 
    fun = strtrim_local_function(fun); 
    % Check for non-alphanumeric characters that must be part of an 
    % expression. 
    if isempty(fun), 
        f = inline('[]'); 
    elseif ~vectorizing && isidentifier_local_function(fun) 
        f = str2func(fun); % Must be a function name only 
        % Note that we avoid collision of f = str2func(fun) with any local 
        % function named fun, by uglifying the local function's name 
 
        if isequal('x',fun) 
            warning('MATLAB:fcnchk:AmbiguousX', ... 
                ['Ambiguous expression or function input.\n The string ''x''
will be ',... 
                    'interpreted as the name of a ',... 
                'function called ''x'' \n (e.g., x.m) and not as the
mathematical expression ''x'' (i.e., f(x)=x). \n ',... 
                'Use inline(''x'') if you meant the mathematical expression
''x''.']); 
        end 
    else 
        if vectorizing 
            f = inline(vectorize(fun),varargin{1:nin-1}); 
            var = argnames(f); 
            f = inline([formula(f) '.*ones(size(' var{1} '))'],var{1:end}); 
        else 
            f = inline(fun,varargin{1:nin-1}); 
        end  
    end 
elseif isa(fun,'function_handle')  
    f = fun;  
    % is it a matlab object with a feval method? 
elseif isobject(fun) 
    % delay the methods call unless we know it is an object to avoid runtime
error for compiler 
    meths = methods(class(fun)); 
    if any(strmatch('feval',meths,'exact')) 
       if vectorizing && any(strmatch('vectorize',meths,'exact')) 
          f = vectorize(fun); 
       else 
          f = fun; 
       end 
    else % no feval method 
        f = ''; 
        message = 'If FUN is a MATLAB object, it must have an feval
method.'; 
        msgident = 'MATLAB:fcnchk:objectMissingFevalMethod'; 
    end 
else 
    f = ''; 
    message = ['FUN must be a function, a valid string expression, ', ... 
            sprintf('\n'),'or an inline function object.']; 
    msgident = 'MATLAB:fcnchk:invalidFunctionSpecifier'; 
end 
 
% If no errors and nothing to report then we are done. 
if nargout < 2 && isempty(message) 
    return 
end 
 
% compute MSG 
if isempty(message) 
    msg.message = ''; 
    msg.identifier = ''; 
    msg = msg(zeros(0,1)); % make sure msg is the right dimension 
else 
    msg.message = message; 
    msg.identifier = msgident; 
end 
 
if nargout < 2, error(msg); end 
 
 
%------------------------------------------ 
function s1 = strtrim_local_function(s) 
%STRTRIM_LOCAL_FUNCTION Trim spaces from string. 
% Note that we avoid collision with line 45: f = str2func('strtrim') 
% by uglifying the local function's name 
 
if isempty(s) 
    s1 = s; 
else 
    % remove leading and trailing blanks (including nulls) 
    c = find(s ~= ' ' & s ~= 0); 
    s1 = s(min(c):max(c)); 
end 
 
%------------------------------------------- 
function tf = isidentifier_local_function(str) 
% Note that we avoid collision with line 45: f = str2func('isidentifier') 
% by uglifying the local function's name 
 
 
tf = false; 
 
if ~isempty(str) 
    first = str(1); 
    if (isletter(first)) 
        letters = isletter(str); 
        numerals = (48 <= str) & (str <= 57); 
        underscore = (95 == str); 
        tf = all(letters | numerals | underscore); 
    end 
end 
 


when I run the program I get the same result as in the example

octave-3.2.4.exe:1> f1=@(x) x.^(2*cos(x))-x.^3-sin(x)+1;

octave-3.2.4.exe:3> rt=newtzero(f1,[],[],1e-15)
rt =

   0.41924 - 0.72287i
   0.41924 + 0.72287i
   1.05841 + 0.00000i
  -0.42211 - 1.26440i
  -0.42211 + 1.26440i
as you see the function "fcnchk" does the job.

if you have suggestions ,they are welcome

Cheers

John




--
View this message in context: 
http://octave.1599824.n4.nabble.com/rootfinding-newtzero-tp4631423.html
Sent from the Octave - General mailing list archive at Nabble.com.