bessel funtion drawing problem part two
Show older comments
I am trying to draw the Neumann Bessel Function. I understand why I need to use harmonic function to draw the Y0(x), but why can't I draw Y1(x) and Y2(x).
%Neumann Bessel Function
x=(0:0.01:15).';
i=0:50;
hold on
set(gca,'YLim',[-1,1]); %range
set(gca,'XLim',[0,15]); %domain
for m=0:2
j=sum((((-1).^i)./(factorial(i).*factorial(i+m))).*(x./2).^(2*i+m),2);
if m ==0
y = bessely(m, x);
else
i_minus=m:50;
j_minus=sum((((-1).^i_minus)./(factorial(i_minus).*factorial(i_minus-m))).*(x./2).^(2*i_minus-m),2);
y=(j.*cos(m*pi)-j_minus)./(sin(m*pi));
end
plot(x,y)
end
legend('Y0','Y1','Y2')
5 Comments
The mathematical reason is that when m = 1:2, the computation of j.*cos(m*pi)-j_minus is effectively equal to zero—absolute zero!
%Neumann Bessel Function
x=(0:0.01:15).';
i=0:50;
hold on
set(gca,'YLim',[-1,1]); %range
set(gca,'XLim',[0,15]); %domain
for m=0:2
j=sum((((-1).^i)./(factorial(i).*factorial(i+m))).*(x./2).^(2*i+m),2);
if m == 0
y(:,m+1)= bessely(m, x);
else
i_minus = m:50;
j_minus = sum((((-1).^i_minus)./(factorial(i_minus).*factorial(i_minus-m))).*(x./2).^(2*i_minus-m),2)
Term1 = j.*cos(m*pi)
Term2 = sin(m*pi)
Term3 = Term1 - j_minus
y(:,m+1)= (Term1 - j_minus)./Term2;
end
plot(x, y(:,m+1))
end
legend('Y0','Y1','Y2')
You divide by 0 when you divide by sin(m*pi). Thus for m being an integer, Y is defined as the limit for m -> 0,1,2. Direct insertion of m in the formula for Y will not work. And when you try it with m near to 0,1 or 2, you have to replace factorial(i+m) by gamma(i+1+m).
Walter Roberson
on 30 Mar 2025
Suppose i = 0 and m = 0, then factorial(i+m) would be factorial(0) . factorial(0) is well-defined, being 1.
i and m are both non-negative integers, so their sum must be a non-negative integer, and factorial() is defined for the non-negative integers.
Inserting m as an integer directly in the definition of "bessely" leads to a division by zero because of the sin(m*pi) in the denominator. So I thought of approaching the integer m-values by inserting m + eps with a small value for eps in the function definition. This would make it necessary to use the gamma function.
David Goodmanson
on 30 Mar 2025
Edited: David Goodmanson
on 30 Mar 2025
sometimes I wish people would look at some of the other already-posted answers.
Answers (1)
David Goodmanson
on 29 Mar 2025
Hello Zeyuan,
You are using the expression
Y(nu,z) = (J(nu,z)*cos(nu*pi) - J(-nu,z))/sin(nu(z))
which does not work when nu is an integer. In those cases the denominator is zero, and you have a 0/0 form. The only reason you got 0 for an answer is that
sin(pi) = 1.2246e-16
in double precision arithmetic, not zero. Otherwise you would have gotten 0/0 = NaN everywhere.
The expression does work as a limit when nu --> integer. Following code uses nu = 1+1e-6
to find Y(1,z), and it calculates J by power series, similar to what you are doing:
kterms = 0:40;
zarray = .08:.01:22;
[z k] = meshgrid(zarray,kterms);
nu = 1;
delt = 1e-6;
nua = nu + delt;
Jmatp= (z/2).^nua.*(-z.^2/4).^k./(gamma(k+1).*gamma(nua+k+1));
Jp = sum(Jmatp);
Jmatm = (z/2).^(-nua).*(-z.^2/4).^k./(gamma(k+1).*gamma(-nua+k+1));
Jm = sum(Jmatm);
Y = (Jp*cos(nua*pi)-Jm)/sin(nua*pi);
figure(1)
plot(zz,Y,zz,bessely(nu,zz))
grid on
The calculated Y in blue finally starts to differ from the bessely function at around z = 21. Making delt smaller, say 1e-7, does not help; it actually makes things worse. You can experiment around to see what works. Not coincidentally the power series to compute J in the first place starts to go bad up around this region. When z gets large enough it's time to go to a different expression for J.
Categories
Find more on Special Functions in Help Center and File Exchange
on 30 Mar 2025
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
