Boost C++ Libraries

This example uses Boost.Multiprecision to implement a high-precision Gauss-Laguerre quadrature integration. The quadrature is used to calculate the airy_ai(x) function for real-valued x on the positive axis with x >= 1.

In this way, the integral representation could be seen as part of a scheme to calculate real-valued Airy functions on the positive axis for medium to large argument. A Taylor series or hypergeometric function (not part of this example) could be used for smaller arguments.

This example has been tested with decimal digits counts ranging from 21...301, by adjusting the parameter local::my_digits10 at compile time.

The quadrature integral representaion of airy_ai(x) used in this example can be found in:

A. Gil, J. Segura, N.M. Temme, "Numerical Methods for Special Functions" (SIAM Press 2007), ISBN 9780898717822, Sect. 5.3.3, in particular Eq. 5.110, page 145.

Subsequently, Gil et al's book cites the another work: W. Gautschi, "Computation of Bessel and Airy functions and of related Gaussian quadrature formulae", BIT, 42 (2002), pp. 110-118, https://doi.org/10.1023/A:1021974203359 that is also available as Gautschi_169.pdf.

This Gauss-Laguerre quadrature is designed for airy_ai(x) with real-valued x >= 1.

The example uses Gauss-Laguerre quadrature integration to compute airy_ai(x / 7) with 7 <= x <= 120 and where x is incremented in steps of 1.

During development of this example, we have empirically found the numbers of Gauss-Laguerre coefficients needed for convergence when using various counts of base-10 digits.

Let's calibrate, for instance, the number of coefficients needed at the point x = 1.

Empirical data were used with Wolfram Alpha :

Fit[{{21.0, 3.5}, {51.0, 11.1}, {101.0, 22.5}, {201.0, 46.8}}, {1, d, d^2}, d]FullSimplify[%]
  0.0000178915 d^2 + 0.235487 d - 1.28301
  or
  -1.28301 + (0.235487 + 0.0000178915 d) d

We need significantly more coefficients at smaller real values than are needed at larger real values because the slope derivative of airy_ai(x) gets more steep as x approaches zero. laguerre_order is calculated using this equation.

Snippets from (copious) output from a progress bar during calculation of approximate root estimates followed by calculation of accurate abscissa and weights is:

std::numeric_limits<local::float_type>::digits10: 101
laguerre_order: 2291

Finding the approximate roots...
root_estimates.size(): 1, 0.0%
root_estimates.size(): 8, 0.3%
root_estimates.size(): 16, 0.7%
...
root_estimates.size(): 2288, 99.9%
root_estimates.size(): 2291, 100.0%


Calculating abscissas and weights. Processed 1, 0.0%
Calculating abscissas and weights. Processed 9, 0.4%
...
Calculating abscissas and weights. Processed 2289, 99.9%
Calculating abscissas and weights. Processed 2291, 100.0%

Finally the result using Gauss-Laguerre quadrature is compared with a calculation using cyl_bessel_k, and both are listed, finally confirming that none differ more than a small tolerance.

airy_ai_value  : 0.13529241631288141552414742351546630617494414298833070600910205475763353480226572366348710990874867334
airy_ai_control: 0.13529241631288141552414742351546630617494414298833070600910205475763353480226572366348710990874868323
airy_ai_value  : 0.11392302126009621102904231059693500086750049240884734708541630001378825889924647699516200868335286103
airy_ai_control: 0.1139230212600962110290423105969350008675004924088473470854163000137882588992464769951620086833528582
...
airy_ai_value  : 3.8990420982303275013276114626640705170145070824317976771461533035231088620152288641360519429331427451e-22
airy_ai_control: 3.8990420982303275013276114626640705170145070824317976771461533035231088620152288641360519429331426481e-22

Total... result_is_ok: true

For more detail see comments in the source code for this example at gauss_laguerre_quadrature.cpp.