recreating roots of a derivative of bessel funtion of first order
30 次查看(过去 30 天)
显示 更早的评论
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.
回答(1 个)
David Goodmanson
2023-5-21
编辑:David Goodmanson
2023-5-22
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 个评论
David Goodmanson
2023-5-22
编辑:David Goodmanson
2023-5-22
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Bessel functions 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!