123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152 |
- [/
- Copyright (c) 2017 Nick Thompson
- Use, modification and distribution are subject to the
- Boost Software License, Version 1.0. (See accompanying file
- LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
- ]
- [section:trapezoidal Trapezoidal Quadrature]
- [heading Synopsis]
- #include <boost/math/quadrature/trapezoidal.hpp>
- namespace boost{ namespace math{ namespace quadrature {
- template<class F, class Real>
- auto trapezoidal(F f, Real a, Real b,
- Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
- size_t max_refinements = 12,
- Real* error_estimate = nullptr,
- Real* L1 = nullptr);
- template<class F, class Real, class ``__Policy``>
- auto trapezoidal(F f, Real a, Real b, Real tol, size_t max_refinements,
- Real* error_estimate, Real* L1, const ``__Policy``& pol);
- }}} // namespaces
- [heading Description]
- The functional `trapezoidal` calculates the integral of a function /f/ 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 [bigo](h[super 2]).
- 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.
- 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 /h/.
- This can be seen by examination of the [@https://en.wikipedia.org/wiki/Euler-Maclaurin_formula Euler-Maclaurin summation formula],
- 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.
- Other classes of integrands which are integrated efficiently by this method are the C[sub 0][super \u221E](\u221D) [@https://en.wikipedia.org/wiki/Bump_function bump functions] and bell-shaped integrals over the infinite interval.
- For details, see [@http://epubs.siam.org/doi/pdf/10.1137/130932132 Trefethen's] SIAM review.
- In its simplest form, an integration can be performed by the following code
- using boost::math::quadrature::trapezoidal;
- auto f = [](double x) { return 1/(5 - 4*cos(x)); };
- double I = trapezoidal(f, 0.0, boost::math::constants::two_pi<double>());
- 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:
- auto bessel_integrand = [&n, &z](double theta)->std::complex<double>
- {
- std::complex<double> z{2, 3};
- using std::cos;
- using std::sin;
- return cos(z*sin(theta) - 2*theta)/pi<double>();
- };
- using boost::math::quadrature::trapezoidal;
- std::complex<double> Jnz = trapezoidal(bessel_integrand, 0.0, pi<Real>());
- 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.
- 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
- double I = trapezoidal(f, 0.0, two_pi<double>(), 1e-6);
- The routine stops when successive estimates of the integral `I1` and `I0` differ by less than the tolerance multiplied by the estimated L[sub 1] norm of the function.
- A good choice for the tolerance is [radic][epsilon], 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 [radic][epsilon], then the actual error should be ~[epsilon].
- If the integrand is *not* periodic, then reducing the error to [radic][epsilon] takes much longer, but is nonetheless possible without becoming a major performance bug.
- 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 `max_refinements` 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:
- size_t max_refinements = 15;
- long double I = trapezoidal(f, 0.0L, two_pi<long double>(), 1e-9L, max_refinements);
- Note that the maximum allowed compute time grows exponentially with `max_refinements`.
- 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:
- double error_estimate;
- double L1;
- double I = trapezoidal(f, 0.0, two_pi<double>(), tolerance, max_refinements, &error_estimate, &L1);
- if (error_estimate > tolerance*L1)
- {
- double I = some_other_quadrature_method(f, 0, two_pi<double>());
- }
- The final argument is the L[sub 1] norm of the integral.
- This is computed along with the integral, and is an essential component of the algorithm.
- First, the L[sub 1] norm establishes a scale against which the error can be measured.
- Second, the L[sub 1] norm can be used to evaluate the stability of the computation.
- This can be formulated in a rigorous manner by defining the *condition number of summation*.
- The condition number of summation is defined by
- [expression ['[kappa](S[sub n]) := [Sigma][sub i][super n] |x[sub i]|/|[Sigma][sub i][super n] x[sub i]|]]
- If this number of ~10[super k],
- then /k/ 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 \u222B |f| dt/|\u222B 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.
- As an example, we consider evaluation of Bessel functions by trapezoidal quadrature.
- The Bessel function of the first kind is defined via
- [expression ['J[sub n](x) = 1/2\u03A0 \u222B[sub -\u03A0][super \u03A0] cos(n t - x sin(t)) dt]]
- 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 /n/ the integrand oscillates rapidly,
- taking on positive and negative values,
- and hence the trapezoidal sums become ill-conditioned.
- In double precision, /x = 17/ and /n = 25/ gives a sum which is so poorly conditioned that zero correct digits are obtained.
- [optional_policy]
- References:
- Trefethen, Lloyd N., Weideman, J.A.C., ['The Exponentially Convergent Trapezoidal Rule], SIAM Review, Vol. 56, No. 3, 2014.
- Stoer, Josef, and Roland Bulirsch. ['Introduction to numerical analysis. Vol. 12.], Springer Science & Business Media, 2013.
- Higham, Nicholas J. ['Accuracy and stability of numerical algorithms.] Society for industrial and applied mathematics, 2002.
- Lyness, James N., and Cleve B. Moler. ['Numerical differentiation of analytic functions.] SIAM Journal on Numerical Analysis 4.2 (1967): 202-210.
- Gil, Amparo, Javier Segura, and Nico M. Temme. ['Computing special functions by using quadrature rules.] Numerical Algorithms 33.1-4 (2003): 265-275.
- [endsect] [/section:trapezoidal Trapezoidal Quadrature]
|