Sign up ×
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 have a Bessel function of the first kind given by the equation $$J_\alpha (\beta) = \sum_{m=0}^{\infty}\frac{(-1)^m}{m!\Gamma(m+\alpha +1)} \left(\frac{\beta}{2}\right)^{2m+\alpha}$$

I am trying to write a simple matlab code calculate the values of $J_{\alpha}(\beta)$. So far I have the following

beta = 4;
alpha = 1;
iteration = 3;

format long

for m = 1:iteration
    J(m) = (((-1)^(m-1))/(factorial(m-1)*gamma(m+alpha)))*(beta/2)^(2*m+alpha-2);

out = sum(J) 

For some reason, when I compare my results to the besselj() function in MS Excel, I do not get the correct answer. I have tried also making the number of summing terms higher (~30), but I'd like to find the fewest number of summing terms as possible. Is there something I am missing here in my code?

I suppose my other question is this. Intuitively, I thought there was an error in this formula, as I expected J to be a large number (in fact infinite) for $\beta$ values larger than 1. Why does this infinite summation work?

I am also aware there are other "forms" of this equation, in the form of integrals and some approximations, but I think this form is the simplest for me to work with right now.

Thanks in advance to all posters.

UPDATE: I used the provided documentation in Matlab, and tried using the formula of their implementation. However, I actually get the same answer as my formula (nothing exciting there), but still not the same answer as the besselj() function call. I'm stumped!

share|cite|improve this question
Any reason why you don't want to use Matlab's built in implementation of Bessel functions? (See It is almost certainly both quicker and more accurate than what you could write yourself without working incredibly hard. – mrf Jun 7 '12 at 21:53
@mrf What you say is true and correct, but I'd like to write the whole summation as just a addition of a few terms if possible. eg from m = 0 to m = 3. Right now, I'm just trying to figure out how many terms would be considered enough to get a reasonable amount of accuracy, but I can't quite figure out what is wrong with the Matlab code that I have implemented above. To me, it seems like it should be possible to estimate the value of J with a finite number of terms, and hopefully a relatively small number of terms at that. – suzu Jun 7 '12 at 23:03
You forget to include the $m=0$ term in your sum. – mrf Jun 9 '12 at 20:59
@mrf, thanks. I just figured that out late last night and forgot to post here. It's the little things.... – suzu Jun 11 '12 at 4:42

Your Answer


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

Browse other questions tagged or ask your own question.