Bessel function of the first kind

My problem is to write a program which calculates a Bessel function of the first kind using the formula:

Jn+1(x) + Jn−1(x) = (2n/x)*Jn(x)

This is to be computed enough times to attain all Jn(x) values up to n = 20. What we know is that x = 1, J0 = 0.76519768655796655145 and J1 = 0.44005058574493351596.

The code I have has a bug which I am unable to figure out. Instead of J20 ~ zero, half way through the iterations the value of Jn(x) begins increasing.

a = 0.76519768655796655145;

b = 0.44005058574493351596;

for n = 1:20

c = abs(2*b*n - a);

disp(c)

temp = b;

part = c;

a = temp;

b = part;

end

I cannot see what is wrong with the code, so I would appreciate any help.

Answer 1

[Torsten]

For an explanation, google "Bessel's maze" and take the first hit.

Best wishes

Torsten.

Answer 2

[Ced]

Hi

interesting question, thanks.

The short answer is: there is (almost) nothing wrong with the code, except the abs, which I don't think should be there.

The reason for the divergence however is numerical.

1) Your initial values, despite having a high precision, are not precise enough.

2) after a few iterations, even with precise initial values, you run into numerical problems.

This code will let you see how the values evolve:

N_max = 20;
bessel_val = zeros(N_max+1,1);
bessel_val_real = zeros(N_max+1,1);

bessel_val(1) = besselj(0,1); %0.76519768655796655145;
bessel_val(2) = besselj(1,1); %0.44005058574493351596;

bessel_val_real(1) = besselj(0,1);
bessel_val_real(2) = besselj(1,1);

for n = 2:N_max
  bessel_val(n+1) = (2*(n-1)*bessel_val(n) - bessel_val(n-1));
  bessel_val_real(n+1) = 2*(n-1)*besselj(n-1,1) - besselj(n-2,1);
  fprintf('recursive, naive: J%i(1) = %14.10g\n',n,bessel_val(n+1))
  fprintf('recursive, true:  J%i(1) = %14.10g\n',n,bessel_val_real(n+1))
end

Here, the "naive" version computes the recursion using only the original J0 and J1 values. The "true" version computes the recursion for a single step, i.e. starting with the actual function values returned by the matlab implementation.

Cheers

PS: Part of the output of the code above:

recursive, naive: J2(1) =   0.1149034849 
recursive, true:  J2(1) =   0.1149034849
recursive, naive: J3(1) =  0.01956335398 
recursive, true:  J3(1) =  0.01956335398
recursive, naive: J4(1) = 0.002476638964 
recursive, true:  J4(1) = 0.002476638964
recursive, naive: J5(1) = 0.0002497577302 
recursive, true:  J5(1) = 0.0002497577302
recursive, naive: J6(1) = 2.093833802e-05 % <-- error becomes visible
recursive, true:  J6(1) =  2.0938338e-05
recursive, naive: J7(1) = 1.502326053e-06
recursive, true:  J7(1) = 1.502325817e-06
recursive, naive: J8(1) = 9.422672065e-08 
recursive, true:  J8(1) = 9.422344173e-08
recursive, naive: J9(1) = 5.301477313e-09 
recursive, true:  J9(1) = 5.24925018e-09
recursive, naive: J10(1) = 1.19987098e-09 
recursive, true:  J10(1) = 2.630615124e-10
...

[Richard Atacador]

Thank you very much Ced!

Given your insight I am definitely motivated to reconstruct my code to this problem with more thought and consideration.

I appreciate your time invested to help me out.

Cheers

Answer 3

[Jolanda Müller]

The issue is numerical, as already explained in Ced's answer. A solution to your problem is to start with precise values for the biggest desired N, and work your way backwards. This way, the relative error, stays below 10^-13 if for all nmax <= 142. [Sidenote: bessel(142,1) = 6.6430*10^-289 is the biggest n, for which the return value is non-zero. For nmax bigger than 142, besselj(>142,1) returns a true 0, thus breaking the possibility of getting the values for smaller n through backward recursion.]

nmax = 142;
bessel_val = zeros(nmax+1,1);
bessel_val(nmax)   = besselj(nmax-1,1);
bessel_val(nmax-1) = besselj(nmax-2,1);

for n = flip(1:nmax-2)
    bessel_val(n)   =  2*(n)*bessel_val(n+1) - bessel_val(n+2)    ;
end

relative_error = zeros(1,nmax);
for n = 1:nmax
    comp = bessel_val(n+1);
    realb = besselj(n,1);
    diff = comp-realb;
    relative_error(n)=abs(diff)/abs(realb);
end
semilogy(relative_error);