boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-23 05:42:15 +00:00
Code Issues Projects Releases Wiki Activity
Files
develop
math/doc/html/math_toolkit/double_exponential/de_tanh_sinh.html
jzmaddock 4c2cdc2cf6 Regenerate docs.
2024-08-06 13:18:32 +01:00

541 lines
37 KiB
HTML
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<title>tanh_sinh</title>
<link rel="stylesheet" href="../../math.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
<link rel="home" href="../../index.html" title="Math Toolkit 4.2.1">
<link rel="up" href="../double_exponential.html" title="Double-exponential quadrature">
<link rel="prev" href="de_overview.html" title="Overview">
<link rel="next" href="de_tanh_sinh_2_arg.html" title="Handling functions with large features near an endpoint with tanh-sinh quadrature">
<meta name="viewport" content="width=device-width, initial-scale=1">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../boost.png"></td>
<td align="center"><a href="../../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<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>
</div>
<div class="section">
<div class="titlepage"><div><div><h3 class="title">
<a name="math_toolkit.double_exponential.de_tanh_sinh"></a><a class="link" href="de_tanh_sinh.html" title="tanh_sinh">tanh_sinh</a>
</h3></div></div></div>
<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>
<span class="keyword">class</span> <span class="identifier">tanh_sinh</span>
<span class="special">{</span>
<span class="keyword">public</span><span class="special">:</span>
<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>
<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>
<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>
<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>
<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>
<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>
<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>
<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>
<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">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>
<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>
<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>
<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>
<span class="special">};</span>
</pre>
<p>
The <code class="computeroutput"><span class="identifier">tanh</span><span class="special">-</span><span class="identifier">sinh</span></code> quadrature routine provided by boost
is a rapidly convergent numerical integration scheme for holomorphic integrands.
By this we mean that the integrand is the restriction to the real line of
a complex-differentiable function which is bounded on the interior of the
unit disk <span class="emphasis"><em>|z| &lt; 1</em></span>, so that it lies within the so-called
<a href="https://en.wikipedia.org/wiki/Hardy_space" target="_top">Hardy space</a>.
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>
integration is optimal, in the sense that it requires the fewest function
evaluations for a given accuracy of any quadrature algorithm for a random
element from the Hardy space.
</p>
<p>
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
is shown below:
</p>
<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>
<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>
<span class="comment">// Integrate over native bounds of (-1,1):</span>
<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="comment">// Integrate over (0,1.1) instead:</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.1</span><span class="special">);</span>
</pre>
<p>
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
can cause the endpoint derivatives to decay rapidly. When the derivatives
at the endpoints decay much faster than the Bernoulli numbers grow, the Euler-Maclaurin
summation formula tells us that simple trapezoidal quadrature converges faster
than any power of <span class="emphasis"><em>h</em></span>. That means the number of correct
digits of the result should roughly double with each new level of integration
(halving of <span class="emphasis"><em>h</em></span>), Hence the default termination condition
for integration is usually set to the square root of machine epsilon. Most
well-behaved integrals should converge to full machine precision with this
termination condition, and in 6 or fewer levels at double precision, or 7
or fewer levels for quad precision.
</p>
<p>
One very nice property of tanh-sinh quadrature is that it can handle singularities
at the endpoints of the integration domain. For instance, the following integrand,
singular at both endpoints, can be efficiently evaluated to 100 binary digits:
</p>
<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>
<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>
</pre>
<p>
Now onto the caveats: As stated before, the integrands must lie in a Hardy
space to ensure rapid convergence. Attempting to integrate a function which
is not bounded on the unit disk by tanh-sinh can lead to very slow convergence.
For example, take the Runge function:
</p>
<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>
<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>
</pre>
<p>
This function has poles at ± ⅈ/5, and as such it is not bounded
on the unit disk. However, the related function
</p>
<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>
<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>
</pre>
<p>
has poles outside the unit disk (at ± 5ⅈ), and is therefore in
the Hardy space. Our benchmarks show that the second integration is performed
22x faster than the first! If you do not understand the structure of your
integrand in the complex plane, do performance testing before deployment.
</p>
<p>
Like the trapezoidal quadrature, the tanh-sinh quadrature produces an estimate
of the L<sub>1</sub> norm of the integral along with the requested integral. This is
to establish a scale against which to measure the tolerance, and to provide
an estimate of the condition number of the summation. This can be queried
as follows:
</p>
<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>
<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>
<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>
<span class="keyword">double</span> <span class="identifier">error</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">L1</span><span class="special">;</span>
<span class="identifier">size_t</span> <span class="identifier">levels</span><span class="special">;</span>
<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>
<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>
</pre>
<p>
If the condition number is large, the computed integral is worthless: typically
one can assume that Q has lost one digit of precision when the condition
number of O(10^Q). The returned error term is not the actual error in the
result, but merely an a posteriori error estimate. It is the absolute difference
between the last two approximations, and for well behaved integrals, the
actual error should be very much smaller than this. The following table illustrates
how the errors and conditioning vary for few sample integrals, in each case
the termination condition was set to the square root of epsilon, and all
tests were conducted in double precision:
</p>
<div class="informaltable"><table class="table">
<colgroup>
<col>
<col>
<col>
<col>
<col>
<col>
<col>
</colgroup>
<thead><tr>
<th>
<p>
Integral
</p>
</th>
<th>
<p>
Range
</p>
</th>
<th>
<p>
Error
</p>
</th>
<th>
<p>
Actual measured error
</p>
</th>
<th>
<p>
Levels
</p>
</th>
<th>
<p>
Condition Number
</p>
</th>
<th>
<p>
Comments
</p>
</th>
</tr></thead>
<tbody>
<tr>
<td>
<p>
<code class="computeroutput"><span class="number">5</span> <span class="special">*</span>
<span class="identifier">x</span> <span class="special">+</span>
<span class="number">7</span></code>
</p>
</td>
<td>
<p>
(0,1)
</p>
</td>
<td>
<p>
3.5e-15
</p>
</td>
<td>
<p>
0
</p>
</td>
<td>
<p>
5
</p>
</td>
<td>
<p>
1
</p>
</td>
<td>
<p>
This trivial case shows just how accurate these methods can be.
</p>
</td>
</tr>
<tr>
<td>
<p>
<code class="computeroutput"><span class="identifier">log</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span>
<span class="special">*</span> <span class="identifier">log</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span></code>
</p>
</td>
<td>
<p>
0, 1)
</p>
</td>
<td>
<p>
0
</p>
</td>
<td>
<p>
0
</p>
</td>
<td>
<p>
5
</p>
</td>
<td>
<p>
1
</p>
</td>
<td>
<p>
This is an example of an integral that Gaussian integrators fail
to handle.
</p>
</td>
</tr>
<tr>
<td>
<p>
<code class="computeroutput"><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">sqrt</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span></code>
</p>
</td>
<td>
<p>
(0,+∞)
</p>
</td>
<td>
<p>
8.0e-10
</p>
</td>
<td>
<p>
1.1e-15
</p>
</td>
<td>
<p>
5
</p>
</td>
<td>
<p>
1
</p>
</td>
<td>
<p>
Gaussian integrators typically fail to handle the singularities
at the endpoints of this one.
</p>
</td>
</tr>
<tr>
<td>
<p>
<code class="computeroutput"><span class="identifier">x</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">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>
</p>
</td>
<td>
<p>
(-1,1)
</p>
</td>
<td>
<p>
7.2e-16
</p>
</td>
<td>
<p>
4.9e-17
</p>
</td>
<td>
<p>
9
</p>
</td>
<td>
<p>
1.89
</p>
</td>
<td>
<p>
This is a truly horrible integral that oscillates wildly and unpredictably
with some very sharp "spikes" in it's graph. The higher
number of levels used reflects the difficulty of sampling the more
extreme features.
</p>
</td>
</tr>
<tr>
<td>
<p>
<code class="computeroutput"><span class="identifier">x</span> <span class="special">==</span>
<span class="number">0</span> <span class="special">?</span>
<span class="number">1</span> <span class="special">:</span>
<span class="identifier">sin</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span>
<span class="special">/</span> <span class="identifier">x</span></code>
</p>
</td>
<td>
<p>
(-∞, ∞)
</p>
</td>
<td>
<p>
3.0e-1
</p>
</td>
<td>
<p>
4.0e-1
</p>
</td>
<td>
<p>
15
</p>
</td>
<td>
<p>
159
</p>
</td>
<td>
<p>
This highly oscillatory integral isn't handled at all well by tanh-sinh
quadrature: there is so much cancellation in the sum that the result
is essentially worthless. The argument transformation of the infinite
integral behaves somewhat badly as well, in fact we do <span class="emphasis"><em>slightly</em></span>
better integrating over 2 symmetrical and large finite limits.
</p>
</td>
</tr>
<tr>
<td>
<p>
<code class="computeroutput"><span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">x</span> <span class="special">/</span>
<span class="special">(</span><span class="number">1</span>
<span class="special">-</span> <span class="identifier">x</span>
<span class="special">*</span> <span class="identifier">x</span><span class="special">))</span></code>
</p>
</td>
<td>
<p>
(0,1)
</p>
</td>
<td>
<p>
1e-8
</p>
</td>
<td>
<p>
1e-8
</p>
</td>
<td>
<p>
5
</p>
</td>
<td>
<p>
1
</p>
</td>
<td>
<p>
This an example of an integral that has all its area close to a
non-zero endpoint, the problem here is that the function being
integrated returns "garbage" values for x very close
to 1. We can easily fix this issue by passing a 2 argument functor
to the integrator: the second argument gives the distance to the
nearest endpoint, and we can use that information to return accurate
values, and thus fix the integral calculation.
</p>
</td>
</tr>
<tr>
<td>
<p>
<code class="computeroutput"><span class="identifier">x</span> <span class="special">&lt;</span>
<span class="number">0.5</span> <span class="special">?</span>
<span class="identifier">sqrt</span><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="number">1</span> <span class="special">-</span> <span class="identifier">x</span>
<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>
<span class="special">((</span><span class="identifier">x</span>
<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>
</p>
</td>
<td>
<p>
(0,1)
</p>
</td>
<td>
<p>
0
</p>
</td>
<td>
<p>
0
</p>
</td>
<td>
<p>
5
</p>
</td>
<td>
<p>
1
</p>
</td>
<td>
<p>
This is the 2-argument version of the previous integral, the second
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>
in this case, and we use 1-x<sup>2</sup> == (1-x)(1+x) to calculate 1-x<sup>2</sup> with
greater accuracy.
</p>
</td>
</tr>
</tbody>
</table></div>
<p>
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
infinite domains by variable transformations, these transformations can create
a very poorly behaved integrand. For this reason, double-exponential variable
transformations have been provided that allow stable computation over infinite
domains; these being the exp-sinh and sinh-sinh quadrature.
</p>
<h5>
<a name="math_toolkit.double_exponential.de_tanh_sinh.h0"></a>
<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
integrals</a>
</h5>
<p>
The <code class="computeroutput"><span class="identifier">tanh_sinh</span></code> integrator
supports integration of functions which return complex results, for example
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
the integral representation:
</p>
<div class="blockquote"><blockquote class="blockquote"><p>
<span class="inlinemediaobject"><img src="../../../equations/sine_integral.svg"></span>
</p></blockquote></div>
<p>
Which we can code up directly as:
</p>
<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>
<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>
<span class="special">{</span>
<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>
<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>
<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>
<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>
<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>
<span class="special">}</span>
</pre>
</div>
<div class="copyright-footer">Copyright © 2006-2021 Nikhar Agrawal, Anton Bikineev, Matthew Borland,
Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno
Lalande, John Maddock, Evan Miller, Jeremy Murphy, Matthew Pulver, Johan Råde,
Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, Daryle
Walker and Xiaogang Zhang<p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
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>)
</p>
</div>
<hr>
<div class="spirit-nav">
<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>
</div>
</body>
</html>
Reference in New Issue View Git Blame Copy Permalink