libs/math/doc/html/math_toolkit/double_exponential/de_caveats.html

   1 <html>
   2 <head>
   3 <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
   4 <title>Caveats</title>
   5 <link rel="stylesheet" href="../../math.css" type="text/css">
   6 <meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
   7 <link rel="home" href="../../index.html" title="Math Toolkit 2.11.0">
   8 <link rel="up" href="../double_exponential.html" title="Double-exponential quadrature">
   9 <link rel="prev" href="de_thread.html" title="Thread Safety">
  10 <link rel="next" href="de_refes.html" title="References">
  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="de_thread.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../double_exponential.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="de_refes.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
  24 </div>
  25 <div class="section">
  26 <div class="titlepage"><div><div><h3 class="title">
  27 <a name="math_toolkit.double_exponential.de_caveats"></a><a class="link" href="de_caveats.html" title="Caveats">Caveats</a>
  28 </h3></div></div></div>
  29 <p>
  30         A few things to keep in mind while using the tanh-sinh, exp-sinh, and sinh-sinh
  31         quadratures:
  32       </p>
  33 <p>
  34         These routines are <span class="bold"><strong>very</strong></span> aggressive about
  35         approaching the endpoint singularities. This allows lots of significant digits
  36         to be extracted, but also has another problem: Roundoff error can cause the
  37         function to be evaluated at the endpoints. A few ways to avoid this: Narrow
  38         up the bounds of integration to say, [a + &#949;, b - &#949;], make sure (a+b)/2 and
  39         (b-a)/2 are representable, and finally, if you think the compromise between
  40         accuracy an usability has gone too far in the direction of accuracy, file
  41         a ticket.
  42       </p>
  43 <p>
  44         Both exp-sinh and sinh-sinh quadratures evaluate the functions they are passed
  45         at <span class="bold"><strong>very</strong></span> large argument. You might understand
  46         that x<sup>12</sup>exp(-x) is should be zero when x<sup>12</sup> overflows, but IEEE floating point
  47         arithmetic does not. Hence <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">pow</span><span class="special">(</span><span class="identifier">x</span><span class="special">,</span> <span class="number">12</span><span class="special">)*</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">exp</span><span class="special">(-</span><span class="identifier">x</span><span class="special">)</span></code> is an indeterminate form whenever <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">pow</span><span class="special">(</span><span class="identifier">x</span><span class="special">,</span>
  48         <span class="number">12</span><span class="special">)</span></code>
  49         overflows. So make sure your functions have the correct limiting behavior;
  50         for example
  51       </p>
  52 <pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">{</span>
  53     <span class="keyword">double</span> <span class="identifier">t</span> <span class="special">=</span> <span class="identifier">exp</span><span class="special">(-</span><span class="identifier">x</span><span class="special">);</span>
  54     <span class="keyword">if</span><span class="special">(</span><span class="identifier">t</span> <span class="special">==</span> <span class="number">0</span><span class="special">)</span>
  55     <span class="special">{</span>
  56         <span class="keyword">return</span> <span class="number">0</span><span class="special">;</span>
  57     <span class="special">}</span>
  58     <span class="keyword">return</span> <span class="identifier">t</span><span class="special">*</span><span class="identifier">pow</span><span class="special">(</span><span class="identifier">x</span><span class="special">,</span> <span class="number">12</span><span class="special">);</span>
  59 <span class="special">};</span>
  60 </pre>
  61 <p>
  62         has the correct behavior for large <span class="emphasis"><em>x</em></span>, but <code class="computeroutput"><span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span>
  63         <span class="identifier">x</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="identifier">exp</span><span class="special">(-</span><span class="identifier">x</span><span class="special">)*</span><span class="identifier">pow</span><span class="special">(</span><span class="identifier">x</span><span class="special">,</span> <span class="number">12</span><span class="special">);</span> <span class="special">};</span></code> does
  64         not.
  65       </p>
  66 <p>
  67         Oscillatory integrals, such as the sinc integral, are poorly approximated
  68         by double-exponential quadrature. Fortunately the error estimates and L1
  69         norm are massive for these integrals, but nonetheless, oscillatory integrals
  70         require different techniques.
  71       </p>
  72 <p>
  73         A special mention should be made about integrating through zero: while our
  74         range adaptors preserve precision when one endpoint is zero, things get harder
  75         when the origin is neither in the center of the range, nor at an endpoint.
  76         Consider integrating:
  77       </p>
  78 <div class="blockquote"><blockquote class="blockquote"><p>
  79           <span class="serif_italic">1 / (1 +x^2)</span>
  80         </p></blockquote></div>
  81 <p>
  82         Over (a, &#8734;). As long as <code class="computeroutput"><span class="identifier">a</span> <span class="special">&gt;=</span> <span class="number">0</span></code> both
  83         the tanh_sinh and the exp_sinh integrators will handle this just fine: in
  84         fact they provide a rather efficient method for this kind of integral. However,
  85         if we have <code class="computeroutput"><span class="identifier">a</span> <span class="special">&lt;</span>
  86         <span class="number">0</span></code> then we are forced to adapt the range
  87         in a way that produces abscissa values near zero that have an absolute error
  88         of &#949;, and since all of the area of the integral is near zero both integrators
  89         thrash around trying to reach the target accuracy, but never actually get
  90         there for <code class="computeroutput"><span class="identifier">a</span> <span class="special">&lt;&lt;</span>
  91         <span class="number">0</span></code>. On the other hand, the simple expedient
  92         of breaking the integral into two domains: (a, 0) and (0, b) and integrating
  93         each seperately using the tanh-sinh integrator, works just fine.
  94       </p>
  95 <p>
  96         Finally, some endpoint singularities are too strong to be handled by <code class="computeroutput"><span class="identifier">tanh_sinh</span></code> or equivalent methods, for example
  97         consider integrating the function:
  98       </p>
  99 <pre class="programlisting"><span class="keyword">double</span> <span class="identifier">p</span> <span class="special">=</span> <span class="identifier">some_value</span><span class="special">;</span>
 100 <span class="identifier">tanh_sinh</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;</span> <span class="identifier">integrator</span><span class="special">;</span>
 101 <span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[&amp;](</span><span class="keyword">double</span> <span class="identifier">x</span><span class="special">){</span> <span class="keyword">return</span> <span class="identifier">pow</span><span class="special">(</span><span class="identifier">tan</span><span class="special">(</span><span class="identifier">x</span><span class="special">),</span> <span class="identifier">p</span><span class="special">);</span> <span class="special">};</span>
 102 <span class="keyword">auto</span> <span class="identifier">Q</span> <span class="special">=</span> <span class="identifier">integrator</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0</span><span class="special">,</span> <span class="identifier">constants</span><span class="special">::</span><span class="identifier">half_pi</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;());</span>
 103 </pre>
 104 <p>
 105         The first problem with this function, is that the singularity is at &#960;/2, so
 106         if we're integrating over (0, &#960;/2) then we can never approach closer to the
 107         singularity than &#949;, and for p less than but close to 1, we need to get <span class="emphasis"><em>very</em></span>
 108         close to the singularity to find all the area under the function. If we recall
 109         the identity <code class="literal">tan(&#960;/2 - x) == 1/tan(x)</code> then we can rewrite
 110         the function like this:
 111       </p>
 112 <pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[&amp;](</span><span class="keyword">double</span> <span class="identifier">x</span><span class="special">){</span> <span class="keyword">return</span> <span class="identifier">pow</span><span class="special">(</span><span class="identifier">tan</span><span class="special">(</span><span class="identifier">x</span><span class="special">),</span> <span class="special">-</span><span class="identifier">p</span><span class="special">);</span> <span class="special">};</span>
 113 </pre>
 114 <p>
 115         And now the singularity is at the origin and we can get much closer to it
 116         when evaluating the integral: all we have done is swap the integral endpoints
 117         over.
 118       </p>
 119 <p>
 120         This actually works just fine for p &lt; 0.95, but after that the <code class="computeroutput"><span class="identifier">tanh_sinh</span></code> integrator starts thrashing around
 121         and is unable to converge on the integral. The problem is actually a lack
 122         of exponent range: if we simply swap type double for something with a greater
 123         exponent range (an 80-bit long double or a quad precision type), then we
 124         can get to at least p = 0.99. If we want to go beyond that, or stick with
 125         type double, then we have to get smart.
 126       </p>
 127 <p>
 128         The easiest method is to notice that for small x, then <code class="literal">tan(x) &#8773; x</code>,
 129         and so we are simply integrating x<sup>-p</sup>. Therefore we can use this approximation
 130         over (0, small), and integrate numerically from (small, &#960;/2), using &#949; as a suitable
 131         crossover point seems sensible:
 132       </p>
 133 <pre class="programlisting"><span class="keyword">double</span> <span class="identifier">p</span> <span class="special">=</span> <span class="identifier">some_value</span><span class="special">;</span>
 134 <span class="keyword">double</span> <span class="identifier">crossover</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</span><span class="identifier">epsilon</span><span class="special">();</span>
 135 <span class="identifier">tanh_sinh</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;</span> <span class="identifier">integrator</span><span class="special">;</span>
 136 <span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[&amp;](</span><span class="keyword">double</span> <span class="identifier">x</span><span class="special">){</span> <span class="keyword">return</span> <span class="identifier">pow</span><span class="special">(</span><span class="identifier">tan</span><span class="special">(</span><span class="identifier">x</span><span class="special">),</span> <span class="identifier">p</span><span class="special">);</span> <span class="special">};</span>
 137 <span class="keyword">auto</span> <span class="identifier">Q</span> <span class="special">=</span> <span class="identifier">integrator</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="identifier">crossover</span><span class="special">,</span> <span class="identifier">constants</span><span class="special">::</span><span class="identifier">half_pi</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;())</span> <span class="special">+</span> <span class="identifier">pow</span><span class="special">(</span><span class="identifier">crossover</span><span class="special">,</span> <span class="number">1</span> <span class="special">-</span> <span class="identifier">p</span><span class="special">)</span> <span class="special">/</span> <span class="special">(</span><span class="number">1</span> <span class="special">-</span> <span class="identifier">p</span><span class="special">);</span>
 138 </pre>
 139 <p>
 140         There is an alternative, more complex method, which is applicable when we
 141         are dealing with expressions which can be simplified by evaluating by logs.
 142         Let's suppose that as in this case, all the area under the graph is infinitely
 143         close to zero, now inagine that we could expand that region out over a much
 144         larger range of abscissa values: that's exactly what happens if we perform
 145         argument substitution, replacing <code class="computeroutput"><span class="identifier">x</span></code>
 146         by <code class="computeroutput"><span class="identifier">exp</span><span class="special">(-</span><span class="identifier">x</span><span class="special">)</span></code> (note
 147         that we must also multiply by the derivative of <code class="computeroutput"><span class="identifier">exp</span><span class="special">(-</span><span class="identifier">x</span><span class="special">)</span></code>).
 148         Now the singularity at zero is moved to +&#8734;, and the &#960;/2 bound to -log(&#960;/2).
 149         Initially our argument substituted function looks like:
 150       </p>
 151 <pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[&amp;](</span><span class="keyword">double</span> <span class="identifier">x</span><span class="special">){</span> <span class="keyword">return</span> <span class="identifier">exp</span><span class="special">(-</span><span class="identifier">x</span><span class="special">)</span> <span class="special">*</span> <span class="identifier">pow</span><span class="special">(</span><span class="identifier">tan</span><span class="special">(</span><span class="identifier">exp</span><span class="special">(-</span><span class="identifier">x</span><span class="special">)),</span> <span class="special">-</span><span class="identifier">p</span><span class="special">);</span> <span class="special">};</span>
 152 </pre>
 153 <p>
 154         Which is hardly any better, as we still run out of exponent range just as
 155         before. However, if we replace <code class="computeroutput"><span class="identifier">tan</span><span class="special">(</span><span class="identifier">exp</span><span class="special">(-</span><span class="identifier">x</span><span class="special">))</span></code> by
 156         <code class="computeroutput"><span class="identifier">exp</span><span class="special">(-</span><span class="identifier">x</span><span class="special">)</span></code> for
 157         suitably small <code class="computeroutput"><span class="identifier">exp</span><span class="special">(-</span><span class="identifier">x</span><span class="special">)</span></code>, and
 158         therefore <code class="literal">x &gt; -log(&#949;)</code>, we can greatly simplify the expression
 159         and evaluate by logs:
 160       </p>
 161 <pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[&amp;](</span><span class="keyword">double</span> <span class="identifier">x</span><span class="special">)</span>
 162 <span class="special">{</span>
 163    <span class="keyword">static</span> <span class="keyword">const</span> <span class="keyword">double</span> <span class="identifier">crossover</span> <span class="special">=</span> <span class="special">-</span><span class="identifier">log</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</span><span class="identifier">epsilon</span><span class="special">());</span>
 164    <span class="keyword">return</span> <span class="identifier">x</span> <span class="special">&gt;</span> <span class="identifier">crossover</span> <span class="special">?</span> <span class="identifier">exp</span><span class="special">((</span><span class="identifier">p</span> <span class="special">-</span> <span class="number">1</span><span class="special">)</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">:</span> <span class="identifier">exp</span><span class="special">(-</span><span class="identifier">x</span><span class="special">)</span> <span class="special">*</span> <span class="identifier">pow</span><span class="special">(</span><span class="identifier">tan</span><span class="special">(</span><span class="identifier">exp</span><span class="special">(-</span><span class="identifier">x</span><span class="special">)),</span> <span class="special">-</span><span class="identifier">p</span><span class="special">);</span>
 165 <span class="special">};</span>
 166 </pre>
 167 <p>
 168         This form integrates just fine over (-log(&#960;/2), +&#8734;) using either the <code class="computeroutput"><span class="identifier">tanh_sinh</span></code> or <code class="computeroutput"><span class="identifier">exp_sinh</span></code>
 169         classes.
 170       </p>
 171 </div>
 172 <table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
 173 <td align="left"></td>
 174 <td align="right"><div class="copyright-footer">Copyright &#169; 2006-2019 Nikhar
 175       Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos,
 176       Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan
 177       R&#229;de, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg,
 178       Daryle Walker and Xiaogang Zhang<p>
 179         Distributed under the Boost Software License, Version 1.0. (See accompanying
 180         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>)
 181       </p>
 182 </div></td>
 183 </tr></table>
 184 <hr>
 185 <div class="spirit-nav">
 186 <a accesskey="p" href="de_thread.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../double_exponential.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="de_refes.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
 187 </div>
 188 </body>
 189 </html>