Guess update for Bessel function

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!

[David Goodmanson]

Hi JS,

what is the modes function?

[JS]

Hi David,

Sorry about that, I forgot to include it in my original post. Modes is a vector.

modes = [19,18,17,16,15,14];

Answer 1

[David Goodmanson]

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

[JS]

Thanks David! It worked perfectly. I implemented it in my main code, but in an inefficient manner by repeating the same section of your code where the roots are solved for as j changes with if statements. I ran into sizing issues or exceeding the indexes so I left it as is using if statements for every different j value. I wasn't sure where in your code the function is being converted into a square wave. Is it the section where you find the approximate zeros? Also, when I plot the zero crossings outside of the loop, they don't match up to the plots when plotting within the loop. This is what ended up doing:

R_i = 3;
R_o = 13;
modes = [19,18,17,16,15,14];

for j = 1:4
    
        m = modes(j);
        k = modes(j+2);
        
        % bessel function
        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)));
        
        % interval search
        mu = 0.01:0.01:4;
        
        if j == 1
            
            % estimate where zero crossings occur
            f = find(abs(diff(sign(fun(mu))))>1);
            mu_approx = mu(f);
            mu0 = [1 (mu_approx(1:end-1)+mu_approx(2:end))/2];
            
            % find the roots
            mu1 = [];
            
            for n = 1:length(mu0)-1
                
                mu1 = [mu1 fzero(fun,mu0(n:n+1))];
                
            end
            
            first = mu1(:,1:5)';
            
            % plot the roots
            mu = 0:.001:3;
            mu_plot = mu1(mu1 < 3);
            figure(1)
            plot(mu,mu.^15.*fun(mu),mu_plot,zeros(size(mu_plot)),'o')
            grid on
            
        end
        
        if j == 2
            
            % estimate where zero crossings occur
            f = find(abs(diff(sign(fun(mu))))>1);
            mu_approx = mu(f);
            mu0 = [1 (mu_approx(1:end-1)+mu_approx(2:end))/2];
            
            % find the roots
            mu2 = [];
            
            for n = 1:length(mu0)-1
                
                mu2 = [mu2 fzero(fun,mu0(n:n+1))];
                
            end
            
            second = mu2(:,1:5)';
            
            % plot the roots
            mu = 0:.001:3;
            mu_plot = mu2(mu2 < 3);
            figure(2)
            plot(mu,mu.^15.*fun(mu),mu_plot,zeros(size(mu_plot)),'o')
            grid on
            
        end
        
        if j == 3
            
            % estimate where zero crossings occur
            f = find(abs(diff(sign(fun(mu))))>1);
            mu_approx = mu(f);
            mu0 = [1 (mu_approx(1:end-1)+mu_approx(2:end))/2];
            
            % find the roots
            mu3 = [];
            
            for n = 1:length(mu0)-1
                
                mu3 = [mu3 fzero(fun,mu0(n:n+1))];
                
            end
            
            third = mu3(:,1:5)';
            
            % plot the roots
            mu = 0:.001:3;
            mu_plot = mu3(mu3 < 3);
            figure(3)
            plot(mu,mu.^15.*fun(mu),mu_plot,zeros(size(mu_plot)),'o')
            grid on
            
        end
        
        if j == 4
            
            % estimate where zero crossings occur
            f = find(abs(diff(sign(fun(mu))))>1);
            mu_approx = mu(f);
            mu0 = [1 (mu_approx(1:end-1)+mu_approx(2:end))/2];
            
            % find the roots
            mu4 = [];
            
            for n = 1:length(mu0)-1
                
                mu4 = [mu4 fzero(fun,mu0(n:n+1))];
                
            end
            
            fourth = mu4(:,1:5)';
            
            % plot the roots
            mu = 0:.001:3;
            mu_plot = mu4(mu4 < 3);
            figure(4)
            plot(mu,mu.^15.*fun(mu),mu_plot,zeros(size(mu_plot)),'o')
            grid on
            
        end
        
        X = fzero(fun,mu_approx(j)); 
    
end

%Eigenvalues
n = [1;2;3;4];
ev = table(n,first,second,third,fourth);
disp(ev)

[David Goodmanson]

Hello JS,

in this line

in = find(abs(diff(sign(fun(mu))))>1);  

sign creates the square wave, diff detects the square wave jumps between +1 and -1 (or vise versa), and find(abs(... gives the array indices of the approximated zero crossings.

The region 0<mu<3 has a different number of roots for different values of j. The code below stores the first ten roots for each j, and deals with 0<mu<3 afterwards, in the part where the figures are generated. The roots are stored in a 4x10 matrix mu0j with roots values across. j values down. This is a far better approach than creating variable mu1, mu2, etc. (dynamic variable naming). For one thing, you don't have to repeat code four times, you can just fill the matrix by using a single for loop. For another, what if you wanted to do this for 20 different values of j? Dynamic variable naming does not generalize well.

Displaying mu0j shows a nice summary of the roots, which of course you can also do with a table.

R_i = 3;
R_o = 13;
modes = [19,18,17,16,15,14];
nroots = 10;
mu0j = zeros(4,nroots);
for j = 1:4;
  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:5;
  % find approximate zeros, find peak locations between zeros
  fin = find(abs(diff(sign(fun(mu))))>1);      
  fin = fin(1:nroots+1);
  mu0approx = mu(fin);
  mu1 = [1 (mu0approx(1:end-1)+mu0approx(2:end))/2];
  % find mu0 values, the roots
  for n = 1:length(mu1)-1
    mu0 = fzero(fun,mu1(n:n+1));
    mu0j(j,n) = mu0;
  end
  
  % plot some zeros
  mu = 0:.001:3;
  mu0plot = mu0j(j,:);
  mu0plot = mu0plot(mu0plot<3);
  figure(j)
  plot(mu,mu.^15.*fun(mu),mu0plot,zeros(size(mu0plot)),'o')
  grid on
 end