Guess update for Bessel function
2 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
David Goodmanson
2020-12-12
编辑:David Goodmanson
2020-12-15
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 个评论
David Goodmanson
2020-12-15
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
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Bessel functions 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!