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

   1 <html>
   2 <head>
   3 <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
   4 <title>tanh_sinh</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_overview.html" title="Overview">
  10 <link rel="next" href="de_tanh_sinh_2_arg.html" title="Handling functions with large features near an endpoint with tanh-sinh quadrature">
  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_overview.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_tanh_sinh_2_arg.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_tanh_sinh"></a><a class="link" href="de_tanh_sinh.html" title="tanh_sinh">tanh_sinh</a>
  28 </h3></div></div></div>
  29 <pre class="programlisting"><span class="keyword">template</span><span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">Real</span><span class="special">&gt;</span>
  30 <span class="keyword">class</span> <span class="identifier">tanh_sinh</span>
  31 <span class="special">{</span>
  32 <span class="keyword">public</span><span class="special">:</span>
  33     <span class="identifier">tanh_sinh</span><span class="special">(</span><span class="identifier">size_t</span> <span class="identifier">max_refinements</span> <span class="special">=</span> <span class="number">15</span><span class="special">,</span> <span class="keyword">const</span> <span class="identifier">Real</span><span class="special">&amp;</span> <span class="identifier">min_complement</span> <span class="special">=</span> <span class="identifier">tools</span><span class="special">::</span><span class="identifier">min_value</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;()</span> <span class="special">*</span> <span class="number">4</span><span class="special">)</span>
  34
  35     <span class="keyword">template</span><span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">&gt;</span>
  36     <span class="keyword">auto</span> <span class="identifier">integrate</span><span class="special">(</span><span class="keyword">const</span> <span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">a</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">b</span><span class="special">,</span>
  37                    <span class="identifier">Real</span> <span class="identifier">tolerance</span> <span class="special">=</span> <span class="identifier">tools</span><span class="special">::</span><span class="identifier">root_epsilon</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;(),</span>
  38                    <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">error</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">,</span>
  39                    <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">L1</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">,</span>
  40                    <span class="identifier">std</span><span class="special">::</span><span class="identifier">size_t</span><span class="special">*</span> <span class="identifier">levels</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">)-&gt;</span><span class="keyword">decltype</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special">&lt;</span><span class="identifier">F</span><span class="special">&gt;()(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;()))</span> <span class="keyword">const</span><span class="special">;</span>
  41
  42     <span class="keyword">template</span><span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">&gt;</span>
  43     <span class="keyword">auto</span> <span class="identifier">integrate</span><span class="special">(</span><span class="keyword">const</span> <span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span>
  44                    <span class="identifier">tolerance</span> <span class="special">=</span> <span class="identifier">tools</span><span class="special">::</span><span class="identifier">root_epsilon</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;(),</span>
  45                    <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">error</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">,</span>
  46                    <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">L1</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">,</span>
  47                    <span class="identifier">std</span><span class="special">::</span><span class="identifier">size_t</span><span class="special">*</span> <span class="identifier">levels</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">)-&gt;</span><span class="keyword">decltype</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special">&lt;</span><span class="identifier">F</span><span class="special">&gt;()(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;()))</span> <span class="keyword">const</span><span class="special">;</span>
  48
  49 <span class="special">};</span>
  50 </pre>
  51 <p>
  52         The <code class="computeroutput"><span class="identifier">tanh</span><span class="special">-</span><span class="identifier">sinh</span></code> quadrature routine provided by boost
  53         is a rapidly convergent numerical integration scheme for holomorphic integrands.
  54         By this we mean that the integrand is the restriction to the real line of
  55         a complex-differentiable function which is bounded on the interior of the
  56         unit disk <span class="emphasis"><em>|z| &lt; 1</em></span>, so that it lies within the so-called
  57         <a href="https://en.wikipedia.org/wiki/Hardy_space" target="_top">Hardy space</a>.
  58         If your integrand obeys these conditions, it can be shown that <code class="computeroutput"><span class="identifier">tanh</span><span class="special">-</span><span class="identifier">sinh</span></code>
  59         integration is optimal, in the sense that it requires the fewest function
  60         evaluations for a given accuracy of any quadrature algorithm for a random
  61         element from the Hardy space.
  62       </p>
  63 <p>
  64         A basic example of how to use the <code class="computeroutput"><span class="identifier">tanh</span><span class="special">-</span><span class="identifier">sinh</span></code> quadrature
  65         is shown below:
  66       </p>
  67 <pre class="programlisting"><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>
  68 <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> <span class="keyword">return</span> <span class="number">5</span><span class="special">*</span><span class="identifier">x</span> <span class="special">+</span> <span class="number">7</span><span class="special">;</span> <span class="special">};</span>
  69 <span class="comment">// Integrate over native bounds of (-1,1):</span>
  70 <span class="keyword">double</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>
  71 <span class="comment">// Integrate over (0,1.1) instead:</span>
  72 <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.0</span><span class="special">,</span> <span class="number">1.1</span><span class="special">);</span>
  73 </pre>
  74 <p>
  75         The basic idea of <code class="computeroutput"><span class="identifier">tanh</span><span class="special">-</span><span class="identifier">sinh</span></code> quadrature is that a variable transformation
  76         can cause the endpoint derivatives to decay rapidly. When the derivatives
  77         at the endpoints decay much faster than the Bernoulli numbers grow, the Euler-Maclaurin
  78         summation formula tells us that simple trapezoidal quadrature converges faster
  79         than any power of <span class="emphasis"><em>h</em></span>. That means the number of correct
  80         digits of the result should roughly double with each new level of integration
  81         (halving of <span class="emphasis"><em>h</em></span>), Hence the default termination condition
  82         for integration is usually set to the square root of machine epsilon. Most
  83         well-behaved integrals should converge to full machine precision with this
  84         termination condition, and in 6 or fewer levels at double precision, or 7
  85         or fewer levels for quad precision.
  86       </p>
  87 <p>
  88         One very nice property of tanh-sinh quadrature is that it can handle singularities
  89         at the endpoints of the integration domain. For instance, the following integrand,
  90         singular at both endpoints, can be efficiently evaluated to 100 binary digits:
  91       </p>
  92 <pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[](</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="identifier">log</span><span class="special">(</span><span class="identifier">x</span><span class="special">)*</span><span class="identifier">log1p</span><span class="special">(-</span><span class="identifier">x</span><span class="special">);</span> <span class="special">};</span>
  93 <span class="identifier">Real</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="special">(</span><span class="identifier">Real</span><span class="special">)</span> <span class="number">0</span><span class="special">,</span> <span class="special">(</span><span class="identifier">Real</span><span class="special">)</span> <span class="number">1</span><span class="special">);</span>
  94 </pre>
  95 <p>
  96         Now onto the caveats: As stated before, the integrands must lie in a Hardy
  97         space to ensure rapid convergence. Attempting to integrate a function which
  98         is not bounded on the unit disk by tanh-sinh can lead to very slow convergence.
  99         For example, take the Runge function:
 100       </p>
 101 <pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">f1</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">t</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="number">1</span><span class="special">/(</span><span class="number">1</span><span class="special">+</span><span class="number">25</span><span class="special">*</span><span class="identifier">t</span><span class="special">*</span><span class="identifier">t</span><span class="special">);</span> <span class="special">};</span>
 102 <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">f1</span><span class="special">,</span> <span class="special">(</span><span class="keyword">double</span><span class="special">)</span> <span class="special">-</span><span class="number">1</span><span class="special">,</span> <span class="special">(</span><span class="keyword">double</span><span class="special">)</span> <span class="number">1</span><span class="special">);</span>
 103 </pre>
 104 <p>
 105         This function has poles at &#177; &#8520;/5, and as such it is not bounded
 106         on the unit disk. However, the related function
 107       </p>
 108 <pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">f2</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">t</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="number">1</span><span class="special">/(</span><span class="number">1</span><span class="special">+</span><span class="number">0.04</span><span class="special">*</span><span class="identifier">t</span><span class="special">*</span><span class="identifier">t</span><span class="special">);</span> <span class="special">};</span>
 109 <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">f2</span><span class="special">,</span> <span class="special">(</span><span class="keyword">double</span><span class="special">)</span> <span class="special">-</span><span class="number">1</span><span class="special">,</span> <span class="special">(</span><span class="keyword">double</span><span class="special">)</span> <span class="number">1</span><span class="special">);</span>
 110 </pre>
 111 <p>
 112         has poles outside the unit disk (at &#177; 5&#8520;), and is therefore in
 113         the Hardy space. Our benchmarks show that the second integration is performed
 114         22x faster than the first! If you do not understand the structure of your
 115         integrand in the complex plane, do performance testing before deployment.
 116       </p>
 117 <p>
 118         Like the trapezoidal quadrature, the tanh-sinh quadrature produces an estimate
 119         of the L<sub>1</sub> norm of the integral along with the requested integral. This is
 120         to establish a scale against which to measure the tolerance, and to provide
 121         an estimate of the condition number of the summation. This can be queried
 122         as follows:
 123       </p>
 124 <pre class="programlisting"><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>
 125 <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> <span class="keyword">return</span> <span class="number">5</span><span class="special">*</span><span class="identifier">x</span> <span class="special">+</span> <span class="number">7</span><span class="special">;</span> <span class="special">};</span>
 126 <span class="keyword">double</span> <span class="identifier">termination</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">sqrt</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>
 127 <span class="keyword">double</span> <span class="identifier">error</span><span class="special">;</span>
 128 <span class="keyword">double</span> <span class="identifier">L1</span><span class="special">;</span>
 129 <span class="identifier">size_t</span> <span class="identifier">levels</span><span class="special">;</span>
 130 <span class="keyword">double</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.0</span><span class="special">,</span> <span class="number">1.0</span><span class="special">,</span> <span class="identifier">termination</span><span class="special">,</span> <span class="special">&amp;</span><span class="identifier">error</span><span class="special">,</span> <span class="special">&amp;</span><span class="identifier">L1</span><span class="special">,</span> <span class="special">&amp;</span><span class="identifier">levels</span><span class="special">);</span>
 131 <span class="keyword">double</span> <span class="identifier">condition_number</span> <span class="special">=</span> <span class="identifier">L1</span><span class="special">/</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">abs</span><span class="special">(</span><span class="identifier">Q</span><span class="special">);</span>
 132 </pre>
 133 <p>
 134         If the condition number is large, the computed integral is worthless: typically
 135         one can assume that Q has lost one digit of precision when the condition
 136         number of O(10^Q). The returned error term is not the actual error in the
 137         result, but merely an a posteriori error estimate. It is the absolute difference
 138         between the last two approximations, and for well behaved integrals, the
 139         actual error should be very much smaller than this. The following table illustrates
 140         how the errors and conditioning vary for few sample integrals, in each case
 141         the termination condition was set to the square root of epsilon, and all
 142         tests were conducted in double precision:
 143       </p>
 144 <div class="informaltable"><table class="table">
 145 <colgroup>
 146 <col>
 147 <col>
 148 <col>
 149 <col>
 150 <col>
 151 <col>
 152 <col>
 153 </colgroup>
 154 <thead><tr>
 155 <th>
 156                 <p>
 157                   Integral
 158                 </p>
 159               </th>
 160 <th>
 161                 <p>
 162                   Range
 163                 </p>
 164               </th>
 165 <th>
 166                 <p>
 167                   Error
 168                 </p>
 169               </th>
 170 <th>
 171                 <p>
 172                   Actual measured error
 173                 </p>
 174               </th>
 175 <th>
 176                 <p>
 177                   Levels
 178                 </p>
 179               </th>
 180 <th>
 181                 <p>
 182                   Condition Number
 183                 </p>
 184               </th>
 185 <th>
 186                 <p>
 187                   Comments
 188                 </p>
 189               </th>
 190 </tr></thead>
 191 <tbody>
 192 <tr>
 193 <td>
 194                 <p>
 195                   <code class="computeroutput"><span class="number">5</span> <span class="special">*</span>
 196                   <span class="identifier">x</span> <span class="special">+</span>
 197                   <span class="number">7</span></code>
 198                 </p>
 199               </td>
 200 <td>
 201                 <p>
 202                   (0,1)
 203                 </p>
 204               </td>
 205 <td>
 206                 <p>
 207                   3.5e-15
 208                 </p>
 209               </td>
 210 <td>
 211                 <p>
 212                   0
 213                 </p>
 214               </td>
 215 <td>
 216                 <p>
 217                   5
 218                 </p>
 219               </td>
 220 <td>
 221                 <p>
 222                   1
 223                 </p>
 224               </td>
 225 <td>
 226                 <p>
 227                   This trivial case shows just how accurate these methods can be.
 228                 </p>
 229               </td>
 230 </tr>
 231 <tr>
 232 <td>
 233                 <p>
 234                   <code class="computeroutput"><span class="identifier">log</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span>
 235                   <span class="special">*</span> <span class="identifier">log</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span></code>
 236                 </p>
 237               </td>
 238 <td>
 239                 <p>
 240                   0, 1)
 241                 </p>
 242               </td>
 243 <td>
 244                 <p>
 245                   0
 246                 </p>
 247               </td>
 248 <td>
 249                 <p>
 250                   0
 251                 </p>
 252               </td>
 253 <td>
 254                 <p>
 255                   5
 256                 </p>
 257               </td>
 258 <td>
 259                 <p>
 260                   1
 261                 </p>
 262               </td>
 263 <td>
 264                 <p>
 265                   This is an example of an integral that Gaussian integrators fail
 266                   to handle.
 267                 </p>
 268               </td>
 269 </tr>
 270 <tr>
 271 <td>
 272                 <p>
 273                   <code class="computeroutput"><span class="identifier">exp</span><span class="special">(-</span><span class="identifier">x</span><span class="special">)</span>
 274                   <span class="special">/</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span></code>
 275                 </p>
 276               </td>
 277 <td>
 278                 <p>
 279                   (0,+&#8734;)
 280                 </p>
 281               </td>
 282 <td>
 283                 <p>
 284                   8.0e-10
 285                 </p>
 286               </td>
 287 <td>
 288                 <p>
 289                   1.1e-15
 290                 </p>
 291               </td>
 292 <td>
 293                 <p>
 294                   5
 295                 </p>
 296               </td>
 297 <td>
 298                 <p>
 299                   1
 300                 </p>
 301               </td>
 302 <td>
 303                 <p>
 304                   Gaussian integrators typically fail to handle the singularities
 305                   at the endpoints of this one.
 306                 </p>
 307               </td>
 308 </tr>
 309 <tr>
 310 <td>
 311                 <p>
 312                   <code class="computeroutput"><span class="identifier">x</span> <span class="special">*</span>
 313                   <span class="identifier">sin</span><span class="special">(</span><span class="number">2</span> <span class="special">*</span> <span class="identifier">exp</span><span class="special">(</span><span class="number">2</span> <span class="special">*</span> <span class="identifier">sin</span><span class="special">(</span><span class="number">2</span> <span class="special">*</span> <span class="identifier">exp</span><span class="special">(</span><span class="number">2</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">))))</span></code>
 314                 </p>
 315               </td>
 316 <td>
 317                 <p>
 318                   (-1,1)
 319                 </p>
 320               </td>
 321 <td>
 322                 <p>
 323                   7.2e-16
 324                 </p>
 325               </td>
 326 <td>
 327                 <p>
 328                   4.9e-17
 329                 </p>
 330               </td>
 331 <td>
 332                 <p>
 333                   9
 334                 </p>
 335               </td>
 336 <td>
 337                 <p>
 338                   1.89
 339                 </p>
 340               </td>
 341 <td>
 342                 <p>
 343                   This is a truely horrible integral that oscillates wildly and unpredictably
 344                   with some very sharp "spikes" in it's graph. The higher
 345                   number of levels used reflects the difficulty of sampling the more
 346                   extreme features.
 347                 </p>
 348               </td>
 349 </tr>
 350 <tr>
 351 <td>
 352                 <p>
 353                   <code class="computeroutput"><span class="identifier">x</span> <span class="special">==</span>
 354                   <span class="number">0</span> <span class="special">?</span>
 355                   <span class="number">1</span> <span class="special">:</span>
 356                   <span class="identifier">sin</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span>
 357                   <span class="special">/</span> <span class="identifier">x</span></code>
 358                 </p>
 359               </td>
 360 <td>
 361                 <p>
 362                   (-&#8734;, &#8734;)
 363                 </p>
 364               </td>
 365 <td>
 366                 <p>
 367                   3.0e-1
 368                 </p>
 369               </td>
 370 <td>
 371                 <p>
 372                   4.0e-1
 373                 </p>
 374               </td>
 375 <td>
 376                 <p>
 377                   15
 378                 </p>
 379               </td>
 380 <td>
 381                 <p>
 382                   159
 383                 </p>
 384               </td>
 385 <td>
 386                 <p>
 387                   This highly oscillatory integral isn't handled at all well by tanh-sinh
 388                   quadrature: there is so much cancellation in the sum that the result
 389                   is essentially worthless. The argument transformation of the infinite
 390                   integral behaves somewhat badly as well, in fact we do <span class="emphasis"><em>slightly</em></span>
 391                   better integrating over 2 symmetrical and large finite limits.
 392                 </p>
 393               </td>
 394 </tr>
 395 <tr>
 396 <td>
 397                 <p>
 398                   <code class="computeroutput"><span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">x</span> <span class="special">/</span>
 399                   <span class="special">(</span><span class="number">1</span>
 400                   <span class="special">-</span> <span class="identifier">x</span>
 401                   <span class="special">*</span> <span class="identifier">x</span><span class="special">))</span></code>
 402                 </p>
 403               </td>
 404 <td>
 405                 <p>
 406                   (0,1)
 407                 </p>
 408               </td>
 409 <td>
 410                 <p>
 411                   1e-8
 412                 </p>
 413               </td>
 414 <td>
 415                 <p>
 416                   1e-8
 417                 </p>
 418               </td>
 419 <td>
 420                 <p>
 421                   5
 422                 </p>
 423               </td>
 424 <td>
 425                 <p>
 426                   1
 427                 </p>
 428               </td>
 429 <td>
 430                 <p>
 431                   This an example of an integral that has all its area close to a
 432                   non-zero endpoint, the problem here is that the function being
 433                   integrated returns "garbage" values for x very close
 434                   to 1. We can easily fix this issue by passing a 2 argument functor
 435                   to the integrator: the second argument gives the distance to the
 436                   nearest endpoint, and we can use that information to return accurate
 437                   values, and thus fix the integral calculation.
 438                 </p>
 439               </td>
 440 </tr>
 441 <tr>
 442 <td>
 443                 <p>
 444                   <code class="computeroutput"><span class="identifier">x</span> <span class="special">&lt;</span>
 445                   <span class="number">0.5</span> <span class="special">?</span>
 446                   <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span>
 447                   <span class="special">/</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="number">1</span> <span class="special">-</span> <span class="identifier">x</span>
 448                   <span class="special">*</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">:</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">x</span> <span class="special">/</span>
 449                   <span class="special">((</span><span class="identifier">x</span>
 450                   <span class="special">+</span> <span class="number">1</span><span class="special">)</span> <span class="special">*</span> <span class="special">(</span><span class="identifier">xc</span><span class="special">)))</span></code>
 451                 </p>
 452               </td>
 453 <td>
 454                 <p>
 455                   (0,1)
 456                 </p>
 457               </td>
 458 <td>
 459                 <p>
 460                   0
 461                 </p>
 462               </td>
 463 <td>
 464                 <p>
 465                   0
 466                 </p>
 467               </td>
 468 <td>
 469                 <p>
 470                   5
 471                 </p>
 472               </td>
 473 <td>
 474                 <p>
 475                   1
 476                 </p>
 477               </td>
 478 <td>
 479                 <p>
 480                   This is the 2-argument version of the previous integral, the second
 481                   argument <span class="emphasis"><em>xc</em></span> is <code class="computeroutput"><span class="number">1</span><span class="special">-</span><span class="identifier">x</span></code>
 482                   in this case, and we use 1-x<sup>2</sup> == (1-x)(1+x) to calculate 1-x<sup>2</sup> with
 483                   greater accuracy.
 484                 </p>
 485               </td>
 486 </tr>
 487 </tbody>
 488 </table></div>
 489 <p>
 490         Although the <code class="computeroutput"><span class="identifier">tanh</span><span class="special">-</span><span class="identifier">sinh</span></code> quadrature can compute integral over
 491         infinite domains by variable transformations, these transformations can create
 492         a very poorly behaved integrand. For this reason, double-exponential variable
 493         transformations have been provided that allow stable computation over infinite
 494         domains; these being the exp-sinh and sinh-sinh quadrature.
 495       </p>
 496 <h5>
 497 <a name="math_toolkit.double_exponential.de_tanh_sinh.h0"></a>
 498         <span class="phrase"><a name="math_toolkit.double_exponential.de_tanh_sinh.complex_integrals"></a></span><a class="link" href="de_tanh_sinh.html#math_toolkit.double_exponential.de_tanh_sinh.complex_integrals">Complex
 499         integrals</a>
 500       </h5>
 501 <p>
 502         The <code class="computeroutput"><span class="identifier">tanh_sinh</span></code> integrator
 503         supports integration of functions which return complex results, for example
 504         the sine-integral <code class="computeroutput"><span class="identifier">Si</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span></code> has
 505         the integral representation:
 506       </p>
 507 <div class="blockquote"><blockquote class="blockquote"><p>
 508           <span class="inlinemediaobject"><img src="../../../equations/sine_integral.svg"></span>
 509
 510         </p></blockquote></div>
 511 <p>
 512         Which we can code up directly as:
 513       </p>
 514 <pre class="programlisting"><span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">Complex</span><span class="special">&gt;</span>
 515 <span class="identifier">Complex</span> <span class="identifier">Si</span><span class="special">(</span><span class="identifier">Complex</span> <span class="identifier">z</span><span class="special">)</span>
 516 <span class="special">{</span>
 517    <span class="keyword">typedef</span> <span class="keyword">typename</span> <span class="identifier">Complex</span><span class="special">::</span><span class="identifier">value_type</span> <span class="identifier">value_type</span><span class="special">;</span>
 518    <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">sin</span><span class="special">;</span>  <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cos</span><span class="special">;</span> <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">exp</span><span class="special">;</span>
 519    <span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[&amp;</span><span class="identifier">z</span><span class="special">](</span><span class="identifier">value_type</span> <span class="identifier">t</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="special">-</span><span class="identifier">exp</span><span class="special">(-</span><span class="identifier">z</span> <span class="special">*</span> <span class="identifier">cos</span><span class="special">(</span><span class="identifier">t</span><span class="special">))</span> <span class="special">*</span> <span class="identifier">cos</span><span class="special">(</span><span class="identifier">z</span> <span class="special">*</span> <span class="identifier">sin</span><span class="special">(</span><span class="identifier">t</span><span class="special">));</span> <span class="special">};</span>
 520    <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">quadrature</span><span class="special">::</span><span class="identifier">tanh_sinh</span><span class="special">&lt;</span><span class="identifier">value_type</span><span class="special">&gt;</span> <span class="identifier">integrator</span><span class="special">;</span>
 521    <span class="keyword">return</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">boost</span><span class="special">::</span><span class="identifier">math</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="identifier">value_type</span><span class="special">&gt;())</span> <span class="special">+</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</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="identifier">value_type</span><span class="special">&gt;();</span>
 522 <span class="special">}</span>
 523 </pre>
 524 </div>
 525 <table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
 526 <td align="left"></td>
 527 <td align="right"><div class="copyright-footer">Copyright &#169; 2006-2019 Nikhar
 528       Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos,
 529       Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan
 530       R&#229;de, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg,
 531       Daryle Walker and Xiaogang Zhang<p>
 532         Distributed under the Boost Software License, Version 1.0. (See accompanying
 533         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>)
 534       </p>
 535 </div></td>
 536 </tr></table>
 537 <hr>
 538 <div class="spirit-nav">
 539 <a accesskey="p" href="de_overview.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_tanh_sinh_2_arg.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
 540 </div>
 541 </body>
 542 </html>