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, \u221E).
+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, \u221E) and(-\u221E, b).
+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.
[[sinh_sinh] [(-[infin],[infin])] [][]]
]
-[endsect]
+[endsect] [/section:de_overview Overview]
[section:de_tanh_sinh tanh_sinh]
};
-The tanh-sinh quadrature routine provided by boost is a rapidly convergent numerical integration scheme for holomorphic integrands.
+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,
+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:
+
+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 (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.
+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/),
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:
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());
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.
+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:
+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]
return integrator.integrate(f, 0, boost::math::constants::half_pi<value_type>()) + boost::math::constants::half_pi<value_type>();
}
-[endsect]
+[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]
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.
+Not only is this form accurate to full machine-precision, but it converges to the result faster as well.
-[endsect]
+[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]
double L1;
double Q = integrator.integrate(f, &error, &L1);
-Note that the limits of integration are understood to be (-\u221E, \u221E).
+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:
}
-[endsect]
+[endsect] [/section:de_sinh_sinh sinh_sinh]
[section:de_exp_sinh exp_sinh]
size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
};
-For half-infinite intervals, the exp-sinh quadrature is provided:
+For half-infinite intervals, the `exp-sinh` quadrature is provided:
exp_sinh<double> integrator;
auto f = [](double x) { return exp(-3*x); };
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.
+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.
+Endpoint singularities and complex-valued integrands are supported by `exp-sinh`.
For example, the modified Bessel function K can be represented via:
`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]
+[endsect] [/section:de_exp_sinh exp_sinh]
[section:de_tol Setting the Termination Condition for Integration]
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,
+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]
+[endsect] [/section:de_tol Setting the Termination Condition for Integration]
[section:de_levels Setting the Maximum Interval Halvings and Memory Requirements]
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]
+[endsect] [/section:de_levels Setting the Maximum Interval Halvings and Memory Requirements]
[section:de_thread Thread Safety]
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]
+[endsect] [/section:de_thread Thread Safety]
[section:de_caveats Caveats]
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:
- 1 / (1 +x^2)
+[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
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
+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;
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
+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.
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.
+This form integrates just fine over (-log([pi]/2), +[infin]) using either the `tanh_sinh` or `exp_sinh` classes.
-
-[endsect]
+[endsect] [/section:de_caveats Caveats]
[section:de_refes References]
* 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]
-[endsect]
+[endsect] [/section:de_refes References]
+
+[endsect] [/section:double_exponential Double-exponential quadrature]
projects / platform / upstream / boost.git / blobdiff