318 lines
34 KiB
HTML
318 lines
34 KiB
HTML
<html>
|
|
<head>
|
|
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
|
|
<title>Fourier Integrals</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 2.11.0">
|
|
<link rel="up" href="../quadrature.html" title="Chapter 13. Quadrature and Differentiation">
|
|
<link rel="prev" href="double_exponential/de_refes.html" title="References">
|
|
<link rel="next" href="naive_monte_carlo.html" title="Naive Monte Carlo Integration">
|
|
</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="double_exponential/de_refes.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../quadrature.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="naive_monte_carlo.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a>
|
|
</div>
|
|
<div class="section">
|
|
<div class="titlepage"><div><div><h2 class="title" style="clear: both">
|
|
<a name="math_toolkit.fourier_integrals"></a><a class="link" href="fourier_integrals.html" title="Fourier Integrals">Fourier Integrals</a>
|
|
</h2></div></div></div>
|
|
<h4>
|
|
<a name="math_toolkit.fourier_integrals.h0"></a>
|
|
<span class="phrase"><a name="math_toolkit.fourier_integrals.synopsis"></a></span><a class="link" href="fourier_integrals.html#math_toolkit.fourier_integrals.synopsis">Synopsis</a>
|
|
</h4>
|
|
<pre class="programlisting"><span class="preprocessor">#include</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">ooura_fourier_integrals</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span>
|
|
|
|
<span class="keyword">namespace</span> <span class="identifier">boost</span> <span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">math</span> <span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">quadrature</span> <span class="special">{</span>
|
|
|
|
<span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">Real</span><span class="special">></span>
|
|
<span class="keyword">class</span> <span class="identifier">ooura_fourier_sin</span> <span class="special">{</span>
|
|
<span class="keyword">public</span><span class="special">:</span>
|
|
<span class="identifier">ooura_fourier_sin</span><span class="special">(</span><span class="keyword">const</span> <span class="identifier">Real</span> <span class="identifier">relative_error_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> <span class="identifier">size_t</span> <span class="identifier">levels</span> <span class="special">=</span> <span class="keyword">sizeof</span><span class="special">(</span><span class="identifier">Real</span><span class="special">));</span>
|
|
|
|
<span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">></span>
|
|
<span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="identifier">Real</span><span class="special">,</span> <span class="identifier">Real</span><span class="special">></span> <span class="identifier">integrate</span><span class="special">(</span><span class="identifier">F</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">omega</span><span class="special">);</span>
|
|
|
|
<span class="special">};</span>
|
|
|
|
|
|
<span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">Real</span><span class="special">></span>
|
|
<span class="keyword">class</span> <span class="identifier">ooura_fourier_cos</span> <span class="special">{</span>
|
|
<span class="keyword">public</span><span class="special">:</span>
|
|
<span class="identifier">ooura_fourier_cos</span><span class="special">(</span><span class="keyword">const</span> <span class="identifier">Real</span> <span class="identifier">relative_error_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> <span class="identifier">size_t</span> <span class="identifier">levels</span> <span class="special">=</span> <span class="keyword">sizeof</span><span class="special">(</span><span class="identifier">Real</span><span class="special">))</span>
|
|
|
|
<span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">></span>
|
|
<span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="identifier">Real</span><span class="special">,</span> <span class="identifier">Real</span><span class="special">></span> <span class="identifier">integrate</span><span class="special">(</span><span class="identifier">F</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">omega</span><span class="special">);</span>
|
|
<span class="special">};</span>
|
|
|
|
<span class="special">}}}</span> <span class="comment">// namespaces</span>
|
|
</pre>
|
|
<p>
|
|
Ooura's method for Fourier integrals computes
|
|
</p>
|
|
<div class="blockquote"><blockquote class="blockquote"><p>
|
|
<span class="serif_italic">∫<sub>0</sub><sup>∞</sup> f(t)sin(ω t) dt</span>
|
|
</p></blockquote></div>
|
|
<p>
|
|
and
|
|
</p>
|
|
<div class="blockquote"><blockquote class="blockquote"><p>
|
|
<span class="serif_italic">∫<sub>0</sub><sup>∞</sup> f(t)cos(ω t) dt</span>
|
|
</p></blockquote></div>
|
|
<p>
|
|
by a double exponentially decaying transformation. These integrals arise when
|
|
computing continuous Fourier transform of odd and even functions, respectively.
|
|
Oscillatory integrals are known to cause trouble for standard quadrature methods,
|
|
so these routines are provided to cope with the most common oscillatory use
|
|
case.
|
|
</p>
|
|
<p>
|
|
The basic usage is shown below:
|
|
</p>
|
|
<pre class="programlisting"><span class="identifier">ooura_fourier_sin</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span><span class="identifier">integrator</span> <span class="special">=</span> <span class="identifier">ooura_fourier_sin</span><span class="special"><</span><span class="keyword">double</span><span class="special">>();</span>
|
|
<span class="comment">// Use the default tolerance root_epsilon and eight levels for type double.</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="comment">// Simple reciprocal function for sinc.</span>
|
|
<span class="keyword">return</span> <span class="number">1</span> <span class="special">/</span> <span class="identifier">x</span><span class="special">;</span>
|
|
<span class="special">};</span>
|
|
|
|
<span class="keyword">double</span> <span class="identifier">omega</span> <span class="special">=</span> <span class="number">1</span><span class="special">;</span>
|
|
<span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="identifier">result</span> <span class="special">=</span> <span class="identifier">integrator</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="identifier">omega</span><span class="special">);</span>
|
|
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Integral = "</span> <span class="special"><<</span> <span class="identifier">result</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">", relative error estimate "</span> <span class="special"><<</span> <span class="identifier">result</span><span class="special">.</span><span class="identifier">second</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
|
|
</pre>
|
|
<p>
|
|
and compare with the expected value π/2 of the integral.
|
|
</p>
|
|
<pre class="programlisting"><span class="keyword">constexpr</span> <span class="keyword">double</span> <span class="identifier">expected</span> <span class="special">=</span> <span class="identifier">half_pi</span><span class="special"><</span><span class="keyword">double</span><span class="special">>();</span>
|
|
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"pi/2 = "</span> <span class="special"><<</span> <span class="identifier">expected</span> <span class="special"><<</span> <span class="string">", difference "</span> <span class="special"><<</span> <span class="identifier">result</span><span class="special">.</span><span class="identifier">first</span> <span class="special">-</span> <span class="identifier">expected</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
|
|
</pre>
|
|
<p>
|
|
The output is
|
|
</p>
|
|
<pre class="programlisting"><span class="identifier">integral</span> <span class="special">=</span> <span class="number">1.5707963267948966</span><span class="special">,</span> <span class="identifier">relative</span> <span class="identifier">error</span> <span class="identifier">estimate</span> <span class="number">1.2655356398390254e-11</span>
|
|
<span class="identifier">pi</span><span class="special">/</span><span class="number">2</span> <span class="special">=</span> <span class="number">1.5707963267948966</span><span class="special">,</span> <span class="identifier">difference</span> <span class="number">0</span>
|
|
</pre>
|
|
<div class="note"><table border="0" summary="Note">
|
|
<tr>
|
|
<td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../../doc/src/images/note.png"></td>
|
|
<th align="left">Note</th>
|
|
</tr>
|
|
<tr><td align="left" valign="top"><p>
|
|
This integrator is more insistent about examining the error estimate, than
|
|
(say) tanh-sinh, which just returns the value of the integral.
|
|
</p></td></tr>
|
|
</table></div>
|
|
<p>
|
|
With the macro BOOST_MATH_INSTRUMENT_OOURA defined, we can follow the progress:
|
|
</p>
|
|
<pre class="programlisting"><span class="identifier">ooura_fourier_sin</span> <span class="identifier">with</span> <span class="identifier">relative</span> <span class="identifier">error</span> <span class="identifier">goal</span> <span class="number">1.4901161193847656e-08</span> <span class="special">&</span> <span class="number">8</span> <span class="identifier">levels</span><span class="special">.</span>
|
|
<span class="identifier">h</span> <span class="special">=</span> <span class="number">1.000000000000000</span><span class="special">,</span> <span class="identifier">I_h</span> <span class="special">=</span> <span class="number">1.571890732004545</span> <span class="special">=</span> <span class="number">0x1</span><span class="special">.</span><span class="number">92676e56d</span><span class="number">853500</span><span class="identifier">p</span><span class="special">+</span><span class="number">0</span><span class="special">,</span> <span class="identifier">absolute</span> <span class="identifier">error</span> <span class="identifier">estimate</span> <span class="special">=</span> <span class="identifier">nan</span>
|
|
<span class="identifier">h</span> <span class="special">=</span> <span class="number">0.500000000000000</span><span class="special">,</span> <span class="identifier">I_h</span> <span class="special">=</span> <span class="number">1.570793292491940</span> <span class="special">=</span> <span class="number">0x1</span><span class="special">.</span><span class="number">921f</span><span class="number">825</span><span class="identifier">c076f600p</span><span class="special">+</span><span class="number">0</span><span class="special">,</span> <span class="identifier">absolute</span> <span class="identifier">error</span> <span class="identifier">estimate</span> <span class="special">=</span> <span class="number">1.097439512605325e-03</span>
|
|
<span class="identifier">h</span> <span class="special">=</span> <span class="number">0.250000000000000</span><span class="special">,</span> <span class="identifier">I_h</span> <span class="special">=</span> <span class="number">1.570796326814776</span> <span class="special">=</span> <span class="number">0x1</span><span class="special">.</span><span class="number">921f</span><span class="identifier">b54458acf00p</span><span class="special">+</span><span class="number">0</span><span class="special">,</span> <span class="identifier">absolute</span> <span class="identifier">error</span> <span class="identifier">estimate</span> <span class="special">=</span> <span class="number">3.034322835882008e-06</span>
|
|
<span class="identifier">h</span> <span class="special">=</span> <span class="number">0.125000000000000</span><span class="special">,</span> <span class="identifier">I_h</span> <span class="special">=</span> <span class="number">1.570796326794897</span> <span class="special">=</span> <span class="number">0x1</span><span class="special">.</span><span class="number">921f</span><span class="identifier">b54442d1800p</span><span class="special">+</span><span class="number">0</span><span class="special">,</span> <span class="identifier">absolute</span> <span class="identifier">error</span> <span class="identifier">estimate</span> <span class="special">=</span> <span class="number">1.987898734512328e-11</span>
|
|
<span class="identifier">Integral</span> <span class="special">=</span> <span class="number">1.570796326794897e+00</span><span class="special">,</span> <span class="identifier">relative</span> <span class="identifier">error</span> <span class="identifier">estimate</span> <span class="number">1.265535639839025e-11</span>
|
|
<span class="identifier">pi</span><span class="special">/</span><span class="number">2</span> <span class="special">=</span> <span class="number">1.570796326794897e+00</span><span class="special">,</span> <span class="identifier">difference</span> <span class="number">0.000000000000000e+00</span>
|
|
</pre>
|
|
<p>
|
|
Working code of this example is at <a href="../../../example/ooura_fourier_integrals_example.cpp" target="_top">ooura_fourier_integrals_example.cpp</a>
|
|
</p>
|
|
<p>
|
|
A classical cosine transform is presented below:
|
|
</p>
|
|
<pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">integrator</span> <span class="special">=</span> <span class="identifier">ooura_fourier_cos</span><span class="special"><</span><span class="keyword">double</span><span class="special">>();</span>
|
|
<span class="comment">// Use the default tolerance root_epsilon and eight levels for type double.</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="comment">// More complex example function.</span>
|
|
<span class="keyword">return</span> <span class="number">1</span> <span class="special">/</span> <span class="special">(</span><span class="identifier">x</span> <span class="special">*</span> <span class="identifier">x</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="identifier">omega</span> <span class="special">=</span> <span class="number">1</span><span class="special">;</span>
|
|
|
|
<span class="keyword">auto</span> <span class="special">[</span><span class="identifier">result</span><span class="special">,</span> <span class="identifier">relative_error</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">integrator</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="identifier">omega</span><span class="special">);</span>
|
|
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Integral = "</span> <span class="special"><<</span> <span class="identifier">result</span> <span class="special"><<</span> <span class="string">", relative error estimate "</span> <span class="special"><<</span> <span class="identifier">relative_error</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
|
|
</pre>
|
|
<p>
|
|
The value of this integral should be π/(2e) and can be shown :
|
|
</p>
|
|
<pre class="programlisting"><span class="keyword">constexpr</span> <span class="keyword">double</span> <span class="identifier">expected</span> <span class="special">=</span> <span class="identifier">half_pi</span><span class="special"><</span><span class="keyword">double</span><span class="special">>()</span> <span class="special">/</span> <span class="identifier">e</span><span class="special"><</span><span class="keyword">double</span><span class="special">>();</span>
|
|
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"pi/(2e) = "</span> <span class="special"><<</span> <span class="identifier">expected</span> <span class="special"><<</span> <span class="string">", difference "</span> <span class="special"><<</span> <span class="identifier">result</span> <span class="special">-</span> <span class="identifier">expected</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
|
|
</pre>
|
|
<p>
|
|
or with the macro BOOST_MATH_INSTRUMENT_OOURA defined, we can follow the progress:
|
|
</p>
|
|
<pre class="programlisting">
|
|
ooura_fourier_cos with relative error goal 1.4901161193847656e-08 & 8 levels.
|
|
epsilon for type = 2.2204460492503131e-16
|
|
h = 1.000000000000000, I_h = 0.588268622591776 = 0x1.2d318b7e96dbe00p-1, absolute error estimate = nan
|
|
h = 0.500000000000000, I_h = 0.577871642184837 = 0x1.27decab8f07b200p-1, absolute error estimate = 1.039698040693926e-02
|
|
h = 0.250000000000000, I_h = 0.577863671186883 = 0x1.27ddbf42969be00p-1, absolute error estimate = 7.970997954576120e-06
|
|
h = 0.125000000000000, I_h = 0.577863674895461 = 0x1.27ddbf6271dc000p-1, absolute error estimate = 3.708578555361441e-09
|
|
Integral = 5.778636748954611e-01, relative error estimate 6.417739540441515e-09
|
|
pi/(2e) = 5.778636748954609e-01, difference 2.220446049250313e-16
|
|
|
|
</pre>
|
|
<p>
|
|
Working code of this example is at <a href="../../../example/ooura_fourier_integrals_cosine_example.cpp" target="_top">ooura_fourier_integrals_consine_example.cpp</a>
|
|
</p>
|
|
<h6>
|
|
<a name="math_toolkit.fourier_integrals.h1"></a>
|
|
<span class="phrase"><a name="math_toolkit.fourier_integrals.performance"></a></span><a class="link" href="fourier_integrals.html#math_toolkit.fourier_integrals.performance">Performance</a>
|
|
</h6>
|
|
<p>
|
|
The integrator precomputes nodes and weights, and hence can be reused for many
|
|
different frequencies with good efficiency. The integrator is pimpl'd and hence
|
|
can be shared between threads without a <code class="computeroutput"><span class="identifier">memcpy</span></code>
|
|
of the nodes and weights.
|
|
</p>
|
|
<p>
|
|
Ooura and Mori's paper identifies criteria for rapid convergence based on the
|
|
position of the poles of the integrand in the complex plane. If these poles
|
|
are too close to the real axis the convergence is slow. It is not trivial to
|
|
predict the convergence rate a priori, so if you are interested in figuring
|
|
out if the convergence is rapid, compile with <code class="computeroutput"><span class="special">-</span><span class="identifier">DBOOST_MATH_INSTRUMENT_OOURA</span></code> and some amount
|
|
of printing will give you a good idea of how well this method is performing.
|
|
</p>
|
|
<h6>
|
|
<a name="math_toolkit.fourier_integrals.h2"></a>
|
|
<span class="phrase"><a name="math_toolkit.fourier_integrals.multi_precision"></a></span><a class="link" href="fourier_integrals.html#math_toolkit.fourier_integrals.multi_precision">Higher
|
|
precision</a>
|
|
</h6>
|
|
<p>
|
|
It is simple to extend to higher precision using <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a>.
|
|
</p>
|
|
<pre class="programlisting"><span class="comment">// Use the default parameters for tolerance root_epsilon and eight levels for a type of 8 bytes.</span>
|
|
<span class="comment">//auto integrator = ooura_fourier_cos<Real>();</span>
|
|
<span class="comment">// Decide on a (tight) tolerance.</span>
|
|
<span class="keyword">const</span> <span class="identifier">Real</span> <span class="identifier">tol</span> <span class="special">=</span> <span class="number">2</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="identifier">Real</span><span class="special">>::</span><span class="identifier">epsilon</span><span class="special">();</span>
|
|
<span class="keyword">auto</span> <span class="identifier">integrator</span> <span class="special">=</span> <span class="identifier">ooura_fourier_cos</span><span class="special"><</span><span class="identifier">Real</span><span class="special">>(</span><span class="identifier">tol</span><span class="special">,</span> <span class="number">8</span><span class="special">);</span> <span class="comment">// Loops or gets worse for more than 8.</span>
|
|
|
|
<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="comment">// More complex example function.</span>
|
|
<span class="keyword">return</span> <span class="number">1</span> <span class="special">/</span> <span class="special">(</span><span class="identifier">x</span> <span class="special">*</span> <span class="identifier">x</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="identifier">omega</span> <span class="special">=</span> <span class="number">1</span><span class="special">;</span>
|
|
<span class="keyword">auto</span> <span class="special">[</span><span class="identifier">result</span><span class="special">,</span> <span class="identifier">relative_error</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">integrator</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="identifier">omega</span><span class="special">);</span>
|
|
</pre>
|
|
<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Integral = "</span> <span class="special"><<</span> <span class="identifier">result</span> <span class="special"><<</span> <span class="string">", relative error estimate "</span> <span class="special"><<</span> <span class="identifier">relative_error</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
|
|
|
|
<span class="keyword">const</span> <span class="identifier">Real</span> <span class="identifier">expected</span> <span class="special">=</span> <span class="identifier">half_pi</span><span class="special"><</span><span class="identifier">Real</span><span class="special">>()</span> <span class="special">/</span> <span class="identifier">e</span><span class="special"><</span><span class="identifier">Real</span><span class="special">>();</span> <span class="comment">// Expect integral = 1/(2e)</span>
|
|
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"pi/(2e) = "</span> <span class="special"><<</span> <span class="identifier">expected</span> <span class="special"><<</span> <span class="string">", difference "</span> <span class="special"><<</span> <span class="identifier">result</span> <span class="special">-</span> <span class="identifier">expected</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
|
|
</pre>
|
|
<p>
|
|
with output:
|
|
</p>
|
|
<pre class="programlisting">
|
|
Integral = 0.5778636748954608589550465916563501587, relative error estimate 4.609814684522163895264277312610830278e-17
|
|
pi/(2e) = 0.5778636748954608659545328919193707407, difference -6.999486300263020581921171645255733758e-18
|
|
|
|
</pre>
|
|
<p>
|
|
And with diagnostics on:
|
|
</p>
|
|
<pre class="programlisting">
|
|
ooura_fourier_cos with relative error goal 3.851859888774471706111955885169854637e-34 & 15 levels.
|
|
epsilon for type = 1.925929944387235853055977942584927319e-34
|
|
h = 1.000000000000000000000000000000000, I_h = 0.588268622591776615359568690603776 = 0.5882686225917766153595686906037760, absolute error estimate = nan
|
|
h = 0.500000000000000000000000000000000, I_h = 0.577871642184837461311756940493259 = 0.5778716421848374613117569404932595, absolute error estimate = 1.039698040693915404781175011051656e-02
|
|
h = 0.250000000000000000000000000000000, I_h = 0.577863671186882539559996800783122 = 0.5778636711868825395599968007831220, absolute error estimate = 7.970997954921751760139710137450075e-06
|
|
h = 0.125000000000000000000000000000000, I_h = 0.577863674895460885593491133506723 = 0.5778636748954608855934911335067232, absolute error estimate = 3.708578346033494332723601147051768e-09
|
|
h = 0.062500000000000000000000000000000, I_h = 0.577863674895460858955046591656350 = 0.5778636748954608589550465916563502, absolute error estimate = 2.663844454185037302771663314961535e-17
|
|
h = 0.031250000000000000000000000000000, I_h = 0.577863674895460858955046591656348 = 0.5778636748954608589550465916563484, absolute error estimate = 1.733336949948512267750380148326435e-33
|
|
h = 0.015625000000000000000000000000000, I_h = 0.577863674895460858955046591656348 = 0.5778636748954608589550465916563479, absolute error estimate = 4.814824860968089632639944856462318e-34
|
|
h = 0.007812500000000000000000000000000, I_h = 0.577863674895460858955046591656347 = 0.5778636748954608589550465916563473, absolute error estimate = 6.740754805355325485695922799047246e-34
|
|
h = 0.003906250000000000000000000000000, I_h = 0.577863674895460858955046591656347 = 0.5778636748954608589550465916563475, absolute error estimate = 1.925929944387235853055977942584927e-34
|
|
Integral = 5.778636748954608589550465916563475e-01, relative error estimate 3.332844800697411177051445985473052e-34
|
|
pi/(2e) = 5.778636748954608589550465916563481e-01, difference -6.740754805355325485695922799047246e-34
|
|
|
|
</pre>
|
|
<p>
|
|
Working code of this example is at <a href="../../../example/ooura_fourier_integrals_multiprecision_example.cpp" target="_top">ooura_fourier_integrals_multiprecision_example.cpp</a>
|
|
</p>
|
|
<p>
|
|
For more examples of other functions and tests, see the full test suite at
|
|
<a href="../../../test/ooura_fourier_integral_test.cpp" target="_top">ooura_fourier_integral_test.cpp</a>.
|
|
</p>
|
|
<p>
|
|
Ngyen and Nuyens make use of <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a>
|
|
in their extension to multiple dimensions, showing relative errors reducing
|
|
to ≅ 10<sup>-2000</sup>!
|
|
</p>
|
|
<h6>
|
|
<a name="math_toolkit.fourier_integrals.h3"></a>
|
|
<span class="phrase"><a name="math_toolkit.fourier_integrals.rationale"></a></span><a class="link" href="fourier_integrals.html#math_toolkit.fourier_integrals.rationale">Rationale</a>
|
|
</h6>
|
|
<p>
|
|
This implementation is base on Ooura's 1999 paper rather than the later 2005
|
|
paper. It does not preclude a second future implementation based on the later
|
|
work.
|
|
</p>
|
|
<p>
|
|
Some of the features of the Ooura's 2005 paper that were less appealing were:
|
|
</p>
|
|
<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
|
|
<li class="listitem">
|
|
The advance of that paper is that one can compute <span class="emphasis"><em>both</em></span>
|
|
the Fourier sine transform and Fourier cosine transform in a single shot.
|
|
But there are examples, like sinc integral, where the Fourier sine would
|
|
converge, but the Fourier cosine would diverge.
|
|
</li>
|
|
<li class="listitem">
|
|
It would force users to live in the complex plane, when many potential
|
|
applications only need real.
|
|
</li>
|
|
</ul></div>
|
|
<h5>
|
|
<a name="math_toolkit.fourier_integrals.h4"></a>
|
|
<span class="phrase"><a name="math_toolkit.fourier_integrals.references"></a></span><a class="link" href="fourier_integrals.html#math_toolkit.fourier_integrals.references">References</a>
|
|
</h5>
|
|
<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
|
|
<li class="listitem">
|
|
Ooura, Takuya, and Masatake Mori, <span class="emphasis"><em>A robust double exponential
|
|
formula for Fourier-type integrals.</em></span> Journal of computational
|
|
and applied mathematics, 112.1-2 (1999): 229-241.
|
|
</li>
|
|
<li class="listitem">
|
|
Ooura, Takuya, <span class="emphasis"><em>A Double Exponential Formula for the Fourier Transforms.</em></span>
|
|
Publ. RIMS, Kyoto Univ., 41 (2005), 971-977. <a href="https://pdfs.semanticscholar.org/16ec/a5d76fd6b3d7acaaff0b2a6e8a70caa70190.pdf" target="_top">https://pdfs.semanticscholar.org/16ec/a5d76fd6b3d7acaaff0b2a6e8a70caa70190.pdf</a>
|
|
</li>
|
|
<li class="listitem">
|
|
Khatibi, Arezoo and Khatibi, Omid,<span class="emphasis"><em>Criteria for the Application
|
|
of Double Exponential Transformation.</em></span> (2017) <a href="https://arxiv.org/pdf/1704.05752.pdf" target="_top">1704.05752.pdf</a>.
|
|
</li>
|
|
<li class="listitem">
|
|
Nguyen, Dong T.P. and Nuyens, Dirk, <span class="emphasis"><em>Multivariate integration
|
|
over Reals with exponential rate of convergence.</em></span> (2016) <a href="https://core.ac.uk/download/pdf/80799199.pdf" target="_top">https://core.ac.uk/download/pdf/80799199.pdf</a>.
|
|
</li>
|
|
</ul></div>
|
|
</div>
|
|
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
|
|
<td align="left"></td>
|
|
<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar
|
|
Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos,
|
|
Hubert Holin, Bruno Lalande, John Maddock, 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></td>
|
|
</tr></table>
|
|
<hr>
|
|
<div class="spirit-nav">
|
|
<a accesskey="p" href="double_exponential/de_refes.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../quadrature.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="naive_monte_carlo.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a>
|
|
</div>
|
|
</body>
|
|
</html>
|