Bessel function of the first kind

3 次查看（过去 30 天）
Richard Atacador
Richard Atacador 2016-3-26
编辑： Jolanda Müller 2021-3-26
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.

请先登录，再回答此问题。

回答（3 个）

Torsten
Torsten 2016-3-29
For an explanation, google "Bessel's maze" and take the first hit.
Best wishes
Torsten.
Ced
Ced 2016-3-26
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
...
  1 个评论
Richard Atacador
Richard Atacador 2016-3-26
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

请先登录，再进行评论。

Jolanda Müller
Jolanda Müller 2021-3-26
编辑：Jolanda Müller 2021-3-26
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);

类别

在 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!

Translated by