2
$\begingroup$

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?

  • 2
    Shouldn't it be for `m=0:loop_var` (I'm not good Matlab programmer, its just a hunch).2012-09-04
  • 0
    @Fabian Good catch. I've fixed that.2012-09-04
  • 0
    @suzu I've addressed your edited question in my answer.2012-09-04
  • 0
    @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 :)2012-09-04
  • 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.2012-09-04
  • 0
    @EdGorcenski In lieu of starting a new question, perhaps you can provide some comments to this (related) question. Since it appears simpler to use the besselj call in MATLAB, is there a nice way to force one of the input parameters to be a transfer function? I.E. for $besselj(nu, Z)$ Z is an integrator function $(1/s)$? Matlab only seems to want numeric input parameters, but I would like to know if there is a clever way around this. Thanks.2012-09-05
  • 0
    @suzu I think it would be best to make that another question. I'm not even sure what it would mean to call `besselj(nu,1/s)`. The function `besselj(nu,Z)` computes one of the solutions to Bessel's differential equation; I'm not sure it makes any sense at all to replace an independent variable with a transfer function in the analytic evaluation of an ODE.2012-09-05

1 Answers 1

5
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.