libs/math/doc/html/math_toolkit/fourier_integrals.html

   1 <html>
   2 <head>
   3 <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
   4 <title>Fourier Integrals</title>
   5 <link rel="stylesheet" href="../math.css" type="text/css">
   6 <meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
   7 <link rel="home" href="../index.html" title="Math Toolkit 2.11.0">
   8 <link rel="up" href="../quadrature.html" title="Chapter&#160;13.&#160;Quadrature and Differentiation">
   9 <link rel="prev" href="double_exponential/de_refes.html" title="References">
  10 <link rel="next" href="naive_monte_carlo.html" title="Naive Monte Carlo Integration">
  11 </head>
  12 <body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
  13 <table cellpadding="2" width="100%"><tr>
  14 <td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td>
  15 <td align="center"><a href="../../../../../index.html">Home</a></td>
  16 <td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td>
  17 <td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
  18 <td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
  19 <td align="center"><a href="../../../../../more/index.htm">More</a></td>
  20 </tr></table>
  21 <hr>
  22 <div class="spirit-nav">
  23 <a accesskey="p" href="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>
  24 </div>
  25 <div class="section">
  26 <div class="titlepage"><div><div><h2 class="title" style="clear: both">
  27 <a name="math_toolkit.fourier_integrals"></a><a class="link" href="fourier_integrals.html" title="Fourier Integrals">Fourier Integrals</a>
  28 </h2></div></div></div>
  29 <h4>
  30 <a name="math_toolkit.fourier_integrals.h0"></a>
  31       <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>
  32     </h4>
  33 <pre class="programlisting"><span class="preprocessor">#include</span> <span class="special">&lt;</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">&gt;</span>
  34
  35 <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>
  36
  37 <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>
  38 <span class="keyword">class</span> <span class="identifier">ooura_fourier_sin</span> <span class="special">{</span>
  39 <span class="keyword">public</span><span class="special">:</span>
  40     <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">&lt;</span><span class="identifier">Real</span><span class="special">&gt;(),</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>
  41
  42     <span class="keyword">template</span><span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">&gt;</span>
  43     <span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">,</span> <span class="identifier">Real</span><span class="special">&gt;</span> <span class="identifier">integrate</span><span class="special">(</span><span class="identifier">F</span> <span class="keyword">const</span> <span class="special">&amp;</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">omega</span><span class="special">);</span>
  44
  45 <span class="special">};</span>
  46
  47
  48 <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>
  49 <span class="keyword">class</span> <span class="identifier">ooura_fourier_cos</span> <span class="special">{</span>
  50 <span class="keyword">public</span><span class="special">:</span>
  51     <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">&lt;</span><span class="identifier">Real</span><span class="special">&gt;(),</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>
  52
  53     <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>
  54     <span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">,</span> <span class="identifier">Real</span><span class="special">&gt;</span> <span class="identifier">integrate</span><span class="special">(</span><span class="identifier">F</span> <span class="keyword">const</span> <span class="special">&amp;</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">omega</span><span class="special">);</span>
  55 <span class="special">};</span>
  56
  57 <span class="special">}}}</span> <span class="comment">// namespaces</span>
  58 </pre>
  59 <p>
  60       Ooura's method for Fourier integrals computes
  61     </p>
  62 <div class="blockquote"><blockquote class="blockquote"><p>
  63         <span class="serif_italic">&#8747;<sub>0</sub><sup>&#8734;</sup> f(t)sin(&#969; t) dt</span>
  64       </p></blockquote></div>
  65 <p>
  66       and
  67     </p>
  68 <div class="blockquote"><blockquote class="blockquote"><p>
  69         <span class="serif_italic">&#8747;<sub>0</sub><sup>&#8734;</sup> f(t)cos(&#969; t) dt</span>
  70       </p></blockquote></div>
  71 <p>
  72       by a double exponentially decaying transformation. These integrals arise when
  73       computing continuous Fourier transform of odd and even functions, respectively.
  74       Oscillatory integrals are known to cause trouble for standard quadrature methods,
  75       so these routines are provided to cope with the most common oscillatory use
  76       case.
  77     </p>
  78 <p>
  79       The basic usage is shown below:
  80     </p>
  81 <pre class="programlisting"><span class="identifier">ooura_fourier_sin</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="identifier">ooura_fourier_sin</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;();</span>
  82 <span class="comment">// Use the default tolerance root_epsilon and eight levels for type double.</span>
  83
  84 <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>
  85 <span class="special">{</span> <span class="comment">// Simple reciprocal function for sinc.</span>
  86   <span class="keyword">return</span> <span class="number">1</span> <span class="special">/</span> <span class="identifier">x</span><span class="special">;</span>
  87 <span class="special">};</span>
  88
  89 <span class="keyword">double</span> <span class="identifier">omega</span> <span class="special">=</span> <span class="number">1</span><span class="special">;</span>
  90 <span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</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>
  91 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Integral = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">result</span><span class="special">.</span><span class="identifier">first</span> <span class="special">&lt;&lt;</span> <span class="string">", relative error estimate "</span> <span class="special">&lt;&lt;</span> <span class="identifier">result</span><span class="special">.</span><span class="identifier">second</span> <span class="special">&lt;&lt;</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
  92 </pre>
  93 <p>
  94       and compare with the expected value &#960;/2 of the integral.
  95     </p>
  96 <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">&lt;</span><span class="keyword">double</span><span class="special">&gt;();</span>
  97 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"pi/2 =     "</span> <span class="special">&lt;&lt;</span> <span class="identifier">expected</span> <span class="special">&lt;&lt;</span> <span class="string">", difference "</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
  98 </pre>
  99 <p>
 100       The output is
 101     </p>
 102 <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>
 103 <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>
 104 </pre>
 105 <div class="note"><table border="0" summary="Note">
 106 <tr>
 107 <td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../../doc/src/images/note.png"></td>
 108 <th align="left">Note</th>
 109 </tr>
 110 <tr><td align="left" valign="top"><p>
 111         This integrator is more insistent about examining the error estimate, than
 112         (say) tanh-sinh, which just returns the value of the integral.
 113       </p></td></tr>
 114 </table></div>
 115 <p>
 116       With the macro BOOST_MATH_INSTRUMENT_OOURA defined, we can follow the progress:
 117     </p>
 118 <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">&amp;</span> <span class="number">8</span> <span class="identifier">levels</span><span class="special">.</span>
 119 <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>
 120 <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>
 121 <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>
 122 <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>
 123 <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>
 124 <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>
 125 </pre>
 126 <p>
 127       Working code of this example is at <a href="../../../example/ooura_fourier_integrals_example.cpp" target="_top">ooura_fourier_integrals_example.cpp</a>
 128     </p>
 129 <p>
 130       A classical cosine transform is presented below:
 131     </p>
 132 <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">&lt;</span><span class="keyword">double</span><span class="special">&gt;();</span>
 133 <span class="comment">// Use the default tolerance root_epsilon and eight levels for type double.</span>
 134
 135 <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>
 136 <span class="special">{</span> <span class="comment">// More complex example function.</span>
 137   <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>
 138 <span class="special">};</span>
 139
 140 <span class="keyword">double</span> <span class="identifier">omega</span> <span class="special">=</span> <span class="number">1</span><span class="special">;</span>
 141
 142 <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>
 143 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Integral = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">result</span> <span class="special">&lt;&lt;</span> <span class="string">", relative error estimate "</span> <span class="special">&lt;&lt;</span> <span class="identifier">relative_error</span> <span class="special">&lt;&lt;</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
 144 </pre>
 145 <p>
 146       The value of this integral should be &#960;/(2e) and can be shown :
 147     </p>
 148 <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">&lt;</span><span class="keyword">double</span><span class="special">&gt;()</span> <span class="special">/</span> <span class="identifier">e</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;();</span>
 149 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"pi/(2e) =  "</span> <span class="special">&lt;&lt;</span> <span class="identifier">expected</span> <span class="special">&lt;&lt;</span> <span class="string">", difference "</span> <span class="special">&lt;&lt;</span> <span class="identifier">result</span> <span class="special">-</span> <span class="identifier">expected</span> <span class="special">&lt;&lt;</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
 150 </pre>
 151 <p>
 152       or with the macro BOOST_MATH_INSTRUMENT_OOURA defined, we can follow the progress:
 153     </p>
 154 <pre class="programlisting">
 155 ooura_fourier_cos with relative error goal 1.4901161193847656e-08 &amp; 8 levels.
 156 epsilon for type = 2.2204460492503131e-16
 157 h = 1.000000000000000, I_h = 0.588268622591776 = 0x1.2d318b7e96dbe00p-1, absolute error estimate = nan
 158 h = 0.500000000000000, I_h = 0.577871642184837 = 0x1.27decab8f07b200p-1, absolute error estimate = 1.039698040693926e-02
 159 h = 0.250000000000000, I_h = 0.577863671186883 = 0x1.27ddbf42969be00p-1, absolute error estimate = 7.970997954576120e-06
 160 h = 0.125000000000000, I_h = 0.577863674895461 = 0x1.27ddbf6271dc000p-1, absolute error estimate = 3.708578555361441e-09
 161 Integral = 5.778636748954611e-01, relative error estimate 6.417739540441515e-09
 162 pi/(2e)  = 5.778636748954609e-01, difference 2.220446049250313e-16
 163
 164 </pre>
 165 <p>
 166       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>
 167     </p>
 168 <h6>
 169 <a name="math_toolkit.fourier_integrals.h1"></a>
 170       <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>
 171     </h6>
 172 <p>
 173       The integrator precomputes nodes and weights, and hence can be reused for many
 174       different frequencies with good efficiency. The integrator is pimpl'd and hence
 175       can be shared between threads without a <code class="computeroutput"><span class="identifier">memcpy</span></code>
 176       of the nodes and weights.
 177     </p>
 178 <p>
 179       Ooura and Mori's paper identifies criteria for rapid convergence based on the
 180       position of the poles of the integrand in the complex plane. If these poles
 181       are too close to the real axis the convergence is slow. It is not trivial to
 182       predict the convergence rate a priori, so if you are interested in figuring
 183       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
 184       of printing will give you a good idea of how well this method is performing.
 185     </p>
 186 <h6>
 187 <a name="math_toolkit.fourier_integrals.h2"></a>
 188       <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
 189       precision</a>
 190     </h6>
 191 <p>
 192       It is simple to extend to higher precision using <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a>.
 193     </p>
 194 <pre class="programlisting"><span class="comment">// Use the default parameters for tolerance root_epsilon and eight levels for a type of 8 bytes.</span>
 195 <span class="comment">//auto integrator = ooura_fourier_cos&lt;Real&gt;();</span>
 196 <span class="comment">// Decide on a (tight) tolerance.</span>
 197 <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">&lt;</span><span class="identifier">Real</span><span class="special">&gt;::</span><span class="identifier">epsilon</span><span class="special">();</span>
 198 <span class="keyword">auto</span> <span class="identifier">integrator</span> <span class="special">=</span> <span class="identifier">ooura_fourier_cos</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;(</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>
 199
 200 <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>
 201 <span class="special">{</span> <span class="comment">// More complex example function.</span>
 202   <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>
 203 <span class="special">};</span>
 204
 205 <span class="keyword">double</span> <span class="identifier">omega</span> <span class="special">=</span> <span class="number">1</span><span class="special">;</span>
 206 <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>
 207 </pre>
 208 <pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Integral = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">result</span> <span class="special">&lt;&lt;</span> <span class="string">", relative error estimate "</span> <span class="special">&lt;&lt;</span> <span class="identifier">relative_error</span> <span class="special">&lt;&lt;</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
 209
 210 <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">&lt;</span><span class="identifier">Real</span><span class="special">&gt;()</span> <span class="special">/</span> <span class="identifier">e</span><span class="special">&lt;</span><span class="identifier">Real</span><span class="special">&gt;();</span> <span class="comment">// Expect integral = 1/(2e)</span>
 211 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"pi/(2e)  = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">expected</span> <span class="special">&lt;&lt;</span> <span class="string">", difference "</span> <span class="special">&lt;&lt;</span> <span class="identifier">result</span> <span class="special">-</span> <span class="identifier">expected</span> <span class="special">&lt;&lt;</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
 212 </pre>
 213 <p>
 214       with output:
 215     </p>
 216 <pre class="programlisting">
 217 Integral = 0.5778636748954608589550465916563501587, relative error estimate 4.609814684522163895264277312610830278e-17
 218 pi/(2e) = 0.5778636748954608659545328919193707407, difference -6.999486300263020581921171645255733758e-18
 219
 220 </pre>
 221 <p>
 222       And with diagnostics on:
 223     </p>
 224 <pre class="programlisting">
 225 ooura_fourier_cos with relative error goal 3.851859888774471706111955885169854637e-34 &amp; 15 levels.
 226 epsilon for type = 1.925929944387235853055977942584927319e-34
 227 h = 1.000000000000000000000000000000000, I_h = 0.588268622591776615359568690603776 = 0.5882686225917766153595686906037760, absolute error estimate = nan
 228 h = 0.500000000000000000000000000000000, I_h = 0.577871642184837461311756940493259 = 0.5778716421848374613117569404932595, absolute error estimate = 1.039698040693915404781175011051656e-02
 229 h = 0.250000000000000000000000000000000, I_h = 0.577863671186882539559996800783122 = 0.5778636711868825395599968007831220, absolute error estimate = 7.970997954921751760139710137450075e-06
 230 h = 0.125000000000000000000000000000000, I_h = 0.577863674895460885593491133506723 = 0.5778636748954608855934911335067232, absolute error estimate = 3.708578346033494332723601147051768e-09
 231 h = 0.062500000000000000000000000000000, I_h = 0.577863674895460858955046591656350 = 0.5778636748954608589550465916563502, absolute error estimate = 2.663844454185037302771663314961535e-17
 232 h = 0.031250000000000000000000000000000, I_h = 0.577863674895460858955046591656348 = 0.5778636748954608589550465916563484, absolute error estimate = 1.733336949948512267750380148326435e-33
 233 h = 0.015625000000000000000000000000000, I_h = 0.577863674895460858955046591656348 = 0.5778636748954608589550465916563479, absolute error estimate = 4.814824860968089632639944856462318e-34
 234 h = 0.007812500000000000000000000000000, I_h = 0.577863674895460858955046591656347 = 0.5778636748954608589550465916563473, absolute error estimate = 6.740754805355325485695922799047246e-34
 235 h = 0.003906250000000000000000000000000, I_h = 0.577863674895460858955046591656347 = 0.5778636748954608589550465916563475, absolute error estimate = 1.925929944387235853055977942584927e-34
 236 Integral = 5.778636748954608589550465916563475e-01, relative error estimate 3.332844800697411177051445985473052e-34
 237 pi/(2e)  = 5.778636748954608589550465916563481e-01, difference -6.740754805355325485695922799047246e-34
 238
 239 </pre>
 240 <p>
 241       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>
 242     </p>
 243 <p>
 244       For more examples of other functions and tests, see the full test suite at
 245       <a href="../../../test/ooura_fourier_integral_test.cpp" target="_top">ooura_fourier_integral_test.cpp</a>.
 246     </p>
 247 <p>
 248       Ngyen and Nuyens make use of <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a>
 249       in their extension to multiple dimensions, showing relative errors reducing
 250       to &#8773; 10<sup>-2000</sup>!
 251     </p>
 252 <h6>
 253 <a name="math_toolkit.fourier_integrals.h3"></a>
 254       <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>
 255     </h6>
 256 <p>
 257       This implementation is base on Ooura's 1999 paper rather than the later 2005
 258       paper. It does not preclude a second future implementation based on the later
 259       work.
 260     </p>
 261 <p>
 262       Some of the features of the Ooura's 2005 paper that were less appealing were:
 263     </p>
 264 <div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
 265 <li class="listitem">
 266           The advance of that paper is that one can compute <span class="emphasis"><em>both</em></span>
 267           the Fourier sine transform and Fourier cosine transform in a single shot.
 268           But there are examples, like sinc integral, where the Fourier sine would
 269           converge, but the Fourier cosine would diverge.
 270         </li>
 271 <li class="listitem">
 272           It would force users to live in the complex plane, when many potential
 273           applications only need real.
 274         </li>
 275 </ul></div>
 276 <h5>
 277 <a name="math_toolkit.fourier_integrals.h4"></a>
 278       <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>
 279     </h5>
 280 <div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
 281 <li class="listitem">
 282           Ooura, Takuya, and Masatake Mori, <span class="emphasis"><em>A robust double exponential
 283           formula for Fourier-type integrals.</em></span> Journal of computational
 284           and applied mathematics, 112.1-2 (1999): 229-241.
 285         </li>
 286 <li class="listitem">
 287           Ooura, Takuya, <span class="emphasis"><em>A Double Exponential Formula for the Fourier Transforms.</em></span>
 288           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>
 289         </li>
 290 <li class="listitem">
 291           Khatibi, Arezoo and Khatibi, Omid,<span class="emphasis"><em>Criteria for the Application
 292           of Double Exponential Transformation.</em></span> (2017) <a href="https://arxiv.org/pdf/1704.05752.pdf" target="_top">1704.05752.pdf</a>.
 293         </li>
 294 <li class="listitem">
 295           Nguyen, Dong T.P. and Nuyens, Dirk, <span class="emphasis"><em>Multivariate integration
 296           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>.
 297         </li>
 298 </ul></div>
 299 </div>
 300 <table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
 301 <td align="left"></td>
 302 <td align="right"><div class="copyright-footer">Copyright &#169; 2006-2019 Nikhar
 303       Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos,
 304       Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan
 305       R&#229;de, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg,
 306       Daryle Walker and Xiaogang Zhang<p>
 307         Distributed under the Boost Software License, Version 1.0. (See accompanying
 308         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>)
 309       </p>
 310 </div></td>
 311 </tr></table>
 312 <hr>
 313 <div class="spirit-nav">
 314 <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>
 315 </div>
 316 </body>
 317 </html>