1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Gauss-Laguerre quadrature</title> 5<link rel="stylesheet" href="../../../../multiprecision.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../../../../index.html" title="Chapter 1. Boost.Multiprecision"> 8<link rel="up" href="../fp_eg.html" title="Examples"> 9<link rel="prev" href="variable_precision.html" title="Variable-Precision Newton Evaluation"> 10<link rel="next" href="../../interval.html" title="Interval Number Types"> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../../../boost.png"></td> 15<td align="center"><a href="../../../../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<a accesskey="p" href="variable_precision.html"><img src="../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../fp_eg.html"><img src="../../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../../index.html"><img src="../../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="../../interval.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h5 class="title"> 27<a name="boost_multiprecision.tut.floats.fp_eg.gauss_lagerre_quadrature"></a><a class="link" href="gauss_lagerre_quadrature.html" title="Gauss-Laguerre quadrature">Gauss-Laguerre 28 quadrature</a> 29</h5></div></div></div> 30<p> 31 This example uses <a href="https://www.boost.org/doc/libs/release/libs/multiprecision/doc/index.html" target="_top">Boost.Multiprecision</a> 32 to implement a high-precision Gauss-Laguerre quadrature integration. 33 The quadrature is used to calculate the <code class="computeroutput"><span class="identifier">airy_ai</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span></code> function for real-valued <code class="computeroutput"><span class="identifier">x</span></code> on the positive axis with <code class="computeroutput"><span class="identifier">x</span> <span class="special">>=</span> <span class="number">1</span></code>. 34 </p> 35<p> 36 In this way, the integral representation could be seen as part of a scheme 37 to calculate real-valued Airy functions on the positive axis for medium 38 to large argument. A Taylor series or hypergeometric function (not part 39 of this example) could be used for smaller arguments. 40 </p> 41<p> 42 This example has been tested with decimal digits counts ranging from 43 21...301, by adjusting the parameter <code class="computeroutput"><span class="identifier">local</span><span class="special">::</span><span class="identifier">my_digits10</span></code> 44 at compile time. 45 </p> 46<p> 47 The quadrature integral representaion of <code class="computeroutput"><span class="identifier">airy_ai</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span></code> used in this example can be found in: 48 </p> 49<p> 50 A. Gil, J. Segura, N.M. Temme, "Numerical Methods for Special Functions" 51 (SIAM Press 2007), ISBN 9780898717822, Sect. 5.3.3, in particular Eq. 52 5.110, page 145. 53 </p> 54<p> 55 Subsequently, Gil et al's book cites the another work: W. Gautschi, "Computation 56 of Bessel and Airy functions and of related Gaussian quadrature formulae", 57 BIT, 42 (2002), pp. 110-118, <a href="https://doi.org/10.1023/A:1021974203359" target="_top">https://doi.org/10.1023/A:1021974203359</a> 58 that is also available as <a href="https://www.cs.purdue.edu/homes/wxg/selected_works/section_02/169.pdf" target="_top">Gautschi_169.pdf</a>. 59 </p> 60<p> 61 This Gauss-Laguerre quadrature is designed for <code class="computeroutput"><span class="identifier">airy_ai</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span></code> with real-valued <code class="computeroutput"><span class="identifier">x</span> 62 <span class="special">>=</span> <span class="number">1</span></code>. 63 </p> 64<p> 65 The example uses Gauss-Laguerre quadrature integration to compute <code class="computeroutput"><span class="identifier">airy_ai</span><span class="special">(</span><span class="identifier">x</span> <span class="special">/</span> <span class="number">7</span><span class="special">)</span></code> with 66 <code class="computeroutput"><span class="number">7</span> <span class="special"><=</span> 67 <span class="identifier">x</span> <span class="special"><=</span> 68 <span class="number">120</span></code> and where <code class="computeroutput"><span class="identifier">x</span></code> 69 is incremented in steps of 1. 70 </p> 71<p> 72 During development of this example, we have empirically found the numbers 73 of Gauss-Laguerre coefficients needed for convergence when using various 74 counts of base-10 digits. 75 </p> 76<p> 77 Let's calibrate, for instance, the number of coefficients needed at the 78 point <code class="computeroutput"><span class="identifier">x</span> <span class="special">=</span> 79 <span class="number">1</span></code>. 80 </p> 81<p> 82 Empirical data were used with <a href="http://www.wolframalpha.com/" target="_top">Wolfram 83 Alpha</a> : 84 </p> 85<pre class="programlisting"><span class="identifier">Fit</span><span class="special">[{{</span><span class="number">21.0</span><span class="special">,</span> <span class="number">3.5</span><span class="special">},</span> <span class="special">{</span><span class="number">51.0</span><span class="special">,</span> <span class="number">11.1</span><span class="special">},</span> <span class="special">{</span><span class="number">101.0</span><span class="special">,</span> <span class="number">22.5</span><span class="special">},</span> <span class="special">{</span><span class="number">201.0</span><span class="special">,</span> <span class="number">46.8</span><span class="special">}},</span> <span class="special">{</span><span class="number">1</span><span class="special">,</span> <span class="identifier">d</span><span class="special">,</span> <span class="identifier">d</span><span class="special">^</span><span class="number">2</span><span class="special">},</span> <span class="identifier">d</span><span class="special">]</span><span class="identifier">FullSimplify</span><span class="special">[%]</span> 86 <span class="number">0.0000178915</span> <span class="identifier">d</span><span class="special">^</span><span class="number">2</span> <span class="special">+</span> <span class="number">0.235487</span> <span class="identifier">d</span> <span class="special">-</span> <span class="number">1.28301</span> 87 <span class="keyword">or</span> 88 <span class="special">-</span><span class="number">1.28301</span> <span class="special">+</span> <span class="special">(</span><span class="number">0.235487</span> <span class="special">+</span> <span class="number">0.0000178915</span> <span class="identifier">d</span><span class="special">)</span> <span class="identifier">d</span> 89</pre> 90<p> 91 We need significantly more coefficients at smaller real values than are 92 needed at larger real values because the slope derivative of <code class="computeroutput"><span class="identifier">airy_ai</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span></code> 93 gets more steep as x approaches zero. <code class="computeroutput"><span class="identifier">laguerre_order</span></code> 94 is calculated using this equation. 95 </p> 96<p> 97 Snippets from (copious) output from a progress bar during calculation 98 of approximate root estimates followed by calculation of accurate abscissa 99 and weights is: 100 </p> 101<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">local</span><span class="special">::</span><span class="identifier">float_type</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">:</span> <span class="number">101</span> 102<span class="identifier">laguerre_order</span><span class="special">:</span> <span class="number">2291</span> 103 104<span class="identifier">Finding</span> <span class="identifier">the</span> <span class="identifier">approximate</span> <span class="identifier">roots</span><span class="special">...</span> 105<span class="identifier">root_estimates</span><span class="special">.</span><span class="identifier">size</span><span class="special">():</span> <span class="number">1</span><span class="special">,</span> <span class="number">0.0</span><span class="special">%</span> 106<span class="identifier">root_estimates</span><span class="special">.</span><span class="identifier">size</span><span class="special">():</span> <span class="number">8</span><span class="special">,</span> <span class="number">0.3</span><span class="special">%</span> 107<span class="identifier">root_estimates</span><span class="special">.</span><span class="identifier">size</span><span class="special">():</span> <span class="number">16</span><span class="special">,</span> <span class="number">0.7</span><span class="special">%</span> 108<span class="special">...</span> 109<span class="identifier">root_estimates</span><span class="special">.</span><span class="identifier">size</span><span class="special">():</span> <span class="number">2288</span><span class="special">,</span> <span class="number">99.9</span><span class="special">%</span> 110<span class="identifier">root_estimates</span><span class="special">.</span><span class="identifier">size</span><span class="special">():</span> <span class="number">2291</span><span class="special">,</span> <span class="number">100.0</span><span class="special">%</span> 111 112 113<span class="identifier">Calculating</span> <span class="identifier">abscissas</span> <span class="keyword">and</span> <span class="identifier">weights</span><span class="special">.</span> <span class="identifier">Processed</span> <span class="number">1</span><span class="special">,</span> <span class="number">0.0</span><span class="special">%</span> 114<span class="identifier">Calculating</span> <span class="identifier">abscissas</span> <span class="keyword">and</span> <span class="identifier">weights</span><span class="special">.</span> <span class="identifier">Processed</span> <span class="number">9</span><span class="special">,</span> <span class="number">0.4</span><span class="special">%</span> 115<span class="special">...</span> 116<span class="identifier">Calculating</span> <span class="identifier">abscissas</span> <span class="keyword">and</span> <span class="identifier">weights</span><span class="special">.</span> <span class="identifier">Processed</span> <span class="number">2289</span><span class="special">,</span> <span class="number">99.9</span><span class="special">%</span> 117<span class="identifier">Calculating</span> <span class="identifier">abscissas</span> <span class="keyword">and</span> <span class="identifier">weights</span><span class="special">.</span> <span class="identifier">Processed</span> <span class="number">2291</span><span class="special">,</span> <span class="number">100.0</span><span class="special">%</span> 118</pre> 119<p> 120 Finally the result using Gauss-Laguerre quadrature is compared with a 121 calculation using <code class="computeroutput"><span class="identifier">cyl_bessel_k</span></code>, 122 and both are listed, finally confirming that none differ more than a 123 small tolerance. 124 </p> 125<pre class="programlisting"><span class="identifier">airy_ai_value</span> <span class="special">:</span> <span class="number">0.13529241631288141552414742351546630617494414298833070600910205475763353480226572366348710990874867334</span> 126<span class="identifier">airy_ai_control</span><span class="special">:</span> <span class="number">0.13529241631288141552414742351546630617494414298833070600910205475763353480226572366348710990874868323</span> 127<span class="identifier">airy_ai_value</span> <span class="special">:</span> <span class="number">0.11392302126009621102904231059693500086750049240884734708541630001378825889924647699516200868335286103</span> 128<span class="identifier">airy_ai_control</span><span class="special">:</span> <span class="number">0.1139230212600962110290423105969350008675004924088473470854163000137882588992464769951620086833528582</span> 129<span class="special">...</span> 130<span class="identifier">airy_ai_value</span> <span class="special">:</span> <span class="number">3.8990420982303275013276114626640705170145070824317976771461533035231088620152288641360519429331427451e-22</span> 131<span class="identifier">airy_ai_control</span><span class="special">:</span> <span class="number">3.8990420982303275013276114626640705170145070824317976771461533035231088620152288641360519429331426481e-22</span> 132 133<span class="identifier">Total</span><span class="special">...</span> <span class="identifier">result_is_ok</span><span class="special">:</span> <span class="keyword">true</span> 134</pre> 135<p> 136 For more detail see comments in the source code for this example at 137 <a href="http://www.boost.org/doc/libs/release/libs/multiprecision/doc/html/../../example/gauss_laguerre_quadrature.cpp" target="_top">gauss_laguerre_quadrature.cpp</a>. 138 </p> 139</div> 140<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 141<td align="left"></td> 142<td align="right"><div class="copyright-footer">Copyright © 2002-2020 John 143 Maddock and Christopher Kormanyos<p> 144 Distributed under the Boost Software License, Version 1.0. (See accompanying 145 file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>) 146 </p> 147</div></td> 148</tr></table> 149<hr> 150<div class="spirit-nav"> 151<a accesskey="p" href="variable_precision.html"><img src="../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../fp_eg.html"><img src="../../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../../index.html"><img src="../../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="../../interval.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a> 152</div> 153</body> 154</html> 155
