Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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?

share|improve this question
2  
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. –  Arkamis 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
1  
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. –  Arkamis Sep 4 '12 at 21:42
show 2 more comments

1 Answer

up vote 3 down vote accepted
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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.