551 lines
29 KiB
Plaintext
551 lines
29 KiB
Plaintext
[/
|
||
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:double_exponential Double-exponential quadrature]
|
||
|
||
[section:de_overview Overview]
|
||
|
||
[heading Synopsis]
|
||
|
||
``
|
||
#include <boost/math/quadrature/tanh_sinh.hpp>
|
||
#include <boost/math/quadrature/exp_sinh.hpp>
|
||
#include <boost/math/quadrature/sinh_sinh.hpp>
|
||
|
||
namespace boost{ namespace math{
|
||
|
||
template<class Real>
|
||
class tanh_sinh
|
||
{
|
||
public:
|
||
tanh_sinh(size_t max_refinements = 15, const Real& min_complement = tools::min_value<Real>() * 4)
|
||
|
||
template<class F>
|
||
auto integrate(const F f, Real a, Real b,
|
||
Real tolerance = tools::root_epsilon<Real>(),
|
||
Real* error = nullptr,
|
||
Real* L1 = nullptr,
|
||
std::size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
|
||
|
||
template<class F>
|
||
auto integrate(const F f, Real
|
||
tolerance = tools::root_epsilon<Real>(),
|
||
Real* error = nullptr,
|
||
Real* L1 = nullptr,
|
||
std::size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
|
||
|
||
};
|
||
|
||
template<class Real>
|
||
class exp_sinh
|
||
{
|
||
public:
|
||
exp_sinh(size_t max_refinements = 9);
|
||
|
||
template<class F>
|
||
auto integrate(const F f, Real a, Real b,
|
||
Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
|
||
Real* error = nullptr,
|
||
Real* L1 = nullptr,
|
||
size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
|
||
template<class F>
|
||
auto integrate(const F f,
|
||
Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
|
||
Real* error = nullptr,
|
||
Real* L1 = nullptr,
|
||
size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
|
||
};
|
||
|
||
template<class Real>
|
||
class sinh_sinh
|
||
{
|
||
public:
|
||
sinh_sinh(size_t max_refinements = 9);
|
||
|
||
template<class F>
|
||
auto integrate(const F f,
|
||
Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
|
||
Real* error = nullptr,
|
||
Real* L1 = nullptr,
|
||
size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
|
||
};
|
||
|
||
}}
|
||
``
|
||
|
||
These three integration routines provide robust general purpose quadrature, each having a "native" range over which
|
||
quadrature is performed.
|
||
For example, the `sinh_sinh` quadrature integrates over the entire real line, the `tanh_sinh` over (-1, 1),
|
||
and the `exp_sinh` over (0, [infin]).
|
||
The latter integrators also have auxilliary ranges which are handled via a change of variables on the function being integrated,
|
||
so that the `tanh_sinh` can handle integration over /(a, b)/, and `exp_sinh` over /(a, [infin]) and(-[infin], b)/.
|
||
|
||
Like the other quadrature routines in Boost, these routines support both real and complex-valued integrands.
|
||
|
||
The `integrate` methods which do not specify a range always integrate over the native range of the method, and generally
|
||
are the most efficient and produce the smallest code, on the other hand the methods which do specify the bounds of integration
|
||
are the most general, and use argument transformations which are generally very robust. The following table summarizes
|
||
the ranges supported by each method:
|
||
|
||
[table
|
||
[[Integrator][Native range][Other supported ranges][Comments]]
|
||
[[tanh_sinh] [(-1,1)] [(a,b)[br](a,[infin])[br](-[infin],b)[br](-[infin],[infin])]
|
||
[Special care is taken for endpoints at or near zero to ensure that abscissa values are calculated without the loss of precision
|
||
that would normally occur. Likewise when transforming to an infinite endpoint, the additional information which tanh_sinh has
|
||
internally on abscissa values is used to ensure no loss of precision during the transformation.]]
|
||
[[exp_sinh] [(0,[infin])] [(a,[infin])[br](-[infin],0)[br](-[infin],b)] []]
|
||
[[sinh_sinh] [(-[infin],[infin])] [][]]
|
||
]
|
||
|
||
[endsect] [/section:de_overview Overview]
|
||
|
||
[section:de_tanh_sinh tanh_sinh]
|
||
|
||
template<class Real>
|
||
class tanh_sinh
|
||
{
|
||
public:
|
||
tanh_sinh(size_t max_refinements = 15, const Real& min_complement = tools::min_value<Real>() * 4)
|
||
|
||
template<class F>
|
||
auto integrate(const F f, Real a, Real b,
|
||
Real tolerance = tools::root_epsilon<Real>(),
|
||
Real* error = nullptr,
|
||
Real* L1 = nullptr,
|
||
std::size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
|
||
|
||
template<class F>
|
||
auto integrate(const F f, Real
|
||
tolerance = tools::root_epsilon<Real>(),
|
||
Real* error = nullptr,
|
||
Real* L1 = nullptr,
|
||
std::size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
|
||
|
||
};
|
||
|
||
The `tanh-sinh` quadrature routine provided by boost is a rapidly convergent numerical integration scheme for holomorphic integrands.
|
||
By this we mean that the integrand is the restriction to the real line of a complex-differentiable function which is bounded on the interior of the unit disk /|z| < 1/,
|
||
so that it lies within the so-called [@https://en.wikipedia.org/wiki/Hardy_space Hardy space].
|
||
If your integrand obeys these conditions, it can be shown that `tanh-sinh` integration is optimal,
|
||
in the sense that it requires the fewest function evaluations for a given accuracy of any quadrature algorithm for a random element from the Hardy space.
|
||
|
||
A basic example of how to use the `tanh-sinh` quadrature is shown below:
|
||
|
||
tanh_sinh<double> integrator;
|
||
auto f = [](double x) { return 5*x + 7; };
|
||
// Integrate over native bounds of (-1,1):
|
||
double Q = integrator.integrate(f);
|
||
// Integrate over (0,1.1) instead:
|
||
Q = integrator.integrate(f, 0.0, 1.1);
|
||
|
||
The basic idea of `tanh-sinh` quadrature is that a variable transformation can cause the endpoint derivatives to decay rapidly.
|
||
When the derivatives at the endpoints decay much faster than the Bernoulli numbers grow,
|
||
the Euler-Maclaurin summation formula tells us that simple trapezoidal quadrature converges faster than any power of /h/.
|
||
That means the number of correct digits of the result should roughly double with each new level of integration (halving of /h/),
|
||
Hence the default termination condition for integration is usually set to the square root of machine epsilon.
|
||
Most well-behaved integrals should converge to full machine precision with this termination condition,
|
||
and in 6 or fewer levels at double precision, or 7 or fewer levels for quad precision.
|
||
|
||
One very nice property of tanh-sinh quadrature is that it can handle singularities at the endpoints of the integration domain.
|
||
For instance, the following integrand, singular at both endpoints, can be efficiently evaluated to 100 binary digits:
|
||
|
||
auto f = [](Real x) { return log(x)*log1p(-x); };
|
||
Real Q = integrator.integrate(f, (Real) 0, (Real) 1);
|
||
|
||
Now onto the caveats: As stated before, the integrands must lie in a Hardy space to ensure rapid convergence.
|
||
Attempting to integrate a function which is not bounded on the unit disk by tanh-sinh can lead to very slow convergence.
|
||
For example, take the Runge function:
|
||
|
||
auto f1 = [](double t) { return 1/(1+25*t*t); };
|
||
Q = integrator.integrate(f1, (double) -1, (double) 1);
|
||
|
||
This function has poles at \u00B1 \u2148/5, and as such it is not bounded on the unit disk.
|
||
However, the related function
|
||
|
||
auto f2 = [](double t) { return 1/(1+0.04*t*t); };
|
||
Q = integrator.integrate(f2, (double) -1, (double) 1);
|
||
|
||
has poles outside the unit disk (at \u00B1 5\u2148), and is therefore in the Hardy space.
|
||
Our benchmarks show that the second integration is performed 22x faster than the first!
|
||
If you do not understand the structure of your integrand in the complex plane, do performance testing before deployment.
|
||
|
||
Like the trapezoidal quadrature, the tanh-sinh quadrature produces an estimate of the L[sub 1] norm of the integral along with the requested integral.
|
||
This is to establish a scale against which to measure the tolerance, and to provide an estimate of the condition number of the summation.
|
||
This can be queried as follows:
|
||
|
||
tanh_sinh<double> integrator;
|
||
auto f = [](double x) { return 5*x + 7; };
|
||
double termination = std::sqrt(std::numeric_limits<double>::epsilon());
|
||
double error;
|
||
double L1;
|
||
size_t levels;
|
||
double Q = integrator.integrate(f, 0.0, 1.0, termination, &error, &L1, &levels);
|
||
double condition_number = L1/std::abs(Q);
|
||
|
||
If the condition number is large, the computed integral is worthless: typically one can assume that Q has lost one digit of precision
|
||
when the condition number of O(10^Q).
|
||
The returned error term is not the actual error in the result, but merely an a posteriori error estimate.
|
||
It is the absolute difference between the last two approximations, and for well behaved integrals, the actual error should be very much smaller than this.
|
||
The following table illustrates how the errors and conditioning vary for few sample integrals, in each case the termination condition was set
|
||
to the square root of epsilon, and all tests were conducted in double precision:
|
||
|
||
[table
|
||
[[Integral][Range][Error][Actual measured error][Levels][Condition Number][Comments]]
|
||
[[`5 * x + 7`][(0,1)][3.5e-15][0][5][1][This trivial case shows just how accurate these methods can be.]]
|
||
[[`log(x) * log(x)`][0, 1)][0][0][5][1][This is an example of an integral that Gaussian integrators fail to handle.]]
|
||
[[`exp(-x) / sqrt(x)`][(0,+[infin])][8.0e-10][1.1e-15][5][1][Gaussian integrators typically fail to handle the singularities at the endpoints of this one.]]
|
||
[[`x * sin(2 * exp(2 * sin(2 * exp(2 * x))))`][(-1,1)][7.2e-16][4.9e-17][9][1.89][This is a truely horrible integral that oscillates wildly and
|
||
unpredictably with some very sharp "spikes" in it's graph. The higher number of levels used reflects the difficulty of sampling the more extreme features.]]
|
||
[[`x == 0 ? 1 : sin(x) / x`][(-[infin], [infin])][3.0e-1][4.0e-1][15][159][This highly oscillatory integral isn't handled at all well by tanh-sinh quadrature: there is so much
|
||
cancellation in the sum that the result is essentially worthless. The argument transformation of the infinite integral behaves somewhat badly as well, in fact
|
||
we do ['slightly] better integrating over 2 symmetrical and large finite limits.]]
|
||
[[`sqrt(x / (1 - x * x))`][(0,1)][1e-8][1e-8][5][1][This an example of an integral that has all its area close to a non-zero endpoint, the problem here is that
|
||
the function being integrated returns "garbage" values for x very close to 1. We can easily fix this issue by passing a 2 argument functor to the integrator:
|
||
the second argument gives the distance to the nearest endpoint, and we can use that information to return accurate values, and thus fix the integral calculation.]]
|
||
[[`x < 0.5 ? sqrt(x) / sqrt(1 - x * x) : sqrt(x / ((x + 1) * (xc)))`][(0,1)][0][0][5][1][This is the 2-argument version of the previous integral, the second
|
||
argument ['xc] is `1-x` in this case, and we use 1-x[super 2] == (1-x)(1+x) to calculate 1-x[super 2] with greater accuracy.]]
|
||
]
|
||
|
||
Although the `tanh-sinh` quadrature can compute integral over infinite domains by variable transformations, these transformations can create a very poorly behaved integrand.
|
||
For this reason, double-exponential variable transformations have been provided that allow stable computation over infinite domains; these being the exp-sinh and sinh-sinh quadrature.
|
||
|
||
[h4 Complex integrals]
|
||
|
||
The `tanh_sinh` integrator supports integration of functions which return complex results, for example the sine-integral `Si(z)` has the integral representation:
|
||
|
||
[equation sine_integral]
|
||
|
||
Which we can code up directly as:
|
||
|
||
template <class Complex>
|
||
Complex Si(Complex z)
|
||
{
|
||
typedef typename Complex::value_type value_type;
|
||
using std::sin; using std::cos; using std::exp;
|
||
auto f = [&z](value_type t) { return -exp(-z * cos(t)) * cos(z * sin(t)); };
|
||
boost::math::quadrature::tanh_sinh<value_type> integrator;
|
||
return integrator.integrate(f, 0, boost::math::constants::half_pi<value_type>()) + boost::math::constants::half_pi<value_type>();
|
||
}
|
||
|
||
[endsect] [/section:de_tanh_sinh tanh_sinh]
|
||
|
||
[section:de_tanh_sinh_2_arg Handling functions with large features near an endpoint with tanh-sinh quadrature]
|
||
|
||
Tanh-sinh quadrature has a unique feature which makes it well suited to handling integrals with either singularities or large features of interest
|
||
near one or both endpoints, it turns out that when we calculate and store the abscissa values at which we will be evaluating the function, we can equally
|
||
well calculate the difference between the abscissa value and the nearest endpoint.
|
||
This makes it possible to perform quadrature arbitrarily close to an endpoint, without suffering cancellation error.
|
||
Note however, that we never actually reach the endpoint, so any endpoint singularity will always be excluded from the quadrature.
|
||
|
||
The tanh_sinh integration routine will use this additional information internally when performing range transformation, so that for example,
|
||
if one end of the range is zero (or infinite) then our transformations will get arbitrarily close to the endpoint without precision loss.
|
||
|
||
However, there are some integrals which may have all of their area near ['both] endpoints, or else near the non-zero endpoint, and in that situation
|
||
there is a very real risk of loss of precision. For example:
|
||
|
||
tanh_sinh<double> integrator;
|
||
auto f = [](double x) { return sqrt(x / (1 - x * x); };
|
||
double Q = integrator.integrate(f, 0.0, 1.0);
|
||
|
||
Results in very low accuracy as all the area of the integral is near 1, and the `1 - x * x` term suffers from cancellation error here.
|
||
|
||
However, both of tanh_sinh's integration routines will automatically handle 2 argument functors: in this case the first argument is the abscissa value as
|
||
before, while the second is the distance to the nearest endpoint, ie `a - x` or `b - x` if we are integrating over (a,b).
|
||
You can always differentiate between these 2 cases because the second argument will be negative if we are nearer to the left endpoint.
|
||
|
||
Knowing this, we can rewrite our lambda expression to take advantage of this additional information:
|
||
|
||
tanh_sinh<double> integrator;
|
||
auto f = [](double x, double xc) { return x <= 0.5 ? sqrt(x) / sqrt(1 - x * x) : sqrt(x / ((x + 1) * (xc))); };
|
||
double Q = integrator.integrate(f, 0.0, 1.0);
|
||
|
||
Not only is this form accurate to full machine-precision, but it converges to the result faster as well.
|
||
|
||
[endsect] [/section:de_tanh_sinh_2_arg Handling functions with large features near an endpoint with tanh-sinh quadrature]
|
||
|
||
[section:de_sinh_sinh sinh_sinh]
|
||
|
||
template<class Real>
|
||
class sinh_sinh
|
||
{
|
||
public:
|
||
sinh_sinh(size_t max_refinements = 9);
|
||
|
||
template<class F>
|
||
auto integrate(const F f,
|
||
Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
|
||
Real* error = nullptr,
|
||
Real* L1 = nullptr,
|
||
size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;;
|
||
};
|
||
|
||
The sinh-sinh quadrature allows computation over the entire real line, and is called as follows:
|
||
|
||
sinh_sinh<double> integrator;
|
||
auto f = [](double x) { return exp(-x*x); };
|
||
double error;
|
||
double L1;
|
||
double Q = integrator.integrate(f, &error, &L1);
|
||
|
||
Note that the limits of integration are understood to be (-[infin], +[infin]).
|
||
|
||
Complex valued integrands are supported as well, for example the [@https://en.wikipedia.org/wiki/Dirichlet_eta_function Dirichlet Eta function]
|
||
can be represented via:
|
||
|
||
[equation complex_eta_integral]
|
||
|
||
which we can directly code up as:
|
||
|
||
template <class Complex>
|
||
Complex eta(Complex s)
|
||
{
|
||
typedef typename Complex::value_type value_type;
|
||
using std::pow; using std::exp;
|
||
Complex i(0, 1);
|
||
value_type pi = boost::math::constants::pi<value_type>();
|
||
auto f = [&, s, i](value_type t) { return pow(0.5 + i * t, -s) / (exp(pi * t) + exp(-pi * t)); };
|
||
boost::math::quadrature::sinh_sinh<value_type> integrator;
|
||
return integrator.integrate(f);
|
||
}
|
||
|
||
|
||
[endsect] [/section:de_sinh_sinh sinh_sinh]
|
||
|
||
[section:de_exp_sinh exp_sinh]
|
||
|
||
template<class Real>
|
||
class exp_sinh
|
||
{
|
||
public:
|
||
exp_sinh(size_t max_refinements = 9);
|
||
|
||
template<class F>
|
||
auto integrate(const F f, Real a, Real b,
|
||
Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
|
||
Real* error = nullptr,
|
||
Real* L1 = nullptr,
|
||
size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
|
||
template<class F>
|
||
auto integrate(const F f,
|
||
Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
|
||
Real* error = nullptr,
|
||
Real* L1 = nullptr,
|
||
size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
|
||
};
|
||
|
||
For half-infinite intervals, the `exp-sinh` quadrature is provided:
|
||
|
||
exp_sinh<double> integrator;
|
||
auto f = [](double x) { return exp(-3*x); };
|
||
double termination = sqrt(std::numeric_limits<double>::epsilon());
|
||
double error;
|
||
double L1;
|
||
double Q = integrator.integrate(f, termination, &error, &L1);
|
||
|
||
|
||
The native integration range of this integrator is (0, [infin]), but we also support /(a, [infin]), (-[infin], 0)/ and /(-[infin], b)/ via argument transformations.
|
||
|
||
Endpoint singularities and complex-valued integrands are supported by `exp-sinh`.
|
||
|
||
For example, the modified Bessel function K can be represented via:
|
||
|
||
[equation complex_bessel_k_integral]
|
||
|
||
Which we can code up as:
|
||
|
||
template <class Complex>
|
||
Complex bessel_K(Complex alpha, Complex z)
|
||
{
|
||
typedef typename Complex::value_type value_type;
|
||
using std::cosh; using std::exp;
|
||
auto f = [&, alpha, z](value_type t)
|
||
{
|
||
value_type ct = cosh(t);
|
||
if (ct > log(std::numeric_limits<value_type>::max()))
|
||
return Complex(0);
|
||
return exp(-z * ct) * cosh(alpha * t);
|
||
};
|
||
boost::math::quadrature::exp_sinh<value_type> integrator;
|
||
return integrator.integrate(f);
|
||
}
|
||
|
||
The only wrinkle in the above code is the need to check for large `cosh(t)` in which case we assume that
|
||
`exp(-x cosh(t))` tends to zero faster than `cosh(alpha x)` tends to infinity and return `0`. Without that
|
||
check we end up with `0 * Infinity` as the result (a NaN).
|
||
|
||
[endsect] [/section:de_exp_sinh exp_sinh]
|
||
|
||
[section:de_tol Setting the Termination Condition for Integration]
|
||
|
||
The integrate method for all three double-exponential quadratures supports ['tolerance] argument that acts as the
|
||
termination condition for integration.
|
||
|
||
The tolerance is met when two subsequent estimates of the integral have absolute error less than `tolerance*L1`.
|
||
|
||
It is highly recommended that the tolerance be left at the default value of [radic][epsilon], or something similar.
|
||
Since double exponential quadrature converges exponentially fast for functions in Hardy spaces, then once the routine has *proved* that the error is ~[radic][epsilon],
|
||
then the error should in fact be ~[epsilon].
|
||
|
||
If you request that the error be ~[epsilon], this tolerance might never be achieved (as the summation is not stabilized ala Kahan),
|
||
and the routine will simply flounder,
|
||
dividing the interval in half in order to increase the precision of the integrand, only to be thwarted by floating point roundoff.
|
||
|
||
If for some reason, the default value doesn't quite achieve full precision, then you could try something a little smaller such as
|
||
[radic][epsilon]/4 or [epsilon][super 2/3].
|
||
However, more likely, you need to check that your function to be integrated is able to return accurate values, and that there are no other issues with your integration scheme.
|
||
|
||
[endsect] [/section:de_tol Setting the Termination Condition for Integration]
|
||
|
||
[section:de_levels Setting the Maximum Interval Halvings and Memory Requirements]
|
||
|
||
The max interval halvings is the maximum number of times the interval can be cut in half before giving up.
|
||
If the accuracy is not met at that time, the routine returns its best estimate, along with the `error` and `L1`,
|
||
which allows the user to decide if another quadrature routine should be employed.
|
||
|
||
An example of this is
|
||
|
||
double tol = std::sqrt(std::numeric_limits<double>::epsilon());
|
||
size_t max_halvings = 12;
|
||
tanh_sinh<double> integrator(max_halvings);
|
||
auto f = [](double x) { return 5*x + 7; };
|
||
double error, L1;
|
||
double Q = integrator.integrate(f, (double) 0, (double) 1, &error, &L1);
|
||
if (error*L1 > 0.01)
|
||
{
|
||
Q = some_other_quadrature_method(f, (double) 0, (double) 1);
|
||
}
|
||
|
||
It's important to remember that the number of sample points doubles with each new level, as does the memory footprint
|
||
of the integrator object. Further, if the integral is smooth, then the precision will be doubling with each new level,
|
||
so that for example, many integrals can achieve 100 decimal digit precision after just 7 levels. That said, abscissa-weight
|
||
pairs for new levels are computed only when a new level is actually required (see thread safety), none the less,
|
||
you should avoid setting the maximum arbitrarily high "just in case" as the time and space requirements for a large
|
||
number of levels can quickly grow out of control.
|
||
|
||
[endsect] [/section:de_levels Setting the Maximum Interval Halvings and Memory Requirements]
|
||
|
||
[section:de_thread Thread Safety]
|
||
|
||
All three of the double-exponential integrators are thread safe as long as BOOST_MATH_NO_ATOMIC_INT is not set. Since the
|
||
integrators store a large amount of fairly hard to compute data, it is recommended that these objects are stored and reused
|
||
as much as possible.
|
||
|
||
Internally all three of the double-exponential integrators use the same caching strategy: they allocate all the vectors needed
|
||
to store the maximum permitted levels, but only populate the first few levels when constructed. This means a minimal amount of memory
|
||
is actually allocated when the integrator is first constructed, and already populated levels can be accessed via a lockfree
|
||
atomic read, and only populating new levels requires a thread lock.
|
||
|
||
In addition, the three built in types (plus `__float128` when available), have the first 7 levels pre-computed: this is generally sufficient for the vast majority
|
||
of integrals - even at quad precision - and means that integrators for these types are relatively cheap to construct.
|
||
|
||
[endsect] [/section:de_thread Thread Safety]
|
||
|
||
[section:de_caveats Caveats]
|
||
|
||
A few things to keep in mind while using the tanh-sinh, exp-sinh, and sinh-sinh quadratures:
|
||
|
||
These routines are *very* aggressive about approaching the endpoint singularities.
|
||
This allows lots of significant digits to be extracted, but also has another problem: Roundoff error can cause the function to be evaluated at the endpoints.
|
||
A few ways to avoid this: Narrow up the bounds of integration to say, [a + [epsilon], b - [epsilon]], make sure (a+b)/2 and (b-a)/2 are representable, and finally,
|
||
if you think the compromise between accuracy an usability has gone too far in the direction of accuracy, file a ticket.
|
||
|
||
Both exp-sinh and sinh-sinh quadratures evaluate the functions they are passed at *very* large argument.
|
||
You might understand that x[super 12]exp(-x) is should be zero when x[super 12] overflows, but IEEE floating point arithmetic does not.
|
||
Hence `std::pow(x, 12)*std::exp(-x)` is an indeterminate form whenever `std::pow(x, 12)` overflows.
|
||
So make sure your functions have the correct limiting behavior; for example
|
||
|
||
auto f = [](double x) {
|
||
double t = exp(-x);
|
||
if(t == 0)
|
||
{
|
||
return 0;
|
||
}
|
||
return t*pow(x, 12);
|
||
};
|
||
|
||
has the correct behavior for large /x/, but `auto f = [](double x) { return exp(-x)*pow(x, 12); };` does not.
|
||
|
||
Oscillatory integrals, such as the sinc integral, are poorly approximated by double-exponential quadrature.
|
||
Fortunately the error estimates and L1 norm are massive for these integrals, but nonetheless, oscillatory integrals require different techniques.
|
||
|
||
A special mention should be made about integrating through zero: while our range adaptors preserve precision when one endpoint is zero,
|
||
things get harder when the origin is neither in the center of the range, nor at an endpoint. Consider integrating:
|
||
|
||
[expression 1 / (1 +x^2)]
|
||
|
||
Over (a, [infin]). As long as `a >= 0` both the tanh_sinh and the exp_sinh integrators will handle this just fine: in fact they provide
|
||
a rather efficient method for this kind of integral. However, if we have `a < 0` then we are forced to adapt the range in a way that
|
||
produces abscissa values near zero that have an absolute error of [epsilon], and since all of the area of the integral is near zero
|
||
both integrators thrash around trying to reach the target accuracy, but never actually get there for `a << 0`. On the other hand, the
|
||
simple expedient of breaking the integral into two domains: (a, 0) and (0, b) and integrating each seperately using the tanh-sinh
|
||
integrator, works just fine.
|
||
|
||
Finally, some endpoint singularities are too strong to be handled by `tanh_sinh` or equivalent methods, for example consider integrating
|
||
the function:
|
||
|
||
double p = some_value;
|
||
tanh_sinh<double> integrator;
|
||
auto f = [&](double x){ return pow(tan(x), p); };
|
||
auto Q = integrator.integrate(f, 0, constants::half_pi<double>());
|
||
|
||
The first problem with this function, is that the singularity is at [pi]/2, so if we're integrating over (0, [pi]/2) then we can never
|
||
approach closer to the singularity than [epsilon], and for p less than but close to 1, we need to get ['very] close to the singularity
|
||
to find all the area under the function. If we recall the identity [^tan([pi]/2 - x) == 1/tan(x)] then we can rewrite the function like this:
|
||
|
||
auto f = [&](double x){ return pow(tan(x), -p); };
|
||
|
||
And now the singularity is at the origin and we can get much closer to it when evaluating the integral: all we have done is swap the
|
||
integral endpoints over.
|
||
|
||
This actually works just fine for p < 0.95, but after that the `tanh_sinh` integrator starts thrashing around and is unable to
|
||
converge on the integral. The problem is actually a lack of exponent range: if we simply swap type double for something
|
||
with a greater exponent range (an 80-bit long double or a quad precision type), then we can get to at least p = 0.99. If we want to go
|
||
beyond that, or stick with type double, then we have to get smart.
|
||
|
||
The easiest method is to notice that for small x, then [^tan(x) [cong] x], and so we are simply integrating x[super -p]. Therefore we can use
|
||
this approximation over (0, small), and integrate numerically from (small, [pi]/2), using [epsilon] as a suitable crossover point
|
||
seems sensible:
|
||
|
||
double p = some_value;
|
||
double crossover = std::numeric_limits<double>::epsilon();
|
||
tanh_sinh<double> integrator;
|
||
auto f = [&](double x){ return pow(tan(x), p); };
|
||
auto Q = integrator.integrate(f, crossover, constants::half_pi<double>()) + pow(crossover, 1 - p) / (1 - p);
|
||
|
||
There is an alternative, more complex method, which is applicable when we are dealing with expressions which can be simplified
|
||
by evaluating by logs. Let's suppose that as in this case, all the area under the graph is infinitely close to zero, now inagine
|
||
that we could expand that region out over a much larger range of abscissa values: that's exactly what happens if we perform
|
||
argument substitution, replacing `x` by `exp(-x)` (note that we must also multiply by the derivative of `exp(-x)`).
|
||
Now the singularity at zero is moved to +[infin], and the [pi]/2 bound to
|
||
-log([pi]/2). Initially our argument substituted function looks like:
|
||
|
||
auto f = [&](double x){ return exp(-x) * pow(tan(exp(-x)), -p); };
|
||
|
||
Which is hardly any better, as we still run out of exponent range just as before. However, if we replace `tan(exp(-x))` by `exp(-x)` for suitably
|
||
small `exp(-x)`, and therefore [^x > -log([epsilon])], we can greatly simplify the expression and evaluate by logs:
|
||
|
||
auto f = [&](double x)
|
||
{
|
||
static const double crossover = -log(std::numeric_limits<double>::epsilon());
|
||
return x > crossover ? exp((p - 1) * x) : exp(-x) * pow(tan(exp(-x)), -p);
|
||
};
|
||
|
||
This form integrates just fine over (-log([pi]/2), +[infin]) using either the `tanh_sinh` or `exp_sinh` classes.
|
||
|
||
[endsect] [/section:de_caveats Caveats]
|
||
|
||
[section:de_refes References]
|
||
|
||
* Hidetosi Takahasi and Masatake Mori, ['Double Exponential Formulas for Numerical Integration] Publ. Res. Inst. Math. Sci., 9 (1974), pp. 721-741.
|
||
* Masatake Mori, ['An IMT-Type Double Exponential Formula for Numerical Integration], Publ RIMS, Kyoto Univ. 14 (1978), 713-729.
|
||
* David H. Bailey, Karthik Jeyabalan and Xiaoye S. Li ['A Comparison of Three High-Precision Quadrature Schemes] Office of Scientific & Technical Information Technical Reports.
|
||
* Tanaka, Ken’ichiro, et al. ['Function classes for double exponential integration formulas.] Numerische Mathematik 111.4 (2009): 631-655.
|
||
|
||
[endsect] [/section:de_refes References]
|
||
|
||
[endsect] [/section:double_exponential Double-exponential quadrature]
|