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

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.