Calculation of Bessel Functions

Question

I want to calculate the Bessel function, given by

$$J_\alpha (\beta) = \sum_{m=0}^{\infty}\frac{(-1)^m}{m!\Gamma(m+\alpha +1)} \left(\frac{\beta}{2}\right)^{2m}$$

I know there are some tables that exist for this, but I want to keep the $\beta$ variable (i.e. I want a symbolic form in terms of $\beta$). If there is a way to simplify the summation part of the equation and leave an equation only in terms of $\beta$, that would be very helpful. (I see there is a dependence on $2m$, but I would like to see a way to break down the "other half" of the equation.)

Another question I have is: how is this calculated for $\beta$ values that are greater than $1$? It seems to me that this would give an infinite sum.

I am looking for something for $\alpha=1,3,5$ and $\beta=4$.

Thanks in advance.

Answer 1

As I alluded to in the comments, in general one would have to write a book chapter's worth of paragraphs to talk about the evaluation of Bessel functions for various argument ranges. Here, things are easier, since I only have to deal with integer orders of modest size. I shall now demonstrate one of my favorite methods, due to Yudell Luke.

Our starting point here is the pair of integrals

$$J_n(x)=\begin{cases}\frac2{\pi}\int_0^{\pi/2}\cos(x\cos\,u)\cos\,nu\;\mathrm du&n\text{ even}\\\frac2{\pi}\int_0^{\pi/2}\sin(x\sin\,u)\sin\,nu\;\mathrm du&n\text{ odd}\end{cases}$$

Two very useful methods for numerically evaluating these integrals are the trapezoidal rule and the midpoint rule. In a sense, these two are very accurate methods for the job, thanks to the Euler-Maclaurin formula. (See this for a deeper discussion.)

Using the odd order case as a concrete example, there is the following approximation which uses the (sadly lesser-known) midpoint rule:

$$J_n(x)\approx\frac1{m}\sum_{k=0}^{m-1}\sin\left(x\sin\left(\frac{\pi}{2m}\left(k+\frac12\right)\right)\right)\sin\left(\frac{\pi n}{2m}\left(k+\frac12\right)\right)$$

where $m$ is an appropriately chosen integer. For the particular case described in your question, taking $m=8$ gives approximations good to at least ten digits. Increase $m$ as needed.

In the case of even $n$, just replace all sines with cosines.

Again, this method is only suitable for modest integer values of $n$ and modest values of $x$; other methods might be more accurate, more efficient, or both for other argument ranges.

Answer 2

Here's another method you might want to consider, if you're in the business of generating a sequence of Bessel functions of fixed argument and consecutive integer orders. To describe the algorithm, which is due to J.C.P. Miller, we first take for granted the inequality ($x$ here is assumed real)

$$|J_n(x)|\leq\frac1{n!}\left|\frac{x}{2}\right|^n$$

and the series

$$1=J_0(x)+2\sum_{k=1}^\infty J_{2k}(x)$$

as well as the recurrence relation

$$\frac{J_n(x)}{J_{n-1}(x)}=\frac{x}{2n-x\frac{J_{n+1}(x)}{J_n(x)}}$$

Miller's idea is to first use an estimate like the inequality I gave to reckon an integer $n^\ast$ such that $\frac{J_{n^\ast}(x)}{J_{n^\ast-1}(x)}$ is smaller than machine epsilon. Having done so, pick some arbitrary value as a starting point (essentially, $f\,J_{n^\ast}(x)$ with $f$ an unknown constant), apply the recurrence backwards an appropriate number of times, while accumulating an unnormalized sum ($f\,J_0(x)+f\,J_1(x)+\cdots$). Once you've ended your recurrence, you can use the sum to normalize the recurrence values you stored along the way, which yields the Bessel function values you need.

To be more concrete, I shall present a Mathematica implementation of Miller's algorithm (which should be easily translatable to your favorite computing environment). I have chosen $n^\ast=24$ here; using the inequality with $x=4$, we have $|J_{24}(4)|\leq\frac{(4/2)^24}{24!}\approx 2.7\times10^{-17}$

x = N[4, 20];
n = 24;
(*hl accumulates ratios of Bessels h; s is the unnormalized sum*)
h = 0; s = 0; hl = {};
Do[
  h = x/(2 k - x h); (*recurrence relation*)
  hl = {h, hl};
  s = h (s + 2 Boole[EvenQ[k]]); , (*i.e., add 2 if k is even, else 0*)
  {k, n - 1, 1, -1}];
hl = Flatten[{1/(1 + s), hl}]; (*numerator is the value of the series*)
Do[hl[[k]] *= hl[[k - 1]], {k, 2, Length[hl]}];
hl

After executing the snippet, hl holds approximations to $J_0(4),J_1(4),\dots,J_{23}(4)$. When I tested it out, the first nineteen values generated were good to at least ten digits. Adapt the algorithm as needed.

Answer 3

@J.M. I really appreciate you taking the time to share your knowledge in this thread. But I think I'd like to backtrack a little bit here.

I think for my purposes, this equation that you provided is most useful.

$$J_n(x)\approx\frac1{m}\sum_{k=0}^{m-1}\sin\left(x\sin\left(\frac{\pi}{2m}\left(k+\frac12\right)\right)\right)\sin\left(\frac{\pi n}{2m}\left(k+\frac12\right)\right)$$

At the end of the day, I think what I need to know is how $\beta$ affects the $J_o(\beta)$. As you probably know, $\beta$ (in communication theory) is equal to $\frac{\Delta \omega}{\omega_m}$ such that $\Delta \omega$ is constant and $\omega_m$ is the modulation frequency. For my case, I have also defined $s = j\omega_m$. (This is important for later)

I have also taken your suggestion and used m = 8. Therefore, we can simplify the above equation to be

$$J_n(\frac{\Delta \omega}{\omega_m})\approx\frac1{8}\sum_{k=0}^{7}\sin\left(\frac{\Delta \omega}{\omega_m}\sin\left(\frac{\pi}{16}\left(k+\frac12\right)\right)\right)\sin\left(\frac{\pi n}{16}\left(k+\frac12\right)\right)$$

Essentially, since I will be doing a frequency response later, I want to change the $\omega_m$ term to be in terms of $s$. So, we can do the substitution for $\omega_m = -js$. So if, we also choose n = 2 (as an example, we are left with.

$$J_2(\frac{\Delta \omega}{\omega_m})\approx\frac1{8}\sum_{k=0}^{7}\sin\left(\frac{\Delta \omega}{-js}\sin\left(\frac{\pi}{16}\left(k+\frac12\right)\right)\right)\sin\left(\frac{2\pi }{16}\left(k+\frac12\right)\right)$$

The problem I am seeing is that for large values of $\Delta\omega$ (around 4k for this case), and putting that into a cosine, $J_n(\beta)$ returns infinite.

Is there something obvious that I am missing here? It doesn't make sense that I cannot evaluate for this size of $\beta$.

I also see in your most recent post that this equation

$$|J_n(x)|\leq\frac1{n!}\left|\frac{x}{2}\right|^n$$

This actually looks much simpler to use and since $\beta$ is positive, I can ignore the absolute values on the right hand side, and would also not have to deal with the cosine/sines like above. Do you think this is possible?

Sorry in advance for the long read

Answer 4

The overall objective of this work is to create a bode plot/frequency response of $J_n(\beta)*H(j(\omega_{so}+n\omega_m)$. Thus, I have made the definition $s = j\omega_m$, since $\omega_{so}$ is a bias point. We also know that $\beta = \frac{\Delta \omega}{\omega_m}$. Since we want a frequency response, I think it makes sense to rewrite $\beta$ in terms of $s$. Therefore, we are left with $\beta = \frac{\Delta \omega}{-js}$. I also know (from previous calculation) that the value of $\Delta \omega$ is equation to 4120. Right away, I can see a problem might occur at low frequency, as the value of $\beta$ will be large.

The worst case would be to have a $\omega_m$ = 1Hz, but this gives a $\beta$ value in the range of 4000. So, my thinking was to start at $\omega_m$ = 1kHz, which is not the worst case, but it still gives a reasonable number of sidebands to work with. So, if we do a frequency sweep from 1kHz to 100MHz, we should be able to capture a reasonable plot of the frequency response. Although to me, this is difficult now as since we are varying $\omega_m$, so the value of $\beta$ is changing, and thus $J_n(\beta)$ is also changing.

Initially this wasn't a problem as I had assumed that $\beta$ was very small, such that I could approximate $J_n(\beta)$ as simply $\frac{\Delta \omega}{2\omega_m}$, and I could easily do the frequency sweep of this. However, my results did not end up being very accurate, since $\beta$ was not actually small. I want to do the same thing, but don't have a simple expression for $J_n(\beta)$ anymore.

So, I think what I am looking for is an expression for $J_n(\beta)$ for changing values of $\beta$?

I think I can still use

$$J_n(x)\approx\frac1{m}\sum_{k=0}^{m-1}\sin\left(x\sin\left(\frac{\pi}{2m}\left(k+\frac12\right)\right)\right)\sin\left(\frac{\pi n}{2m}\left(k+\frac12\right)\right)$$

as a starting point, but I ran into problems evaluating the sine and cosines, and was not able to easily bypass this.

I would like to have values of $n$ for 1 - 5.

Hopefully this is more clear, please let me know if you have any questions and I will try to reply to them ASAP.

Thanks.