Positive Roots of bessel functions.

i have this type of potential and i want to find the positive roots of the bessel functions(also for the derivative in denomerator)

i want to find out this wave number. How i can find it?

[David Goodmanson]

Hi Javeria, please see the comment posted after my answer.

Answer 1

[Torsten]

Making a Google search with

file exchange & zeros of bessel functions & matlab

will give you some codes you can test for your purpose, e.g.

https://uk.mathworks.com/matlabcentral/fileexchange/48403-bessel-zero-solver

https://uk.mathworks.com/matlabcentral/fileexchange/28001-bessel-derivative-zeros?s_tid=prof_contriblnk

[Javeria]

I search these code but i did not understand how to make use of it for my purpose. Like i want this code for where . Can someone please edit the above codes for this value.

[Walter Roberson]

You are not going to find code that finds the zeros of the bessel function analytically, only code that finds the zeros by doing some variety of zero finder starting from approximate solutions.

Answer 2

[David Goodmanson]

Hi Javeria, here is a function besse0j to calculate the roots of J and J'. (at the origin, Jn(0) = 0 for n>0 but these zeros are not shown).

examples:

% first seven roots of J1
r = bessel0j(1,7)
    r = 3.8317    7.0156   10.1735   13.3237   16.4706   19.6159   22.7601

% check J1 at the roots, should be tiny
chk = besselj(1,r)
ans = 1.0e-15 *
      0.0730    0.1771    0.2805    0.1854   -0.1615    0.0675   -0.0833

% first seven roots of J2'
s =  bessel0j(2,7,'d')
    s = 3.0542    6.7061    9.9695   13.1704   16.3475   19.5129   22.6716

The function can supply a large number of roots and it's fast.

function x = bessel0j(n,q,opt)
% a row vector of the first q roots of bessel function Jn(x), integer order.
% if opt = 'd', row vector of the first q roots of dJn(x)/dx, integer order.
% if opt is not provided, the default is zeros of Jn(x).
% all roots are positive, except when n=0,
% x=0 is included as a root of dJ0(x)/dx (standard convention).
%
% starting point for for zeros of Jn was borrowed from Cleve Moler, 
% but the starting points for both Jn and Jn' can be found in
% Abramowitz and Stegun 9.5.12, 9.5.13. 
%
% David Goodmanson
%
% x = bessel0j(n,q,opt)

k = 1:q;
if nargin==3 & opt=='d'
    beta = (k + n/2 - 3/4)*pi;
    mu = 4*n^2;
    x = beta - (mu+3)./(8*beta) - 4*(7*mu^2+82*mu-9)./(3*(8*beta).^3);
    for j=1:8
        xnew = x - besseljd(n,x)./ ...
            (besselj(n,x).*((n^2./x.^2)-1) -besseljd(n,x)./x);
        x = xnew;    
    end
    if n==0
        x(1) = 0;    % correct a small numerical difference from 0
    end
else
    beta = (k + n/2 - 1/4)*pi;
    mu = 4*n^2;
    x = beta - (mu-1)./(8*beta) - 4*(mu-1)*(7*mu-31)./(3*(8*beta).^3);
    for j=1:8
        xnew = x - besselj(n,x)./besseljd(n,x);
        x = xnew;
    end
end

[Javeria]

My question is that all the above codes will give x. As in my case i have

So i want to modify the above codes for this. I want that the code will give me for N=20(i can change this value according to my desire).

Also then i want to calculate , so can i used the above one for this one case also. because the for the derivatives of bessel functions we have some recurrence formulas.

I hope someone solved my this one problem.

Thanks

[David Goodmanson]

Hi Javeria, if you know R then by middle of the second line that you have above,

chi = x/R

So for a given L (this font does a terrible job of showing lowercase L) the first twenty wave numbers are

chi = bessel0j(L,20)/R

and the values of J', evaluated at the roots of J are

<MODIFIED to correct an error>

Jprime = besseljd(L,chi*R)
% or you could just use
Jprime = besseljd(L,x)
% if x is still around.

but you can't do anything unless you know R.

function y = besseljd(n,x);
% derivative of bessel function of integer order
% 
% function y = besseljd(n,x);

y = -besselj(n+1,x) + n*besselj(n,x)./x;
% get rid of nans at the origin
if n==1
    y(x==0) = 1/2;
else
    y(x==0) = 0;
end

Answer 3

[Torsten]

As I understand your question, given l, you try to find 0 < x_1 < x_2 < x_3 < ... < x_m that makes

J_l(x_i) = 0 for i = 1,2,...,m. That's what the codes I linked to do.

The i-th wavenumber is then x_i/R.

For your specific question:

format long
m = 20;                    % maximum number of roots required
kind = 1;                  % We are dealing with bessel functions of the 1st kind

l0 = 0;                    % order of bessel function = 0
x = besselzero(l0,m,kind); % compute the m positive roots j_{0,1:20}
besselj(0,x)               % Check if J0(x) = 0

ans = 1×20

1.0e-15 * 0 0.072280293343350 0.106532214258324 0.317198790301651 0.126781573163115 0.210495026175707 0.230136870533624 -0.148300789919755 0.246413959240029 0.077521514175126 0.199635584581609 0.121787667715647 -0.033787009113403 -0.103897431962289 -0.106388428011883 -0.270228885935885 -0.193112818680809 0.240140662229920 0.195700929034410 -0.359179341848299

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

x0 = x(m)                  % take the m'th root = j_{0,20}

x0 = 
  62.048469190227166

l1 = 1;                    % order of bessel function = 1
x = besselzero(l1,m,kind); % compute the m positive roots j_{1,1:20}
besselj(1,x)               % Check if J1(x) = 0

ans = 1×20

1.0e-15 * 0.072970654098185 0.177126737870703 -0.405185458628602 0.185384865640674 -0.161462937561244 0.067498797579601 -0.083254253021943 -0.076410056217517 -0.040830239270359 0.217654079300945 -0.437666769753990 -0.092934445411549 0.070048687752665 -0.277836538178080 -0.399356446774912 -0.302261212094386 0.238746153819269 -0.152502037319241 -0.024295541444588 -0.176575262730997

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

x1 = x(m)                 % take the m'th root = j_{1,20}

x1 = 
  63.611356698481231

dJ0dx0 = dbesselj(l0,x0)  % compute derivative J_0'(j_{0,20})

dJ0dx0 = 
   0.101293498933943

(besselj(l0,x0+1e-8)-besselj(l0,x0-1e-8))/2e-8  % check derivative

ans = 
   0.101293507314989

dJ1dx1 = dbesselj(l1,x1)  % compute derivative J_1'(j_{1,20})

dJ1dx1 = 
   0.100035146811523

(besselj(l1,x1+1e-8)-besselj(l1,x1-1e-8))/2e-8  % check derivative

ans = 
   0.100035155087701

function x = besselzero(n, k, kind)
% besselzero calculates the zeros of Bessel function of the first and second kind
%
%   x = besselzero(n)
%   x = besselzero(n, k)
%   x = besselzero(n, k, kind)
%
%% Inputs
% * *n* - The order of the bessel function. n can be a scalar, vector, or
%       matrix.  n can be positive, negative, fractional, or any
%       combinaiton. abs(n) must be less than or equal to
%       146222.16674537213 or 370030.762407380 for first and second kind
%       respectively. Above these values, this algorithm will not find the
%       correct zeros because of the starting values therefore an error is
%       thrown instead.
% * k - The number of postive zeros to calculate.  When k is not supplied,
%       k = 5 is the default. k must be a scalar.
% * kind - kind is either 1 or 2. When kind is not supplied, default is
%          kind = 1.
%
%% Outputs
% * x - The calculated zeros.  size(x) = [size(n) k].
%
%% Description
% besselzero calculates the first k positive zeros of nth order bessel
% function of the first or second kind.  Note, that zero is not included as
% the first zero.
%
%% Algorithm
% the first three roots of any order bessel can be approximated by a simple
% equations.  These equations were generated using a least squares fit of
% the roots from orders of n=0:10000. The approximation is used to start
% the iteration of Halley's method.  The 4th and higher roots can be
% approximated by understanding the roots are regularly spaced for a given
% order.  Once the 2nd and 3rd roots are found, the spacing can be
% approximated by the distance between the 2nd and 3rd root.  Then again
% Halley's method can be applied to precisely locate the root.
%%
% Because the algorithm depends on good guesses of the first three zeros,
% if the guess is to far away then Halley's method will converge to the
% wrong zero which will subsequently cause any other zero to be incorrectly
% located. Therefore, a limit is put on abs(n) of 146222.16674537213 and
% 370030.762407380 for first and second kind respectively.  If n is
% specified above these limits, then an error is thrown.
%
%% Example
%   n = (1:2)';
%   k = 10;
%   kind = 1;
%   z = besselzero(n, k, kind);
%   x = linspace(0, z(end), 1000);
%   y = nan(2, length(x));
%   y(1,:) = besselj(n(1), x);
%   y(2,:) = besselj(n(2), x);
%   nz = nan(size(z));
%   nz(1,:) = besselj(n(1), z(1,:));
%   nz(2,:) = besselj(n(2), z(2,:));
%   plot(x, y, z, nz,'kx')
%

% Originally written by
% Written by: Greg von Winckel - 01/25/05
% Contact: gregvw(at)chtm(dot)unm(dot)edu
%
% Modified, Improved, and Documented by
% Jason Nicholson 2014-Nov-06
% Contact: [e-mail]

%% Change Log
% * Original release. 2005-Jan-25, Greg von Winckel.
% * Updated Documentation and commented algorithm. Fixed bug in finding the
%   the first zero of the bessel function of the second kind. Improved speed by
%   factor of 20. 2014-Nov-06, Jason Nicholson.
%

%% Input checking
assert(nargin >=1 | nargin <=3,'Wrong number of input arguments.');

% Take care of default cases of k and kind
if nargin < 2
    k = 5;
end

if nargin < 3
    kind = 1;
end

assert(isscalar(kind) & any(kind == [1 2]), '''kind''must be a scalar with value 1 or 2 only.');
assert(isscalar(k) & fix(k)==k & k>0, 'k must a positive scalar integer.');
assert(all(isreal(n(:))), 'n must be a real number.');

% negative orders have the same roots as the positive orders
n = abs(n);

% Check for that n is less than the ORDER_MAX
if kind==1
    ORDER_MAX = 146222.16674537213;
    assert(all(n <= ORDER_MAX), 'all n values must be less than or equal %6.10f for kind=1.', ORDER_MAX);
elseif kind==2
    ORDER_MAX = 370030.762407380;
    assert(all(n(:) <= ORDER_MAX), 'all n values must be less than or equal %6.10f for kind=2.', ORDER_MAX);
end
%% Setup Arrays

% output size
nSize = size(n);
if nSize(end) ==1
    outputSize = [nSize(1:end-1) k];
else
    outputSize = [nSize k];
end
% number of orders for each kth root
nOrdersPerRoot = prod(outputSize(1:end-1));
x = nan(outputSize);

%% Solve for Roots
switch kind
    case 1
        % coefficients and exponent are from least squares fitting the k=1,
        % n=0:10000.
        coefficients1j = [0.411557013144507;0.999986723293410;0.698028985524484;1.06977507291468];
        exponent1j = [0.335300369843979,0.339671493811664];
        % guess for k = 1
        x((1:nOrdersPerRoot)') = coefficients1j(1) + coefficients1j(2)*n(:) + coefficients1j(3)*(n(:)+1).^(exponent1j(1)) + coefficients1j(4)*(n(:)+1).^(exponent1j(2));
        % find first zero
        x((1:nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 1, x0, kind), n(:), x((1:nOrdersPerRoot)'));

        if k >= 2
            % coefficients and exponent are from least squares fitting the k=2,
            % n=0:10000.
            coefficients2j = [1.93395115137444;1.00007656297072;-0.805720018377132;3.38764629174694];
            exponent2j = [0.456215294517928,0.388380341189200];
            % guess for k = 2
            x((nOrdersPerRoot+1:2*nOrdersPerRoot)') = coefficients2j(1) + coefficients2j(2)*n(:) + coefficients2j(3)*(n(:)+1).^(exponent2j(1)) + coefficients2j(4)*(n(:)+1).^(exponent2j(2));
            % find second zero
            x((nOrdersPerRoot+1:2*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 2, x0, kind), n(:), x((nOrdersPerRoot+1:2*nOrdersPerRoot)'));
        end

        if k >= 3
            % coefficients and exponent are from least squares fitting the k=3,
            % n=0:10000.
            coefficients3j = [5.40770803992613;1.00093850589418;2.66926179799040;-0.174925559314932];
            exponent3j = [0.429702214054531,0.633480051735955];
            % guess for k = 3
            x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)') = coefficients3j(1) + coefficients3j(2)*n(:) + coefficients3j(3)*(n(:)+1).^(exponent3j(1)) + coefficients3j(4)*(n(:)+1).^(exponent3j(2));
            % find second zero
            x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 3, x0, kind), n(:), x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)'));
        end
    case 2
        % coefficients and exponent are from least squares fitting the k=1,
        % n=0:10000.
        coefficients1y = [0.0795046982450635;0.999998378297752;0.890380645613825;0.0270604048106402];
        exponent1y = [0.335377217953294,0.308720059086699];
        % guess for k = 1
        x((1:nOrdersPerRoot)') = coefficients1y(1) + coefficients1y(2)*n(:) + coefficients1y(3)*(n(:)+1).^(exponent1y(1)) + coefficients1y(4)*(n(:)+1).^(exponent1y(2));
        % find first zero
        x((1:nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 1, x0, kind), n(:), x((1:nOrdersPerRoot)'));

        if k >= 2
            % coefficients and exponent are from least squares fitting the k=2,
            % n=0:10000.
            coefficients2y = [1.04502538172394;1.00002054874161;-0.437921325402985;2.70113114990400];
            exponent2y = [0.434823025111322,0.366245194174671];
            % guess for k = 2
            x((nOrdersPerRoot+1:2*nOrdersPerRoot)') = coefficients2y(1) + coefficients2y(2)*n(:) + coefficients2y(3)*(n(:)+1).^(exponent2y(1)) + coefficients2y(4)*(n(:)+1).^(exponent2y(2));
            % find second zero
            x((nOrdersPerRoot+1:2*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 2, x0, kind), n(:), x((nOrdersPerRoot+1:2*nOrdersPerRoot)'));
        end

        if k >= 3
            % coefficients and exponent are from least squares fitting the k=3,
            % n=0:10000.
            coefficients3y = [3.72777931751914;1.00035294977757;2.68566718444899;-0.112980454967090];
            exponent3y = [0.398247585896959,0.604770035236606];
            % guess for k = 3
            x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)') = coefficients3y(1) + coefficients3y(2)*n(:) + coefficients3y(3)*(n(:)+1).^(exponent3y(1)) + coefficients3y(4)*(n(:)+1).^(exponent3y(2));
            % find second zero
            x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 3, x0, kind), n(:), x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)'));
        end
    otherwise
        error('Code should never get here.');
end

if k >= 4
    for iRoot = 4:k
        % guesses for remaining roots x[k] = rootSpacing + x[k-1]
        x(((iRoot-1)*nOrdersPerRoot+1:iRoot*nOrdersPerRoot)') = 2*x(((iRoot-2)*nOrdersPerRoot+1:(iRoot-1)*nOrdersPerRoot)')- x(((iRoot-3)*nOrdersPerRoot+1:(iRoot-2)*nOrdersPerRoot)');
        % find the remaining zeros
        x(((iRoot-1)*nOrdersPerRoot+1:iRoot*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, k, x0, kind), n, x(((iRoot-1)*nOrdersPerRoot+1:iRoot*nOrdersPerRoot)'));
    end
end

end



function x=findzero(n,k,x0,kind)
% Uses Halley's method to find a zero given the starting point x0
% http://en.wikipedia.org/wiki/Halley's_method

ITERATIONS_MAX = 100;       % Maximum number of iteration
TOLERANCE_RELATIVE = 1e4;   % 16-4 = 12 significant digits

% Setup loop
error = 1;
loopCount = 0;
x = 1; % Initialization value only.  It is is not used.

% Begin loop
while abs(error)>eps(x)*TOLERANCE_RELATIVE && loopCount<ITERATIONS_MAX

    switch kind
        case 1
            a = besselj(n,x0);
            b = besselj((n+1),x0);
        case 2
            a = bessely(n,x0);
            b = bessely((n+1),x0);
    end

    xSquared = x0*x0;

    error = 2*a*x0*(n*a-b*x0)/(2*b*b*xSquared-a*b*x0*(4*n+1)+(n*(n+1)+xSquared)*a*a);

    % Prepare for next loop
    x=x0-error;
    x0=x;
    loopCount=loopCount+1;

end

% Handle maximum iterations
if loopCount>ITERATIONS_MAX-1
    warning('Failed to converge to within relative tolerance of %e for n=%f and k=%d in %d iterations', eps(x)*TOLERANCE_RELATIVE, n, k, ITERATIONS_MAX);
end

end
function dJndx = dbesselj(n,x)
% DBESSELJ      A function that will generically calculate the
%               the derivative of a Bessel function of the first
%               kind of order n for all values of x.
%
% Example usage: dJndx = dbesselj(n,x);
% 
% INPUT ARGUMENTS
% ================
% n             Order of the Bessel function of the first kind
% x             Input variable to Bessel function
% 
% OUTPUT ARGUMENTS
% ================
% dJndx         Derivative of nth order Bessel function of the first 
%               kind at all values of x
    dJndx = n*besselj(n,x)./x - besselj(n+1,x);
    
end

[Torsten]

If you have little experience with MATLAB, you should first pass Matlab Onramp, an introductory course to the new software which is free of costs:

https://matlabacademy.mathworks.com/details/matlab-onramp/gettingstarted

Due to the complexity of your project, I think this could help you a lot.

[Javeria]

I have this following code for finding the roots of bessel functions , now how i can modify this one for my case. Yes R is know value.

order = 0;             % Order of the Bessel function
n_roots = 5;           % Number of positive roots desired
roots = zeros(1, n_roots);

x0 = 1;                % Initial guess
step = 5;              % Step to find next interval
found = 0;

while found < n_roots
    x1 = x0;
    x2 = x0 + step;
    if besselj(order, x1) * besselj(order, x2) < 0
        found = found + 1;
        roots(found) = fzero(@(x) besselj(order, x), [x1 x2]);
    end
    x0 = x2;
end

disp('Positive roots of J0(x):');
disp(roots);

[Torsten]

Your code is not trustworthy because it misses several roots (see below).

now how i can modify this one for my case.

What do you think is unknown in the expression you posted ? I guess everything now is available to evaluate it if you specify h, z, d, r and R.

order = 0;             % Order of the Bessel function
n_roots = 5;           % Number of positive roots desired
roots = zeros(1, n_roots);
x0 = 1;                % Initial guess
step = 5;              % Step to find next interval
found = 0;
while found < n_roots
    x1 = x0;
    x2 = x0 + step;
    if besselj(order, x1) * besselj(order, x2) < 0
        found = found + 1;
        roots(found) = fzero(@(x) besselj(order, x), [x1 x2]);
    end
    x0 = x2;
end
disp('Positive roots of J0(x):');

Positive roots of J0(x):

disp(roots);

    8.6537   18.0711   33.7758   43.1998   58.9070

m = 5;                     % maximum number of roots required
kind = 1;                  % We are dealing with bessel functions of the 1st kind
l0 = 0;                    % order of bessel function = 0
x = besselzero(l0,m,kind)  % compute the m positive roots j_{0,1:20}

x = 1×5

2.4048 5.5201 8.6537 11.7915 14.9309

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function x = besselzero(n, k, kind)
% besselzero calculates the zeros of Bessel function of the first and second kind
%
%   x = besselzero(n)
%   x = besselzero(n, k)
%   x = besselzero(n, k, kind)
%
%% Inputs
% * *n* - The order of the bessel function. n can be a scalar, vector, or
%       matrix.  n can be positive, negative, fractional, or any
%       combinaiton. abs(n) must be less than or equal to
%       146222.16674537213 or 370030.762407380 for first and second kind
%       respectively. Above these values, this algorithm will not find the
%       correct zeros because of the starting values therefore an error is
%       thrown instead.
% * k - The number of postive zeros to calculate.  When k is not supplied,
%       k = 5 is the default. k must be a scalar.
% * kind - kind is either 1 or 2. When kind is not supplied, default is
%          kind = 1.
%
%% Outputs
% * x - The calculated zeros.  size(x) = [size(n) k].
%
%% Description
% besselzero calculates the first k positive zeros of nth order bessel
% function of the first or second kind.  Note, that zero is not included as
% the first zero.
%
%% Algorithm
% the first three roots of any order bessel can be approximated by a simple
% equations.  These equations were generated using a least squares fit of
% the roots from orders of n=0:10000. The approximation is used to start
% the iteration of Halley's method.  The 4th and higher roots can be
% approximated by understanding the roots are regularly spaced for a given
% order.  Once the 2nd and 3rd roots are found, the spacing can be
% approximated by the distance between the 2nd and 3rd root.  Then again
% Halley's method can be applied to precisely locate the root.
%%
% Because the algorithm depends on good guesses of the first three zeros,
% if the guess is to far away then Halley's method will converge to the
% wrong zero which will subsequently cause any other zero to be incorrectly
% located. Therefore, a limit is put on abs(n) of 146222.16674537213 and
% 370030.762407380 for first and second kind respectively.  If n is
% specified above these limits, then an error is thrown.
%
%% Example
%   n = (1:2)';
%   k = 10;
%   kind = 1;
%   z = besselzero(n, k, kind);
%   x = linspace(0, z(end), 1000);
%   y = nan(2, length(x));
%   y(1,:) = besselj(n(1), x);
%   y(2,:) = besselj(n(2), x);
%   nz = nan(size(z));
%   nz(1,:) = besselj(n(1), z(1,:));
%   nz(2,:) = besselj(n(2), z(2,:));
%   plot(x, y, z, nz,'kx')
%

% Originally written by
% Written by: Greg von Winckel - 01/25/05
% Contact: gregvw(at)chtm(dot)unm(dot)edu
%
% Modified, Improved, and Documented by
% Jason Nicholson 2014-Nov-06
% Contact: [e-mail]

%% Change Log
% * Original release. 2005-Jan-25, Greg von Winckel.
% * Updated Documentation and commented algorithm. Fixed bug in finding the
%   the first zero of the bessel function of the second kind. Improved speed by
%   factor of 20. 2014-Nov-06, Jason Nicholson.
%

%% Input checking
assert(nargin >=1 | nargin <=3,'Wrong number of input arguments.');

% Take care of default cases of k and kind
if nargin < 2
    k = 5;
end

if nargin < 3
    kind = 1;
end

assert(isscalar(kind) & any(kind == [1 2]), '''kind''must be a scalar with value 1 or 2 only.');
assert(isscalar(k) & fix(k)==k & k>0, 'k must a positive scalar integer.');
assert(all(isreal(n(:))), 'n must be a real number.');

% negative orders have the same roots as the positive orders
n = abs(n);

% Check for that n is less than the ORDER_MAX
if kind==1
    ORDER_MAX = 146222.16674537213;
    assert(all(n <= ORDER_MAX), 'all n values must be less than or equal %6.10f for kind=1.', ORDER_MAX);
elseif kind==2
    ORDER_MAX = 370030.762407380;
    assert(all(n(:) <= ORDER_MAX), 'all n values must be less than or equal %6.10f for kind=2.', ORDER_MAX);
end
%% Setup Arrays

% output size
nSize = size(n);
if nSize(end) ==1
    outputSize = [nSize(1:end-1) k];
else
    outputSize = [nSize k];
end
% number of orders for each kth root
nOrdersPerRoot = prod(outputSize(1:end-1));
x = nan(outputSize);

%% Solve for Roots
switch kind
    case 1
        % coefficients and exponent are from least squares fitting the k=1,
        % n=0:10000.
        coefficients1j = [0.411557013144507;0.999986723293410;0.698028985524484;1.06977507291468];
        exponent1j = [0.335300369843979,0.339671493811664];
        % guess for k = 1
        x((1:nOrdersPerRoot)') = coefficients1j(1) + coefficients1j(2)*n(:) + coefficients1j(3)*(n(:)+1).^(exponent1j(1)) + coefficients1j(4)*(n(:)+1).^(exponent1j(2));
        % find first zero
        x((1:nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 1, x0, kind), n(:), x((1:nOrdersPerRoot)'));

        if k >= 2
            % coefficients and exponent are from least squares fitting the k=2,
            % n=0:10000.
            coefficients2j = [1.93395115137444;1.00007656297072;-0.805720018377132;3.38764629174694];
            exponent2j = [0.456215294517928,0.388380341189200];
            % guess for k = 2
            x((nOrdersPerRoot+1:2*nOrdersPerRoot)') = coefficients2j(1) + coefficients2j(2)*n(:) + coefficients2j(3)*(n(:)+1).^(exponent2j(1)) + coefficients2j(4)*(n(:)+1).^(exponent2j(2));
            % find second zero
            x((nOrdersPerRoot+1:2*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 2, x0, kind), n(:), x((nOrdersPerRoot+1:2*nOrdersPerRoot)'));
        end

        if k >= 3
            % coefficients and exponent are from least squares fitting the k=3,
            % n=0:10000.
            coefficients3j = [5.40770803992613;1.00093850589418;2.66926179799040;-0.174925559314932];
            exponent3j = [0.429702214054531,0.633480051735955];
            % guess for k = 3
            x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)') = coefficients3j(1) + coefficients3j(2)*n(:) + coefficients3j(3)*(n(:)+1).^(exponent3j(1)) + coefficients3j(4)*(n(:)+1).^(exponent3j(2));
            % find second zero
            x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 3, x0, kind), n(:), x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)'));
        end
    case 2
        % coefficients and exponent are from least squares fitting the k=1,
        % n=0:10000.
        coefficients1y = [0.0795046982450635;0.999998378297752;0.890380645613825;0.0270604048106402];
        exponent1y = [0.335377217953294,0.308720059086699];
        % guess for k = 1
        x((1:nOrdersPerRoot)') = coefficients1y(1) + coefficients1y(2)*n(:) + coefficients1y(3)*(n(:)+1).^(exponent1y(1)) + coefficients1y(4)*(n(:)+1).^(exponent1y(2));
        % find first zero
        x((1:nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 1, x0, kind), n(:), x((1:nOrdersPerRoot)'));

        if k >= 2
            % coefficients and exponent are from least squares fitting the k=2,
            % n=0:10000.
            coefficients2y = [1.04502538172394;1.00002054874161;-0.437921325402985;2.70113114990400];
            exponent2y = [0.434823025111322,0.366245194174671];
            % guess for k = 2
            x((nOrdersPerRoot+1:2*nOrdersPerRoot)') = coefficients2y(1) + coefficients2y(2)*n(:) + coefficients2y(3)*(n(:)+1).^(exponent2y(1)) + coefficients2y(4)*(n(:)+1).^(exponent2y(2));
            % find second zero
            x((nOrdersPerRoot+1:2*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 2, x0, kind), n(:), x((nOrdersPerRoot+1:2*nOrdersPerRoot)'));
        end

        if k >= 3
            % coefficients and exponent are from least squares fitting the k=3,
            % n=0:10000.
            coefficients3y = [3.72777931751914;1.00035294977757;2.68566718444899;-0.112980454967090];
            exponent3y = [0.398247585896959,0.604770035236606];
            % guess for k = 3
            x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)') = coefficients3y(1) + coefficients3y(2)*n(:) + coefficients3y(3)*(n(:)+1).^(exponent3y(1)) + coefficients3y(4)*(n(:)+1).^(exponent3y(2));
            % find second zero
            x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, 3, x0, kind), n(:), x((2*nOrdersPerRoot+1:3*nOrdersPerRoot)'));
        end
    otherwise
        error('Code should never get here.');
end

if k >= 4
    for iRoot = 4:k
        % guesses for remaining roots x[k] = rootSpacing + x[k-1]
        x(((iRoot-1)*nOrdersPerRoot+1:iRoot*nOrdersPerRoot)') = 2*x(((iRoot-2)*nOrdersPerRoot+1:(iRoot-1)*nOrdersPerRoot)')- x(((iRoot-3)*nOrdersPerRoot+1:(iRoot-2)*nOrdersPerRoot)');
        % find the remaining zeros
        x(((iRoot-1)*nOrdersPerRoot+1:iRoot*nOrdersPerRoot)') = arrayfun(@(n, x0) findzero(n, k, x0, kind), n, x(((iRoot-1)*nOrdersPerRoot+1:iRoot*nOrdersPerRoot)'));
    end
end

end



function x=findzero(n,k,x0,kind)
% Uses Halley's method to find a zero given the starting point x0
% http://en.wikipedia.org/wiki/Halley's_method

ITERATIONS_MAX = 100;       % Maximum number of iteration
TOLERANCE_RELATIVE = 1e4;   % 16-4 = 12 significant digits

% Setup loop
error = 1;
loopCount = 0;
x = 1; % Initialization value only.  It is is not used.

% Begin loop
while abs(error)>eps(x)*TOLERANCE_RELATIVE && loopCount<ITERATIONS_MAX

    switch kind
        case 1
            a = besselj(n,x0);
            b = besselj((n+1),x0);
        case 2
            a = bessely(n,x0);
            b = bessely((n+1),x0);
    end

    xSquared = x0*x0;

    error = 2*a*x0*(n*a-b*x0)/(2*b*b*xSquared-a*b*x0*(4*n+1)+(n*(n+1)+xSquared)*a*a);

    % Prepare for next loop
    x=x0-error;
    x0=x;
    loopCount=loopCount+1;

end

% Handle maximum iterations
if loopCount>ITERATIONS_MAX-1
    warning('Failed to converge to within relative tolerance of %e for n=%f and k=%d in %d iterations', eps(x)*TOLERANCE_RELATIVE, n, k, ITERATIONS_MAX);
end
end

Answer 4

[Javeria]

Hello , I modefied the above scripts for my problem, So is my this following code is okay or not?

function x = bessel0j(n,q,opt)
% a row vector of the first q roots of bessel function Jn(x), integer order.
% if opt = 'd', row vector of the first q roots of dJn(x)/dx, integer order.
% if opt is not provided, the default is zeros of Jn(x).
% all roots are positive, except when n=0,
% x=0 is included as a root of dJ0(x)/dx (standard convention).
%
% starting point for for zeros of Jn was borrowed from Cleve Moler, 
% but the starting points for both Jn and Jn' can be found in
% Abramowitz and Stegun 9.5.12, 9.5.13. 
%
% David Goodmanson
%
% x = bessel0j(n,q,opt)

k = 1:q;
if nargin==3 && opt=='d'
    beta = (k + n/2 - 3/4)*pi;
    mu = 4*n^2;
    x = beta - (mu+3)./(8*beta) - 4*(7*mu^2+82*mu-9)./(3*(8*beta).^3);
    for j=1:8
        xnew = x - besseljd(n,x)./ ...
            (besselj(n,x).*((n^2./x.^2)-1) -besseljd(n,x)./x);
        x = xnew;    
    end
    if n==0
        x(1) = 0;    % correct a small numerical difference from 0
    end
else
    beta = (k + n/2 - 1/4)*pi;
    mu = 4*n^2;
    x = beta - (mu-1)./(8*beta) - 4*(mu-1)*(7*mu-31)./(3*(8*beta).^3);
    for j=1:8
        xnew = x - besselj(n,x)./besseljd(n,x);
        x = xnew;
    end
end
end
% --- Local helper function for derivative of Bessel function ---
function dJ = besseljd(n, x)
    dJ = 0.5 * (besselj(n - 1, x) - besselj(n + 1, x));
end
% --- Bessel Root Test Script ---

% Parameters
l = 0;              % Bessel order
q = 20;             % Number of roots
R = 4.0;            % Radius
use_derivative = false;  % Set to true for J_l'(x)

% Compute Bessel roots
if use_derivative
    x = bessel0j(l, q, 'd');        % Derivative roots
    root_type = sprintf("J_%d'(x)", l);
else
    x = bessel0j(l, q);             % Function roots
    root_type = sprintf("J_%d(x)", l);
end

% Compute chi_n^l = x_n / R
chi_nl = x / R;

% Display results
fprintf('First %d chi_n^%d values (R = %.2f):\n', q, l, R);

First 20 chi_n^0 values (R = 4.00):

disp(chi_nl);

  Columns 1 through 18

    0.6012    1.3800    2.1634    2.9479    3.7327    4.5178    5.3029    6.0881    6.8734    7.6587    8.4440    9.2293   10.0146   10.7999   11.5853   12.3707   13.1560   13.9414

  Columns 19 through 20

   14.7267   15.5121

[Torsten]

I trust in @David Goodmanson code :-)

Where do you need the roots of the derivative of the Bessel function of order l ? From what you posted, you only need to evaluate the derivative at R*(roots of the Bessel function of order l). Thus the "use_derivative" part of your code would be superfluous. But maybe it's needed somewhere else for your problem...

[David Goodmanson]

When Torsten [ thanks for your comment :-) ] pointed out not needing the roots of the derivative function J', I realized that I had made a big mistake in one of the comments after my answer. I went back and corrected it, but I think it got Javeria partly off on the wrong track.

The functions themselves are all ok, but I misapplied one of them. Based on the expression in your original question, first you can find

x = bessel0j(l, q);             % Function roots          x = gamma/R

but for the J' part you need the function J' evaluated at the roots you just found for J. That is not

x = bessel0j(l, q, 'd');        % Derivative roots             incorrect 

but rather, with your local helper function

Jd(l, x);        % J' evaluated at the roots of J

and these are the values for the denominator of the original question.