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

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? :/

[Walter Roberson]

I suggest that you plot the function and use what you see for estimates of starting points for fzero

[Jaroslav Zadnik]

I have only:

x=10:0.1:20;

i=1;

J=besselj(i,x)

plot(x,J)

[John D'Errico]

You are inconsistent, since in the title you ask about the interval [12,20], yet in your question, you say [10,20].

Answer 1

[David Goodmanson]

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

[sars sabina]

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]

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.

Answer 2

[John D'Errico]

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.