recreating roots of a derivative of bessel funtion of first order

Hello , i want to recreate the root table which are the derivative of a bessel function of the first order.

i think besselj(nu,Z) and fzero plus derivative is the combination.

but i dont know how exactly to formula the combination so ill get the table bellow?

Thanks.

[Walter Roberson]

https://www.mathworks.com/matlabcentral/answers/93677-how-can-i-evaluate-the-derivatives-of-a-bessel-function-at-different-points#answer_103028

Answer 1

[David Goodmanson]

Hi fv,

Here is a function for the first q zeros of both Jn(x) and dJn(x) /dx. As an example, to find the first 100 zeros of the derivative of J_5(x) takes a couple of milliseconds.

< minor improvements to bessel0j since first posted >
 
function x = bessel0j(n,q,opt)

% first q roots of bessel function Jn(x), integer order.
% if opt = 'd', 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 (a standard convention).
%
% starting point for for zeros of Jn is 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
%
% function 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  % end of function

% -------------------------------------------------------
 
function y = besseljd(n,x,in1,in2);

% derivative of bessel function of integer order
% if type = '+', then J(n,x)' = -J(n+1,x) + (n/x)*J(n,x)
% if type = '-', then J(n,x)' =  J(n-1,x) - (n/x)*J(n,x)
% default is '+' 
% if s = 1, result is scaled by exp(-abs(imag(z))), same as with besselj.
% default is 0, no scaling 
% input order of s and type does not matter, and either
% or both can be omitted, no placeholder required 
% 
% function y = besseljd(n,x,type,s);

type = '+';  s = 0;
if nargin==4
  if ~ischar(in1)
    s = in1; type = in2;
  else
    type = in1; s = in2;    
  end
elseif nargin==3
  if ~ischar(in1)
    s = in1;
  else
    type = in1;
  end
end

if type=='+'
   y = -besselj(n+1,x,s) + n*besselj(n,x,s)./x;
else
   y =  besselj(n-1,x,s) - n*besselj(n,x,s)./x; 
end

% get rid of nans, integer case so far
if n==1
  y(x==0) = 1/2;
else
  y(x==0) = 0;
end

% 'if'check is not required for newer versions, but at one time besselj
% had a bug, for integer n~=0 and real negative x, output was real + 0i
if isint(n) & isreal(x)
  y = real(y);
end 

end   % end of function

[Walter Roberson]

@David Goodmanson

That looks like it has a structure similar to Newton's Method ? Does it improve by 2 bits per iteration so 8 iterations is enough?

besseljd appears to be from https://www.mathworks.com/matlabcentral/fileexchange/27819-optical-fibre-toolbox ?

[David Goodmanson]

Hi Walter,

It does use Newton's method, which is nice since the iterations can be done vectorally. The trick is to get all q of the startiing points right so that for each would-be zero, the method heads for the right location for every zero from 1 to q. To do things really correctly I guess I should put in a convergence criterion, and the code is fast enough that a convergence check for all q of them should not have a huge impact on the run time. But I did some tests and 8 iterations keeps the code simple and seems to work fine in practical cases, although I don't know how it is doing in bits per iteration.

You're right, I forgot to include besseljd which is a function I wrote so I will go back and include it in the answer.

[fima v]

Hello, So basicly in need to use besseljd function.

i need to find zeros of it.basicly i need to recreate my original table.

func(0,1)=3.832 func(1,1)=1.841.

How do create such a function wrapper?

Thanks.

[David Goodmanson]

The last two rows of the table are bessel0j(n,3,'d') for n = 1,2. For the first row with n=0, the code has 0 as the first root, but the table is ignoring the zero, so you can do something like xroots = bessel0j(0,4,'d'); xroots = xroots(2:4); to eliminate the zero. Then concatenate the rows vertically using, say, vertcat to create a 3x3 matrix that's the same as the table.