Calculation of bessel function versus matlab solution

6.6k Views Asked by At

I am looking to calculate the Bessel function of the first kind $J_o(\beta)$. I am using the formula (referenced from wikipedia) to accomplish this.
$$J_\alpha (\beta) = \sum_{m=0}^{\infty}\frac{(-1)^m}{m!\Gamma(m+\alpha +1)} \left(\frac{\beta}{2}\right)^{2m+\alpha}$$

I am also aware of that MATLAB has a function which can calculate the solution as well. The function call is $besselj(nu,Z)$. However, what I am finding is that there is a discrepancy between what my $for$ loop calculates and what MATLAB outputs.

Does anyone see why this might be happening? I have included my code for reference:

loop_var = 100;
beta = 0.120;
alpha = 0;

sum = 0;
amp = 0;

[matl,ierr]  =besselj(alpha,beta)

for m=0:loop_var
    amp=amp+((-1)^m)/(factorial(m)*gamma(m+alpha+1))*(beta/2)^(2*m+alpha);
end
amp

Thanks for all comments and suggestions.

EDIT: Thanks to Fabian and Ed for catching that error.

I suppose as a "followup" question, the wikipedia formula and that MATLAB function seem to only match for "small" values of $\beta$. After $\beta$ is great than approximately 10, there is some error between the two values. Does anyone then know how the Matlab number is obtained?

1

There are 1 best solutions below

0
On BEST ANSWER
for m=loop_var
    amp=amp+((-1)^m)/(factorial(m)*gamma(m+alpha+1))*(beta/2)^(2*m+alpha);
end

This code only executes the loop once.

for m=0:loop_var
    amp=amp+((-1)^m)/(factorial(m)*gamma(m+alpha+1))*(beta/2)^(2*m+alpha);
end

This code is correct (notice the 0: after m=).

For your edited question...

From the MATLAB besselj documentation:

Algorithms

The besselj function uses a Fortran MEX-file to call a library developed by D.E. Amos [3] [4].

[3] Amos, D.E., "A Subroutine Package for Bessel Functions of a Complex Argument and Nonnegative Order," Sandia National Laboratory Report, SAND85-1018, May, 1985.

[4] Amos, D.E., "A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order," Trans. Math. Software, 1986.

The first paper can be found here.