...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
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.