# Calculation of Bessel Functions

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

-
Perhaps you would be interested in an integral representation for the Bessel function? – Antonio Vargas May 6 '12 at 2:18
thanks for fixing up my LaTex code, you beat me to it! Would the integral representation get rid of the infinite sum? What i want to do is to get a function just in terms of $\beta$ so that I can manipulate it. preferably the sigma summation can be simplified to a constant or an approximation. – suzu May 6 '12 at 2:27
Sure thing! Any of those integral representations are equal to the sum you gave (where applicable, look to the right of the formulas on that page to see their domains). If you're interested in an approximation, say for large $\beta$, you may take as many or as few terms as you like from an asymptotic expansion for the Bessel function. I strongly urge you to take a look around the rest of that site's Bessel function category. – Antonio Vargas May 6 '12 at 2:38
The best algorithm to use very much depends on what kind of parameters you are interested in. For $\alpha$: is it an integer, or an arbitrary real number, or an arbitrary complex number? For $\beta$, is it small or large? To give an example: there are what are called asymptotic expansions you can use for $\alpha$ large, $\beta$ large, or both $\alpha$ and $\beta$ large. For other arguments, things are not too simple. – J. M. May 6 '12 at 4:25
Antonio, at first glance, it seems that Eq. 10.17.3 gives even more complicated equations than before. To answer J.M.'s question (and maybe for a little more detail to my problem), I am looking for something when $\alpha$ values of 1, 3, 5 and $\beta$ equal to 4. This is troubling as I cannot really use any of the limiting formulas referred to by Antonio. Please advise. – suzu May 6 '12 at 17:52

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.

-
Would Simpson's rule yield any improvement? – Antonio Vargas May 7 '12 at 1:16
@Antonio: Not much of an improvement in terms of accuracy, considering that the approximation using Simpson's rule is more complicated, and you are only limited to an even number of sub-intervals. Recall that Simpson's rule can be expressed as a linear combination of two trapezoidal rule evaluations, one with twice the number of sub-intervals than the other; since the trapezoidal rule already converges quickly, the improvement afforded by Simpson's is paltry. – J. M. May 7 '12 at 1:25
Any particular reason for using the midpoint rule instead of the trapezoidal rule (which I think of as the rectangle rule in the case of periodic integrands)? – Jitse Niesen May 7 '12 at 9:29
@Jitse: the weights are simpler in the midpoint case (i.e. I don't have to take half of the first and last terms), and the midpoint rule gives slightly more accuracy than the trapezoidal rule, at least in the tests I did. – J. M. May 7 '12 at 10:03
Are you sure your formula is correct? When I did some tests, it appears that $J_2$, $J_6$, $J_{10}$ etc. needs to change sign to match the original. Seems like a $(-1)^{2n}$ factor is needed somewhere. Odd ones work OK. Also, can you elaborate some more on how did you derive those formulas? – SasQ Oct 21 '13 at 7:45

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.

-

@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

-
I thought you said $\beta$ was modestly sized; I did give the caveat that the midpoint rule works nicely only for modestly-sized arguments. This really is information you should have given earlier. – J. M. May 9 '12 at 1:43
As I already said in the comments, one uses vastly different methods depending on the type of argument expected; for large arguments like what you have, one uses asymptotic expansions. I'd write about them, but now I'm not sure if you'll suddenly pull the rug from under me and tell me that you are dealing with a different set of parameters again. – J. M. May 9 '12 at 1:45
Lastly: if you'll look carefully, the last expression in your post is an inequality; it bounds the Bessel function, but is in no way, shape, or form an approximation. – J. M. May 9 '12 at 1:47
Sorry for the confusion. I think $\beta$ would still be of "modest" size for higher frequencies, since the value of $\omega_m$ can be large since we are doing a frequency sweep (ranges from 10 Hz to 100MHz.) So $\beta$ ranges from 0.04 to 4. I suppose that means the frequency range from 0 - 100Hz cannot be easily calculated. – suzu May 9 '12 at 2:10
I didn't notice the $j$ before; again, evaluating at complex values requires different strategy from evaluating at real values. Can you maybe give the full story this time before we proceed further? – J. M. May 9 '12 at 17:42

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.

-