1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 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.12.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"><</span><span class="keyword">class</span> <span class="identifier">Real</span><span class="special">></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">&</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"><</span><span class="identifier">Real</span><span class="special">>()</span> <span class="special">*</span> <span class="number">4</span><span class="special">)</span> 34 35 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">></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"><</span><span class="identifier">Real</span><span class="special">>(),</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">)-></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"><</span><span class="identifier">F</span><span class="special">>()(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special"><</span><span class="identifier">Real</span><span class="special">>()))</span> <span class="keyword">const</span><span class="special">;</span> 41 42 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">></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"><</span><span class="identifier">Real</span><span class="special">>(),</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">)-></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"><</span><span class="identifier">F</span><span class="special">>()(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">declval</span><span class="special"><</span><span class="identifier">Real</span><span class="special">>()))</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| < 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"><</span><span class="keyword">double</span><span class="special">></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 ± ⅈ/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 ± 5ⅈ), 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"><</span><span class="keyword">double</span><span class="special">></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"><</span><span class="keyword">double</span><span class="special">>::</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">&</span><span class="identifier">error</span><span class="special">,</span> <span class="special">&</span><span class="identifier">L1</span><span class="special">,</span> <span class="special">&</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,+∞) 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 truly 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 (-∞, ∞) 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"><</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"><</span><span class="keyword">class</span> <span class="identifier">Complex</span><span class="special">></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">[&</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"><</span><span class="identifier">value_type</span><span class="special">></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"><</span><span class="identifier">value_type</span><span class="special">>())</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"><</span><span class="identifier">value_type</span><span class="special">>();</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 © 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å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> 543
