bessel funtion drawing problem part two

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')

[Sam Chak]

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

j_minus = 1501×1

0 -0.0050 -0.0100 -0.0150 -0.0200 -0.0250 -0.0300 -0.0350 -0.0400 -0.0450

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Term1 = 1501×1

0 -0.0050 -0.0100 -0.0150 -0.0200 -0.0250 -0.0300 -0.0350 -0.0400 -0.0450

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Term2 = 1.2246e-16

Term3 = 1501×1

0 0 0 0 0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

j_minus = 1501×1

0 0.0000 0.0000 0.0001 0.0002 0.0003 0.0004 0.0006 0.0008 0.0010

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Term1 = 1501×1

0 0.0000 0.0000 0.0001 0.0002 0.0003 0.0004 0.0006 0.0008 0.0010

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Term2 = -2.4493e-16

Term3 = 1501×1

0 0 0 0 0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

legend('Y0','Y1','Y2')

[Torsten]

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]

@Torsten

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.

[Torsten]

@Walter Roberson

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]

sometimes I wish people would look at some of the other already-posted answers.

Answer 1

[David Goodmanson]

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.