roots bessel functions(2) answer to Jordi

[john]

Hi Jordi,

I downloaded "findzeros from  a matlab site on the internet.It works perfect
in octave.

Here is the program:
%----------------------------------------------------------------------
% Usage: [Z, err, fZ] = findzeros(x_start, x_step, x_end, fun)
% Find roots of f(x) = 0 in the interval [x_start, x_end].
%----------------------------------------------------------------------
%   fun:  Function handle (e.g. @sin) or the function's name as string.
% Returns:
%    Z :   The zeros.
%  err :   Values abs(a-b), where a, b are the closest bounds such that
%              f(a) <= 0 <= f(b) or
%              f(a) >= 0 >= f(b) and abs(a-b) <= eps.
%   fZ :   Values f(Z)
%-----------------------------------------------------------------------
% The function f(x) is sampled at equally spaced x's of distance x_step,
% then bisection is used to find the zeros of f(x).
% In order to find _all_ zeros, choose x_step small enough, so that:
%  1. There is at most one zero between two samples.
%  2. For every zero x_0 there are samples f(a), f(b) with
%        a < x_0 < b and
%        sign(f(a)) != sign(f(b)).
% Otherwise only some zeros will be found. If all samples are > 0 or < 0
% respectively, the procedure will fail. Remark: Zeros, that are local
% minima/maxima might be "overlooked". This procedure is especially
% _not_ suitable when all zeros are local minima/maxima.
%-----------------------------------------------------------------------
% Example 1: Finding zeros of sin(x) in [0,10].
%    Z = findzeros(0, 0.5*pi(), 10, "sin")
%-----------------------------------------------------------------------
% Example 2: Finding zeros of bessel functions of the first kind. (The
% procedure is suitable for functions with a single argument only, so
% we need helper functions.)
%    function y = J0(x); y = besselj(0,x); endfunction;
%    function y = J1(x); y = besselj(1,x); endfunction;
%    function y = J2(x); y = besselj(2,x); endfunction;
%    [Z0,er0] = findzeros(0, 0.25*pi(), 100, "J0");
%    [Z1,er1] = findzeros(0, 0.25*pi(), 100, "J1");
%    [Z2,er2] = findzeros(0, 0.25*pi(), 100, "J2");
%-----------------------------------------------------------------------
% Warnings: This might also find "zeros" that in reality are poles:
%   > function y = f(x); y = 1 ./ x; end
%   >[z,r,fZ] = findzeros(-1,0.3,1,"f")
%   warning: division by zero
%   z = -4.9407e-324
%   r =  4.9407e-324
%   fZ = -Inf
% Also, precisely "hitting" the zeros with the samples will lead to
% strange results / overlooking of zeros. The following zeros of sin(x)
% are only found due to roundoff errors.
%   > [z,r] = findzeros(0,2*pi(),100,"sin")
%   z =
%       0.00000
%      81.68141
%      84.82300
%   r =
%      4.9407e-324
%      1.4211e-014
%      1.4211e-014
%-----------------------------------------------------------------------
% Tested with Octave 3.2.4; public domain; no warranties whatsoever; not
% built for speed; don't put your life on it; don't sue me :)
% Version: 2010-04-14, S. Reiser (address@hidden).
%-----------------------------------------------------------------------
function [zeros, err, fZ] = findzeros(xA, step, xE, fname)
   numZeros = 0;             % Anzahl bisher gefundener Nullstellen.
   A = NaN;                  % Matrix mit Intervallen um Nullstellen.
   X = (xA:step:xE)';        % Intervallgrenzen festlegen.
   if X(rows(X) < xE)        % Intervall endet mit xE!
      X(rows(X)+1) = xE;
   endif
   Y = arrayfun(fname, X);   % Funktionswerte bestimmen.
   i = 1;                    % Vorzeichenwechsel finden.
   LL = rows(X);
   SGN = sign(Y(1));
   while (i <= LL)
      if (SGN != sign(Y(i))) % Vorzeichenwechsel! Nur für i > 1 möglich.
         SGN = sign(Y(i));
         A(++numZeros,1) = X(i-1);
         A(  numZeros,2) = X(i);
      endif
      i++;
   endwhile
   if (numZeros < 1)         % Keine Vorzeichenwechsel gefunden.
      error "No zeros found. Try smaller x_step."
   endif
   % Matrix A enthält nun je Zeile ein Intervall, in dem sich zumindest
   % eine Nullstelle befindet. Wir finden je Intervall eine Nullstelle
   % per Bisektion.
   i = 1;
   LL = rows(A);
   while (i <= LL)
      a = A(i,1);  b = A(i,2);
      F(i) = feval(fname, A(i,1));
      F_b = feval(fname, A(i,2));
      while true
         c = (a + b) / 2;
         F_c = feval(fname, c);
         if ((c == a) || (c == b))
            break;
         elseif (F_c == 0)
            a = c;  b = c;
            break;
         elseif (sign(F_c) == sign(F(i)))
            a = c;
         else
            b = c;
         endif
      endwhile
      A(i,1) = a;  A(i,2) = b;
      F(i) = F_c;
      i++;
   endwhile
   zeros = A(:,1);
   err = (A(:,2) - A(:,1));
   fZ  = F;
endfunction

example

function y=Dini(x)
y=@(x) x*besselj(1,x)+3*besselj(0,x);
endfunction;




octave-3.2.4.exe:2> findzeros(0,0.25*pi,100,Dini)
ans =

    3.1313
    6.6013
    9.8829
   13.1004
   16.2895
   19.4635
   22.6287
   25.7881
   28.9437
   32.0966
   35.2475
   38.3969
   41.5451
   44.6923
   47.8389
   50.9848
   54.1302
   57.2752
   60.4199
   63.5642
   66.7083
   69.8522
   72.9958
   76.1393
   79.2827
   82.4259
   85.5690
   88.7120
   91.8549
   94.9977
   98.1404
My question was just is there a tool in octave to compute the complex zeros
of nonlinear equation.

Thank you

John

--
View this message in context: 
http://octave.1599824.n4.nabble.com/roots-bessel-functions-2-answer-to-Jordi-tp4620176.html
Sent from the Octave - General mailing list archive at Nabble.com.