1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
|
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<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 4.2.1">
<link rel="up" href="../quadrature.html" title="Chapter 14. 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">
<meta name="viewport" content="width=device-width, initial-scale=1">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td>
<td align="center"><a href="../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="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>
<div class="copyright-footer">Copyright © 2006-2021 Nikhar Agrawal, Anton Bikineev, Matthew Borland,
Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno
Lalande, John Maddock, Evan Miller, Jeremy Murphy, Matthew Pulver, Johan Råde,
Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, Daryle
Walker and Xiaogang Zhang<p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
</p>
</div>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="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>
|
