5

For fixed positive $x$, I would like to accurately (close to full precision) evaluate the ratio of Bessel functions of the first kind $$ R_n(x):= \frac{J_{n+1}(x)}{J_n(x)} $$ as $n$ becomes extremely large. Evaluating the special functions directly fails because both numerator and denominator become very small. However, it is known that this ratio obeys the large $n$ asymptotics $$ R_n(x) \sim \sqrt{\frac{n}{n+1}} \frac{x}{2(n+1)} $$ so in principle it seems possible that there is a stable method for computing this quantity. Does anybody have any suggestions?

Christopher A. Wong
  • 22,445
  • 3
  • 51
  • 82
  • You can use continued fractions, see http://dlmf.nist.gov/10.10 – gammatester Oct 31 '17 at 22:35
  • Related: https://math.stackexchange.com/questions/1871733/convergence-of-a-harmonic-continued-fraction – Jack D'Aurizio Oct 31 '17 at 23:47
  • @gammatester Thank you for your suggestion. I do not understand how to cast the continued fraction in a form that is numerically stable. Even after a few iterations, the error blows up. – Christopher A. Wong Nov 02 '17 at 21:39
  • It cannot be (much) worse than your formula: Using the 10.10.2, first order approximation gives already $\frac{x}{2(1+n)}$ which basically is your asymptotics, the second order approximation is $$R_n(x) \sim \frac{x}{2(1+n)\left(1-\frac{x^2}{4(n+1)(n+2)}\right)}$$ – gammatester Nov 02 '17 at 22:01
  • @gammatester I see. I was doing something else. I was using the first expression to compute the recurrence $r_n = 2n/k - 1/r_{n-1}$, which is highly unstable. – Christopher A. Wong Nov 02 '17 at 22:20

2 Answers2

1

Using $$J_{\nu}(x) \approx \frac{1}{\sqrt{2 \, \pi \, \nu}} \, \left(\frac{e \, x}{2 \, \nu}\right)^{n}$$ for $\nu \to \infty$ leads to $$R_{n}(x) = \frac{J_{n+1}(x)}{J_{n}(x)} \approx \frac{e \, x}{2} \, \frac{1}{n+1} \, \left(\frac{n}{n+1}\right)^{n+ \frac{1}{2}}.$$

If this is not within the realm sought then continued fractions, approximation of integral representations, etc may be applied. Similar results may be obtained though.

Leucippus
  • 26,329
0

Boost.math uses Miller's algorithm for this problem; see J.M.'s description here. Code to compute the ratio to 100 digits or so:

#include <iostream>
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>

using std::sqrt;
using boost::math::cyl_bessel_j;
using namespace boost::multiprecision;

int main(int argc, char** argv)
{
    cpp_dec_float_100 x{5.7};
    for (int i = 10000; i < 50000; ++i)
    {
        std::cout << "J_" << i << "/J_" << i+1 << " = " << cyl_bessel_j<cpp_dec_float_100>(i+1, x)/cyl_bessel_j<cpp_dec_float_100>(i, x) << "\n";
        std::cout << "asymptotic form = " << sqrt((cpp_dec_float_100)i /(cpp_dec_float_100) i+1)*x/(2*i+2) << "\n";
    }
}
user14717
  • 4,902