MATLAB Answers

Bessel function - roots and zeros on interval <12,20>

49 views (last 30 days)
Jaroslav Zadnik
Jaroslav Zadnik on 30 Jan 2020
Commented: David Goodmanson on 21 Dec 2020
Hi, I'm beginner here. I have Bessel function of the first kind J1. I need find roots and zeros on interval <10,20>. Could you help me plase? :/

  3 Comments

Walter Roberson
Walter Roberson on 30 Jan 2020
I suggest that you plot the function and use what you see for estimates of starting points for fzero
Jaroslav Zadnik
Jaroslav Zadnik on 30 Jan 2020
I have only:
x=10:0.1:20;
i=1;
J=besselj(i,x)
plot(x,J)
John D'Errico
John D'Errico on 31 Jan 2020
You are inconsistent, since in the title you ask about the interval [12,20], yet in your question, you say [10,20].

Sign in to comment.

Answers (2)

David Goodmanson
David Goodmanson on 31 Jan 2020
Edited: David Goodmanson on 1 Feb 2020
HI Jaroslav
Do roots and zeros mean the same thing? The code below uses Newton's method to find the roots, which has the advantage that you can apply it to a vector of roots and get them all at once. I don't have a tolerance check for the number of times to iterate Newton, but j = 1:8 seems to work (good up to at least the 14th decimal place for all roots in this case).
function x = besselzeroj(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 Jn(x).
% includes the root x = 0 for dJ0(x)/dx (zeroth order)
% but not for higher orders.
%
% starting point for for zeros of Jn is borrowed from Cleve Moler,
% but starting points for both Jn and Jn' can be found in
% Abramowitz and Stegun 9.5.12, 9.5.13.
%
% David Goodmanson
%
% function x = besselzeroj(n,q,opt)
k = 1:q;
if nargin==3 & opt=='d'
% find starting points for J'
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);
% Newton
for j=1:8
xnew = x - besseljdv(n,x)./ ...
(besselj(n,x).*((n^2./x.^2)-1) -besseljdv(n,x)./x);
x = xnew;
end
else
% find starting points for J
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);
% Newton
for j=1:8
xnew = x - besselj(n,x)./besseljdv(n,x);
x = xnew;
end
end
end
function y = besseljdv(n,x)
% derivative of bessel function of integer order,
% plain vanilla version
%
% function y = besseljdv(n,x)
y = -besselj(n+1,x) + n*besselj(n,x)./x;
% get rid of nans, integer case
if n==1
y(x==0) = 1/2;
else
y(x==0) = 0;
end
end
format long
besselzeroj(1,12)'
ans =
3.831705970207512
7.015586669815619
10.173468135062722
13.323691936314225
16.470630050877634
19.615858510468243
22.760084380592772
25.903672087618382
29.046828534916855
32.189679910974405
35.332307550083868
38.474766234771614

  2 Comments

sars sabina
sars sabina on 20 Dec 2020
Could you please explain more about the starting point that you used in your code?
% find starting points for J
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);
David Goodmanson
David Goodmanson on 21 Dec 2020
Hi s s,
I can't say much more than what is in the answer, that the expression for x is in Abramowitz and Stegun, Handbood of Mathematical Functions 9.5.12. I used just the first three terms in 9.5.12, which is an asymptotic expansion that gets more accurate the larger the roots are. But it turns out that even for small roots, the three-term expression does well at providing a first guess. In the expression, the first term can be derived fairly easily by evaluating a certain integral by the method of steepest descents. I think I could probably get the second term if I worked at it enough. After that, not happening.

Sign in to comment.


John D'Errico
John D'Errico on 31 Jan 2020
This is easier than you think. :) Of course, that often is true, once you know the answer.
By the way, it can be a dangerous idea in general to use i as a variable name, since i is also sqrt(-1) by default in MATLAB. One day that will be important for you.
So what you have is a start:
x = 10:0.1:20;
n = 1;
J = besselj(n,x);
plot(x,J)
Now, do you know of some tool that can find a root, given a start point? Better yet, what if you have an interval in x that brackets a root? Note that fzero does exactly what you need, if you have an interval that brackets the root.
Is there a way you can locate pairs of points that bracket the root? What does that mean? Essentially, you need to find a pair of points x(m) and x(m+1), such that J(m) and J(m+1) have opposite signs. The sign function can help you here. For example:
m = find(diff(sign(J)) ~= 0)
m =
2 34 65 97
[J(m);J(m+1)]
ans =
0.018395515457572 -0.005177480554671 0.013894680686257 -0.002856572403405
-0.006615743297723 0.016599019864009 -0.005764213735631 0.015100612097755
So it appears there will be 4 roots in the interval of interest, because there are 4 zero crossings. Here is the first root:
fun = @(x) besselj(n,x);
[xvl,fval] = fzero(fun,[x(m(1)),x(m(1) + 1)])
xvl =
10.173468135062722
fval =
2.805130098197981e-16
Do the same call for each of the zero crossings identified in m.

  0 Comments

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!