recreating roots of a derivative of bessel funtion of first order
16 views (last 30 days)
Show older comments
Hello , i want to recreate the root table which are the derivative of a bessel function of the first order.
but i dont know how exactly to formula the combination so ill get the table bellow?
Thanks.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1389394/image.jpeg)
Answers (1)
David Goodmanson
on 21 May 2023
Edited: David Goodmanson
on 22 May 2023
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
4 Comments
David Goodmanson
on 22 May 2023
Edited: David Goodmanson
on 22 May 2023
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.
See Also
Categories
Find more on Bessel functions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!