Finding the nth root of the bessel function of the first kind (J0(x)) using bisection

Question

First of all, I would just like to clarify that this is an assignment for school, so I am not looking for a solution. I simply want to be pushed in the right direction.

Now, for the problem.

We have code for finding the root of a polynomial using bisection:

function [root, niter, rlist] = bisection2( func, xint, tol )
% BISECTION2: Bisection algorithm for solving a nonlinear equation
%             (non-recursive).
%
% Sample usage:
%   [root, niter, rlist] = bisection2( func, xint, tol )
% 
%  Input:
%     func    - function to be solved
%     xint    - interval [xleft,xright] bracketing the root
%     tol     - convergence tolerance (OPTIONAL, defaults to 1e-6)
%
%  Output:
%     root    - final estimate of the root
%     niter   - number of iterations needed  
%     rlist   - list of midpoint values obtained in each iteration.

% First, do some error checking on parameters.
if nargin < 2
  fprintf( 1, 'BISECTION2: must be called with at least two arguments' );
  error( 'Usage:  [root, niter, rlist] = bisection( func, xint, [tol])' );
end
if length(xint) ~= 2, error( 'Parameter ''xint'' must be a vector of length 2.' ), end  
if nargin < 3, tol = 1e-6; end

% fcnchk(...) allows a string function to be sent as a parameter, and
% coverts it to the correct type to allow evaluation by feval().
func = fcnchk( func );

done  = 0;
rlist = [xint(1); xint(2)];
niter = 0;

while ~done

 % The next line is a more accurate way of computing
 % xmid = (x(1) + x(2)) / 2 that avoids cancellation error.
 xmid = xint(1) + (xint(2) - xint(1)) / 2;
 fmid = feval(func,xmid);

 if fmid * feval(func,xint(1)) < 0
   xint(2) = xmid;
 else
   xint(1) = xmid;
 end

 rlist = [rlist; xmid];
 niter = niter + 1;

 if abs(xint(2)-xint(1)) < 2*tol || abs(fmid) < tol
   done = 1;
 end
end

root = xmid;
%END bisection2.

We must use this code to find the nth zero of a Bessel function of the first kind (J0(x)). It is quite simple to insert a range and then find the specific root we are looking for. However, we must plot Xn vs. n and for that, we would need to be able to calculate a large number of roots in relation to n. So for that, I wrote this code:

bound = 1000;
x = linspace(0, bound, 1000);
for i=0:bound
    for j=1:bound
        y = bisection2(@(x) besselj(0,x), [i,j], 1e-6)
    end
end

I believed this would work, but the roots it provides are not in order and keep repeating. The problem I believe is my range when I call bisection2. I know [i,j] is not the best way to do it and was hoping someone could lead me in the right direction on how to fix this problem.

Thank you.

Answer 1

Your implementation is in the right direction but it isn't completely correct.

bound = 1000;
% x = linspace(0, bound, 1000);  No need of this line.
x_ini = 0; n =1; 
Root = zeros(bound+1,100);  % Assuming there are 100 roots in [1,1000] range

for i=0:bound                                            
  for j=1:bound                                          
    y = bisection2(@(x) besselj(i,x), [x_ini,j], 1e-6);   % besselj(i,x) = Ji(x); i = 0,1,2,3,...
    if abs(bessel(i,y)) < 1e-6
       x_ini = y;                                         % Finds nth root
       Root(i+1,n) = y;
       n = n+1;
    end
  end
  n = 1;
end

I have replaced besselj(0,x) in your code with besselj(i,x). This gives you roots for not only J0(x) but also for J1(x), J2(x), J3(x), .....so on. (i=0,1,2,... )

Another change I made in your code is replacing [i,j] by [x_ini,j]. Initially x_ini=0 & j=1. This tries to find root in the interval [0,1]. Since the 1st root for J0 occurs at 2.4, the root your bisection function calculates (0.999) isn't actually the first root. The lines between if.....end will check if the root found by Bisection function is actually a root or not. If it is, x_ini will take the value of the root as the next root will occur after x = x_ini (or y).