# bessel functions

Asked by aaa on 7 Jun 2012
Latest activity Answered by aaa on 11 Jun 2012

Hi,

I am trying to write a code to calculate the bessel functions of the first kind. The equation is given here: http://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind_:_J.CE.B1

I also have the following code:

```beta = 4;
alpha = 1;
iteration = 3;
```
```format long
```
```for m = 1:iteration
J(m) = (((-1)^m)/(factorial(m)*gamma(m+alpha+1)))*(beta/2)^(2*m+alpha);
end
```
```out = sum(J)
```

However, when i compare my results with the built in function besselj, I cannot get the same result. Unfortunately, I cannot see where I am going wrong.

All assistance is appreciated.

tia

## Products

No products are associated with this question.

Answer by aaa on 11 Jun 2012

I have solved this issue. I should have started the iterations at 0 instead of 1.

Answer by Honglei Chen on 8 Jun 2012

The iteration should go to infinity but yours stop at 3.

In fact, this kind of implementation has another numerical issue, see the notes below from factorial's reference page

Since double precision numbers only have about 15 digits, the answer is only accurate for n <= 21. For larger n, the answer will have the right magnitude, and is accurate for the first 15 digits.

aaa on 8 Jun 2012

hmm, i probably should have not left it at 3. However, I found that after any number greater than 8 for the number of iterations, produces the same result, so i'm not sure if an infinite sum is needed for an approximation.

Are you saying that for iteration < 21, the answer should still be accurate for the first 15 digits? for my needs, this is more than enough, but i found that even for iteration = 10, that the answer is significantly different than when i use the besselj( ) function call.

Honglei Chen on 8 Jun 2012

This is indeed the equation but I'm sure there are programming tricks to get around the numerical issue. I personally don't know what they are, but looks to me the brute force approach will hit the wall very soon due to numerical issues.

Answer by aaa on 8 Jun 2012

I want to also add that i tried to use the

```doc besselj
```

and then formula that was presented there. However, I still do not get the same result as besselj(), but I do get the same result as the formula i initially used.