123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317 |
- <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>
|