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
|
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<title>Trapezoidal Quadrature</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="../quadrature.html" title="Chapter 14. Quadrature and Differentiation">
<link rel="next" href="gauss.html" title="Gauss-Legendre quadrature">
<meta name="viewport" content="width=device-width, initial-scale=1">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td>
<td align="center"><a href="../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="../quadrature.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="gauss.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.trapezoidal"></a><a class="link" href="trapezoidal.html" title="Trapezoidal Quadrature">Trapezoidal Quadrature</a>
</h2></div></div></div>
<h4>
<a name="math_toolkit.trapezoidal.h0"></a>
<span class="phrase"><a name="math_toolkit.trapezoidal.synopsis"></a></span><a class="link" href="trapezoidal.html#math_toolkit.trapezoidal.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">trapezoidal</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">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Real</span><span class="special">></span>
<span class="keyword">auto</span> <span class="identifier">trapezoidal</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">a</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">b</span><span class="special">,</span>
<span class="identifier">Real</span> <span class="identifier">tol</span> <span class="special">=</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">Real</span><span class="special">>::</span><span class="identifier">epsilon</span><span class="special">()),</span>
<span class="identifier">size_t</span> <span class="identifier">max_refinements</span> <span class="special">=</span> <span class="number">12</span><span class="special">,</span>
<span class="identifier">Real</span><span class="special">*</span> <span class="identifier">error_estimate</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">,</span>
<span class="identifier">Real</span><span class="special">*</span> <span class="identifier">L1</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">);</span>
<span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Real</span><span class="special">,</span> <span class="keyword">class</span> <a class="link" href="../policy.html" title="Chapter 22. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span>
<span class="keyword">auto</span> <span class="identifier">trapezoidal</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">a</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">b</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">tol</span><span class="special">,</span> <span class="identifier">size_t</span> <span class="identifier">max_refinements</span><span class="special">,</span>
<span class="identifier">Real</span><span class="special">*</span> <span class="identifier">error_estimate</span><span class="special">,</span> <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">L1</span><span class="special">,</span> <span class="keyword">const</span> <a class="link" href="../policy.html" title="Chapter 22. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&</span> <span class="identifier">pol</span><span class="special">);</span>
<span class="special">}}}</span> <span class="comment">// namespaces</span>
</pre>
<h4>
<a name="math_toolkit.trapezoidal.h1"></a>
<span class="phrase"><a name="math_toolkit.trapezoidal.description"></a></span><a class="link" href="trapezoidal.html#math_toolkit.trapezoidal.description">Description</a>
</h4>
<p>
The functional <code class="computeroutput"><span class="identifier">trapezoidal</span></code>
calculates the integral of a function <span class="emphasis"><em>f</em></span> using the surprisingly
simple trapezoidal rule. If we assume only that the integrand is twice continuously
differentiable, we can prove that the error of the composite trapezoidal rule
is 𝑶(h<sup>2</sup>). Hence halving the interval only cuts the error by about a fourth,
which in turn implies that we must evaluate the function many times before
an acceptable accuracy can be achieved.
</p>
<p>
However, the trapezoidal rule has an astonishing property: If the integrand
is periodic, and we integrate it over a period, then the trapezoidal rule converges
faster than any power of the step size <span class="emphasis"><em>h</em></span>. This can be
seen by examination of the <a href="https://en.wikipedia.org/wiki/Euler-Maclaurin_formula" target="_top">Euler-Maclaurin
summation formula</a>, which relates a definite integral to its trapezoidal
sum and error terms proportional to the derivatives of the function at the
endpoints and the Bernoulli numbers. If the derivatives at the endpoints are
the same or vanish, then the error very nearly vanishes. Hence the trapezoidal
rule is essentially optimal for periodic integrands.
</p>
<p>
Other classes of integrands which are integrated efficiently by this method
are the C<sub>0</sub><sup>∞</sup>(∝) <a href="https://en.wikipedia.org/wiki/Bump_function" target="_top">bump
functions</a> and bell-shaped integrals over the infinite interval. For
details, see <a href="http://epubs.siam.org/doi/pdf/10.1137/130932132" target="_top">Trefethen's</a>
SIAM review.
</p>
<p>
In its simplest form, an integration can be performed by the following code
</p>
<pre class="programlisting"><span class="keyword">using</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">trapezoidal</span><span class="special">;</span>
<span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="number">1</span><span class="special">/(</span><span class="number">5</span> <span class="special">-</span> <span class="number">4</span><span class="special">*</span><span class="identifier">cos</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">I</span> <span class="special">=</span> <span class="identifier">trapezoidal</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0.0</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">constants</span><span class="special">::</span><span class="identifier">two_pi</span><span class="special"><</span><span class="keyword">double</span><span class="special">>());</span>
</pre>
<p>
The integrand must accept a real number argument, but can return a complex
number. This is useful for contour integrals (which are manifestly periodic)
and high-order numerical differentiation of analytic functions. An example
using the integral definition of the complex Bessel function is shown here:
</p>
<pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">bessel_integrand</span> <span class="special">=</span> <span class="special">[&</span><span class="identifier">n</span><span class="special">,</span> <span class="special">&</span><span class="identifier">z</span><span class="special">](</span><span class="keyword">double</span> <span class="identifier">theta</span><span class="special">)-></span><span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span>
<span class="special">{</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">z</span><span class="special">{</span><span class="number">2</span><span class="special">,</span> <span class="number">3</span><span class="special">};</span>
<span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cos</span><span class="special">;</span>
<span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">sin</span><span class="special">;</span>
<span class="keyword">return</span> <span class="identifier">cos</span><span class="special">(</span><span class="identifier">z</span><span class="special">*</span><span class="identifier">sin</span><span class="special">(</span><span class="identifier">theta</span><span class="special">)</span> <span class="special">-</span> <span class="number">2</span><span class="special">*</span><span class="identifier">theta</span><span class="special">)/</span><span class="identifier">pi</span><span class="special"><</span><span class="keyword">double</span><span class="special">>();</span>
<span class="special">};</span>
<span class="keyword">using</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">trapezoidal</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">Jnz</span> <span class="special">=</span> <span class="identifier">trapezoidal</span><span class="special">(</span><span class="identifier">bessel_integrand</span><span class="special">,</span> <span class="number">0.0</span><span class="special">,</span> <span class="identifier">pi</span><span class="special"><</span><span class="identifier">Real</span><span class="special">>());</span>
</pre>
<p>
Other special functions which are efficiently evaluated in the complex plane
by trapezoidal quadrature are modified Bessel functions and the complementary
error function. Another application of complex-valued trapezoidal quadrature
is computation of high-order numerical derivatives; see Lyness and Moler for
details.
</p>
<p>
Since the routine is adaptive, step sizes are halved continuously until a tolerance
is reached. In order to control this tolerance, simply call the routine with
an additional argument
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">I</span> <span class="special">=</span> <span class="identifier">trapezoidal</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0.0</span><span class="special">,</span> <span class="identifier">two_pi</span><span class="special"><</span><span class="keyword">double</span><span class="special">>(),</span> <span class="number">1e-6</span><span class="special">);</span>
</pre>
<p>
The routine stops when successive estimates of the integral <code class="computeroutput"><span class="identifier">I1</span></code>
and <code class="computeroutput"><span class="identifier">I0</span></code> differ by less than
the tolerance multiplied by the estimated L<sub>1</sub> norm of the function. A good choice
for the tolerance is √ε, which is the default. If the integrand is periodic,
then the number of correct digits should double on each interval halving. Hence,
once the integration routine has estimated that the error is √ε, then the actual
error should be ~ε. If the integrand is <span class="bold"><strong>not</strong></span>
periodic, then reducing the error to √ε takes much longer, but is nonetheless
possible without becoming a major performance bug.
</p>
<p>
A question arises as to what to do when successive estimates never pass below
the tolerance threshold. The stepsize would be halved repeatedly, generating
an exponential explosion in function evaluations. As such, you may pass an
optional argument <code class="computeroutput"><span class="identifier">max_refinements</span></code>
which controls how many times the interval may be halved before giving up.
By default, this maximum number of refinement steps is 12, which should never
be hit in double precision unless the function is not periodic. However, for
higher-precision types, it may be of interest to allow the algorithm to compute
more refinements:
</p>
<pre class="programlisting"><span class="identifier">size_t</span> <span class="identifier">max_refinements</span> <span class="special">=</span> <span class="number">15</span><span class="special">;</span>
<span class="keyword">long</span> <span class="keyword">double</span> <span class="identifier">I</span> <span class="special">=</span> <span class="identifier">trapezoidal</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0.0L</span><span class="special">,</span> <span class="identifier">two_pi</span><span class="special"><</span><span class="keyword">long</span> <span class="keyword">double</span><span class="special">>(),</span> <span class="number">1e-9L</span><span class="special">,</span> <span class="identifier">max_refinements</span><span class="special">);</span>
</pre>
<p>
Note that the maximum allowed compute time grows exponentially with <code class="computeroutput"><span class="identifier">max_refinements</span></code>. The routine will not throw
an exception if the maximum refinements is achieved without the requested tolerance
being attained. This is because the value calculated is more often than not
still usable. However, for applications with high-reliability requirements,
the error estimate should be queried. This is achieved by passing additional
pointers into the routine:
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">error_estimate</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">L1</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">I</span> <span class="special">=</span> <span class="identifier">trapezoidal</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0.0</span><span class="special">,</span> <span class="identifier">two_pi</span><span class="special"><</span><span class="keyword">double</span><span class="special">>(),</span> <span class="identifier">tolerance</span><span class="special">,</span> <span class="identifier">max_refinements</span><span class="special">,</span> <span class="special">&</span><span class="identifier">error_estimate</span><span class="special">,</span> <span class="special">&</span><span class="identifier">L1</span><span class="special">);</span>
<span class="keyword">if</span> <span class="special">(</span><span class="identifier">error_estimate</span> <span class="special">></span> <span class="identifier">tolerance</span><span class="special">*</span><span class="identifier">L1</span><span class="special">)</span>
<span class="special">{</span>
<span class="keyword">double</span> <span class="identifier">I</span> <span class="special">=</span> <span class="identifier">some_other_quadrature_method</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0</span><span class="special">,</span> <span class="identifier">two_pi</span><span class="special"><</span><span class="keyword">double</span><span class="special">>());</span>
<span class="special">}</span>
</pre>
<p>
The final argument is the L<sub>1</sub> norm of the integral. This is computed along with
the integral, and is an essential component of the algorithm. First, the L<sub>1</sub> norm
establishes a scale against which the error can be measured. Second, the L<sub>1</sub> norm
can be used to evaluate the stability of the computation. This can be formulated
in a rigorous manner by defining the <span class="bold"><strong>condition number
of summation</strong></span>. The condition number of summation is defined by
</p>
<div class="blockquote"><blockquote class="blockquote"><p>
<span class="serif_italic"><span class="emphasis"><em>κ(S<sub>n</sub>) := Σ<sub>i</sub><sup>n</sup> |x<sub>i</sub>|/|Σ<sub>i</sub><sup>n</sup> x<sub>i</sub>|</em></span></span>
</p></blockquote></div>
<p>
If this number of ~10<sup>k</sup>, then <span class="emphasis"><em>k</em></span> additional digits are expected
to be lost in addition to digits lost due to floating point rounding error.
As all numerical quadrature methods reduce to summation, their stability is
controlled by the ratio ∫ |f| dt/|∫ f dt |, which is easily seen
to be equivalent to condition number of summation when evaluated numerically.
Hence both the error estimate and the condition number of summation should
be analyzed in applications requiring very high precision and reliability.
</p>
<p>
As an example, we consider evaluation of Bessel functions by trapezoidal quadrature.
The Bessel function of the first kind is defined via
</p>
<div class="blockquote"><blockquote class="blockquote"><p>
<span class="serif_italic"><span class="emphasis"><em>J<sub>n</sub>(x) = 1/2Π ∫<sub>-Π</sub><sup>Π</sup> cos(n
t - x sin(t)) dt</em></span></span>
</p></blockquote></div>
<p>
The integrand is periodic, so the Euler-Maclaurin summation formula guarantees
exponential convergence via the trapezoidal quadrature. Without careful consideration,
it seems this would be a very attractive method to compute Bessel functions.
However, we see that for large <span class="emphasis"><em>n</em></span> the integrand oscillates
rapidly, taking on positive and negative values, and hence the trapezoidal
sums become ill-conditioned. In double precision, <span class="emphasis"><em>x = 17</em></span>
and <span class="emphasis"><em>n = 25</em></span> gives a sum which is so poorly conditioned
that zero correct digits are obtained.
</p>
<p>
The final <a class="link" href="../policy.html" title="Chapter 22. Policies: Controlling Precision, Error Handling etc">Policy</a> argument is optional and can
be used to control the behaviour of the function: how it handles errors, what
level of precision to use etc. Refer to the <a class="link" href="../policy.html" title="Chapter 22. Policies: Controlling Precision, Error Handling etc">policy documentation
for more details</a>.
</p>
<p>
References:
</p>
<p>
Trefethen, Lloyd N., Weideman, J.A.C., <span class="emphasis"><em>The Exponentially Convergent
Trapezoidal Rule</em></span>, SIAM Review, Vol. 56, No. 3, 2014.
</p>
<p>
Stoer, Josef, and Roland Bulirsch. <span class="emphasis"><em>Introduction to numerical analysis.
Vol. 12.</em></span>, Springer Science & Business Media, 2013.
</p>
<p>
Higham, Nicholas J. <span class="emphasis"><em>Accuracy and stability of numerical algorithms.</em></span>
Society for industrial and applied mathematics, 2002.
</p>
<p>
Lyness, James N., and Cleve B. Moler. <span class="emphasis"><em>Numerical differentiation of
analytic functions.</em></span> SIAM Journal on Numerical Analysis 4.2 (1967):
202-210.
</p>
<p>
Gil, Amparo, Javier Segura, and Nico M. Temme. <span class="emphasis"><em>Computing special
functions by using quadrature rules.</em></span> Numerical Algorithms 33.1-4
(2003): 265-275.
</p>
</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="../quadrature.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="gauss.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>
|
