# Calculation of bessel function versus matlab solution

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?

-
Shouldn't it be for m=0:loop_var (I'm not good Matlab programmer, its just a hunch). – Fabian Sep 4 '12 at 21:15
@Fabian Good catch. I've fixed that. – suzu Sep 4 '12 at 21:24
@suzu I've addressed your edited question in my answer. – Emily Sep 4 '12 at 21:29
@EdGorcenski Thanks to you as well! I've now edited my question a bit with a follow-up that I hope has a quick fix as well :) – suzu Sep 4 '12 at 21:32
Since you have a $\beta^m$ term in there (actually $\beta^{2m+\alpha}$, but we can ignore the factor of 2 and the $\alpha$), then what happens as $m$ becomes large is that $\beta^m$ blows up for $\beta > 1$. There is no simple fix for this. The easiest solution might be to solve the Bessel's differential equation with those parameters, using an appropriate ODE method. – Emily Sep 4 '12 at 21:42

## 1 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.

-