Guess update for Bessel function

2 views (last 30 days)
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!
  2 Comments
David Goodmanson
David Goodmanson on 11 Dec 2020
Hi JS,
what is the modes function?
JS
JS on 11 Dec 2020
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];

Sign in to comment.

Accepted Answer

David Goodmanson
David Goodmanson on 12 Dec 2020
Edited: David Goodmanson on 15 Dec 2020
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
  2 Comments
JS
JS on 13 Dec 2020
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
David Goodmanson on 15 Dec 2020
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

Sign in to comment.

More Answers (0)

Products


Release

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!