MATLAB: Guess update for Bessel function

besselbessel functioneigenvaluefzeroguessMATLABzero

Hi,
I'm trying to figure out a way to properly update the guess (from the initial within a range of 1 to 3) for a bessel function (1st and 2nd kind) that I'm solving roots/eigenvalues for as I go through the loop. Below is a section of my code that I'm currently having trouble in trying to initiate the correct second guess where the function is zero:
R_i = 3;
R_o = 13;
x0 = 1;
for i = 1:5
    for j = 1:4
        
        m = modes(j);
        k = modes(j+2);
        
        x = x0;
        fun = @(mu) (0.5*(besselj(m,mu*R_o)-besselj(k,mu*R_o)))*(0.5*(bessely(m,mu*R_i)-bessely(k,mu*R_i))) ...
            - (0.5*(besselj(m,mu*R_i)-besselj(k,mu*R_i)))*(0.5*(bessely(m,mu*R_o)-bessely(k,mu*R_o)));
        eig(i,j) = fzero(fun,x); %solve for the zeros of the function
        mu_guess(i,j) = eig(i,j); %storing all eigenvalues
        
        x0 = i; %update guess per loop
        
    end
end
Appreciate any help I can get. Thanks!

Best Answer

Hi JS,
For j = 1, figure 1 shows the general behavior of the function. Since it's a plot of log(abs(fun)) there can be no zero crossings but you can see the attempts to head for zero. There are a lot of roots but no roots below mu = 1 which is also true for j = 2,3,4.
Since it's an analytic function (no noise) and the roots are very close to equally spaced (which will become evident in figure 2), a fairly simple method can be used to find the approximate locations of the zeros: convert the function to a square wave and find the discontinuities.
With fzero. rather than feed it an approximate location for a zero, it is safer to give it an upper and lower bound where there is sure to be exactly one zero crossing in between. The code finds the points halfway between the approximate zeros, and adds 1 at the beginning since there is no zero to the left of that.
Plot 2 shows a result for j = 1. The function is dropping in amplitude so rapidly in that region that it helps to multiply by a factor of mu.^15 just to make it easy to see. That of course does not affect the location of the zeros.
Results for j = 2,3,4 are similar.
R_i = 3;
R_o = 13;
modes = [19,18,17,16,15,14];
j = 1;
m = modes(j);
k = modes(j+2);
fun = @(mu) (0.5*(besselj(m,mu*R_o)-besselj(k,mu*R_o))).*(0.5*(bessely(m,mu*R_i)-bessely(k,mu*R_i))) ...
          - (0.5*(besselj(m,mu*R_i)-besselj(k,mu*R_i))).*(0.5*(bessely(m,mu*R_o)-bessely(k,mu*R_o)));
mu = .01:.01:50;
figure(1)
loglog(mu,abs(fun(mu)))
grid on
% find approximate zeros, find peak locations between zeros
fin = find(abs(diff(sign(fun(mu))))>1);      
mu0approx = mu(fin);
mu1 = [1 (mu0approx(1:end-1)+mu0approx(2:end))/2];
% find mu0 values, the roots
mu0 = [];
for n = 1:length(mu1)-1
  mu0 = [mu0 fzero(fun,mu1(n:n+1))];
end
% plot some zeros
mu = 0:.001:6;
mu06 = mu0(mu0<6);
figure(2)
plot(mu,mu.^15.*fun(mu),mu06,zeros(size(mu06)),'o')
grid on
Related Question