Explanation of continued fraction for Bessel functions

bessel functionscontinued-fractionsreference-request

Doing some computational searches, I found some nice continued fractions that could be used to compute Bessel functions of the first kind (I and J):

-1 + besseli(0, 2)                 4           
────────────────── = 5 - ──────────────────────
-2 + besseli(0, 2)                    9        
                         10 - ─────────────────
                                        16     
                              17 - ────────────
                                           25  
                                   26 - ───────
                                             36
                                        37 - ──
                                             ..
⎛        2                      2          ⎞
⎝a(n) = n  + 4⋅n + 5, b(n) = - n  - 2⋅n - 1⎠

===============================================

  besseli(1, 2)                    2          
────────────────── = 3 - ─────────────────────
-1 + besseli(1, 2)                   6        
                         7 - ─────────────────
                                       12     
                             13 - ────────────
                                          20  
                                  21 - ───────
                                            30
                                       31 - ──
                                            ..
⎛        2                      2    ⎞
⎝a(n) = n  + 3⋅n + 3, b(n) = - n  - n⎠

==============================================

 -besselj(1, 2)                    2          
────────────────── = 1 + ─────────────────────
-1 + besselj(1, 2)                   6        
                         5 + ─────────────────
                                       12     
                             11 + ────────────
                                          20  
                                  19 + ───────
                                            30
                                       29 + ──
                                            ..
⎛        2                    2    ⎞
⎝a(n) = n  + 3⋅n + 1, b(n) = n  + n⎠

==============================================

 -2⋅besselj(0, 2)                    3           
──────────────────── = 4 + ──────────────────────
-1 + 4⋅besselj(0, 2)                    8        
                           10 + ─────────────────
                                          15     
                                18 + ────────────
                                             24  
                                     28 + ───────
                                               35
                                          40 + ──
                                               ..
⎛        2                    2      ⎞
⎝a(n) = n  + 5⋅n + 4, b(n) = n  + 2⋅n⎠

and others. The left sides are all Mobius transformations of Bessel functions, so that transformation can be inverted to give one extra leading term to the continued fraction. These seem like they might even be useful computationally.

I am wondering if these are known — I would figure they must be! — but all the text I could find about Bessel function GCFs involves only ratios of Bessel functions. For instance, at the bottom https://mathworld.wolfram.com/ContinuedFractionConstants.html it discusses how CFs with terms in arithmetic progress can be solved using ratios of Bessel functions. There's no reference in https://dlmf.nist.gov/10.10 or https://dlmf.nist.gov/10.33 (which would be the right places for it to be listed in DLMF). And https://aip.scitation.org/doi/pdf/10.1063/1.168382 discusses the ratio aspect of the J's but laments the recursive nature of computing them.

Questions:

  1. Are these kinds of expressions known to wider literature?
  2. Can someone provide a proof of their validity? (I have verified them out to ~1000 digits so I have very high confidence, but no proof)
  3. How can I get better at literature searching so I could find sources on these expressions myself, in the future? 😛

Best Answer

Figured this out. The Bessel J function has the Talyor series

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

Call the terms in this sum $t_0, t_1$ etc. The first term $t_0$ is $\frac{1}{\Gamma(\alpha+1)}{\left(\frac{x}{2}\right)}^{2\alpha}$, and the ratio between successive terms is

$$r_m = \frac{t_m}{t_{m-1}} = -\frac{1}{m(m+\alpha)}\frac{x^2}{4}$$

Apply Euler's continued fraction formula, which says that in a general a sum

$$t_0 + t_0r_1 + t_0r_1r_2 + \dots = \cfrac{t_0}{1 - \cfrac{r_1}{1 + r_1 - \cfrac{r_2}{1 + r_2 - \cfrac{r_3}{\ddots}}}}\,$$

With our values, this becomes

$$J_\alpha(x) = \cfrac{\frac{1}{\Gamma(\alpha+1)}{\left(\frac{x}{2}\right)}^{2\alpha}}{1 - \cfrac{\frac{-1}{1(1+\alpha)}\frac{x^2}{4}}{1 - \frac{1}{1(1+\alpha)}\frac{x^2}{4} - \cfrac{\frac{-1}{2(2+\alpha)}\frac{x^2}{4}}{1 - \frac{1}{2(2+\alpha)}\frac{x^2}{4} - \cfrac{\frac{-1}{3(3+\alpha)}\frac{x^2}{4}}{\ddots}}}}$$

Setting $x=2$, move the prefactor to the other side, and apply an equivalence transformation where multiply through by the polynomial $n(n+\alpha)$:

$$J_\alpha(2)\Gamma(\alpha+1) = \cfrac{1}{1 - \cfrac{\frac{-1(1+\alpha)}{1(1+\alpha)}}{1(1+\alpha) - \cfrac{\frac{-1(1+\alpha)2(2+\alpha)}{2(2+\alpha)}}{2(2+\alpha)\left(1 - \frac{1}{2(2+\alpha)}\right) - \cfrac{\frac{-2(2+\alpha)3(3+\alpha)}{3(3+\alpha)}}{\ddots}}}}$$

Simplify the fractions:

$$J_\alpha(2)\Gamma(\alpha+1) = \cfrac{1}{1 - \cfrac{-1}{1(1+\alpha) - 1 - \cfrac{-1(1+\alpha)}{2(2+\alpha) - 1 - \cfrac{-2(2+\alpha)}{\ddots}}}}$$

Take $-1/(1-1/x)$ of each side:

$$\frac{-1}{1-\frac{1}{J_\alpha(2)\Gamma(\alpha+1)}} = 1(1+\alpha) - 1 - \cfrac{-1(1+\alpha)}{2(2+\alpha) - 1 - \cfrac{-2(2+\alpha)}{3(3+\alpha) - 1 \ddots}}$$

Now put the left side as a simple ratio involving $J$, and rewrite the right side to look like polynomials:

$$\frac{-J_\alpha(2)\Gamma(\alpha+1)}{J_\alpha(2)\Gamma(\alpha+1)-1} = (1^2+\alpha 1-1) + \cfrac{1^2+\alpha 1}{(2^2+\alpha 2-1) + \cfrac{2^2+\alpha 2}{(3^2+\alpha 3-1) + \ddots}}$$

Which explains the results in my initial post. For example with $\alpha=1$, this is $$\frac{-J_1(2)}{J_1(2)-1} = a(0) + \cfrac{b(1)}{a(1) + \cfrac{b(2)}{a(2) + \ddots}}$$ with $a(n) = n^2 + 3n + 1$, $b(n) = n^2 + n$, which is the third result above.


If you wanted to use this to evaluate $J$ at integer $x$ and $\alpha$ with only integers in the continued fraction, the relevant expression would be

$$J_\alpha(x) = \frac{1}{\Gamma(\alpha+1)}{\left(\frac{x}{2}\right)}^{2\alpha}\cdot \cfrac{1}{1 + \cfrac{x^2}{4\cdot 1(1+\alpha) - x^2 + \cfrac{4\cdot 1(1+\alpha)x^2}{4\cdot 2(2+\alpha) - x^2 + \cfrac{4\cdot 2(2+\alpha)x^2}{4\cdot 3(3+\alpha) - x^2 \ddots}}}}$$ $$ = \frac{1}{\Gamma(\alpha+1)}{\left(\frac{x}{2}\right)}^{2\alpha}\cdot \cfrac{1}{1 + \cfrac{x^2}{b(0) + \cfrac{a(1)}{b(1) + \cfrac{a(2)}{b(2) +\ddots}}}}$$ $$\textrm{with }b(n) = 4n^2 + (8+4\alpha)n + (4+4\alpha-x^2),\quad a(n) = 4n^2 + 4\alpha n - x^2$$ (This switches $a$ and $b$ to match Wikipedia, not the notation in my original post.)

I ran this in Python and, compared to the mpmath implementation of besselj, which uses essentially the same hypergeometric formula, this runs about ten times faster for integer arguments, by using only integer arithmetic.

import mpmath

def gcf(a_, b_):
    pairs = zip(a_, b_)
    prev_p, prev_q = 1, 0
    p, q = next(pairs)[1], 1
    for an, bn in pairs:
        tmp_a = p
        tmp_b = q
        p = bn * p + an * prev_p
        q = bn * q + an * prev_q
        prev_p = tmp_a
        prev_q = tmp_b
    return mpmath.mpf(p)/q

def gcfBessel(a, x, N=20):
    x2 = x*x
    a_ = [None, x2] + [(4*i**2 + 4*a*i) * x2 for i in range(1,N)]
    b_ = [1] + [4*i**2 + (8+4*a)*i + (4+4*a-x2) for i in range(N)]
    gcfPart = gcf(a_, b_)
    return x**a / (2**a * mpmath.gamma(a+1) * gcfPart)

mpmath.besselj(3.4, 2)
gcfBessel(3.4, 2)

mpmath.besselj(3.4, 3)
gcfBessel(3.4, 3)

import timeit
mpmath.mp.dps = 10000

str(mpmath.besselj(5, 3)) == str(gcfBessel(5, 3, 1900)) #True -- it gets enough precision.

timeit.timeit(lambda: mpmath.besselj(5, 3), number=10) #1.11 seconds
timeit.timeit(lambda: gcfBessel(5, 3, 1900), number=10) #0.122 seconds
Related Question