Template Class GKIntegrator

• Defined in File GKIntegrator.hpp

Class Documentation

template<class F, GKOrder Order>
class GKIntegrator

Adaptively integrate a vector-valued function of one variable.

An adaptation of the GSL’s gsl_integration_qag() to vector-valued integrands, using C++ templating rather than gsl_function’s pointer-to-void mechanism to pass any extra arguments f itself requires beyond the integration variable.

Template Parameters:

• F – Type of the callable integrand: either a standalone callable, invocable as f(x, …), where x is the (scalar) integration variable and … is any number of additional arguments forwarded from integrate()/quadSingle(); or a pointer to a member function, in which case the first of those additional arguments must instead be an instance of the class f is a member of (x is still passed as that member function’s own first ordinary argument, immediately after the instance) – see invokeF()’s own comment for how the two cases are told apart and reordered as needed. Either way, f must return a contiguous container of nInt doubles (e.g. std::vector<double> or std::array<double, nInt>) – the value of the integrand at x.

• Order – Which Gauss-Kronrod pair to integrate each subinterval with – see GKOrder’s own comment; selected at compile time (rather than a runtime option) so the corresponding GKRule<Order> specialization’s own fixed-size xgk/wg/wgk arrays are available with no indirection or runtime dispatch.

Public Functions

inline GKIntegrator(F f, std::size_t nInt, std::size_t maxIter = 0, double absTol = 0.0, double relTol = 1e-6)

Construct a GKIntegrator.

Parameters:

• f – The integrand; see this class’s own

• nInt – Number of quantities f returns per point

• maxIter – Maximum number of integrate()’s own bisection iterations (0 = unlimited) – see its own comment; not a count of raw evaluations of f_ itself

• absTol – Required absolute error

• relTol – Required relative error

Template Parameters:

F –

template<class ...Args>
inline auto quadSingle(const double a, const double b, Args&&... args) const -> GKIntegrationInterval

Apply a single, non-adaptive Gauss-Kronrod rule to one interval.

A direct translation of slug2’s own slug_imf_integrator<T>::integrate_gk (src/utils/ slug_imf_integrator.cpp), generalized from a single scalar-or-vector quantity T to an arbitrary, runtime-sized nInt_ of them, and with the IMF weighting integrate_gk itself folds in (imf_val) dropped entirely – f_ is evaluated on its own here, unweighted; PDFIntegrator is where that weighting belongs instead, layered on top of this class, not inside it.

Builds the symmetric grid of gknum = 2 * ngk - 1 abscissae spanning [a, b] from GKRule<Order>::xgk (see its own comment for the half-plus-center storage convention this mirrors), evaluates f_ at every one of them, then forms two weighted sums over those same nInt_-length values: the full Kronrod sum (all gknum points, weighted by GKRule<Order>::wgk) and the embedded Gauss sum (only the subset of points the lower-order Gauss rule itself uses, weighted by GKRule<Order>::wg) – including the shared center point (x = (a + b) / 2) in the Gauss sum too, whenever GKRule<Order>::ngk is even (see GKRule’s own comment for why that parity is exactly the right test). Both sums are scaled by the interval’s own half-length to convert from the standard [-1, 1] quadrature domain to [a, b], and their absolute difference is returned as this call’s own error estimate – a cheap, standard proxy for the Kronrod estimate’s own true error, since the embedded Gauss rule is enough lower-order that the two estimates should differ substantially only where the Kronrod one is itself unreliable.

Parameters:

• a – Lower limit of the interval

• b – Upper limit of the interval

• args – Any additional arguments f_ requires, forwarded to it unchanged after x at every point evaluated

Returns:

A GKIntegrationInterval with a_/b_ set to a/b, quad_ set to the ngk-point Kronrod rule’s own estimate of the integral (one value per quantity f_ returns), and err_ set to an error estimate for each of those same values – the absolute difference between the Kronrod estimate and the embedded (lower-order) Gauss rule’s own estimate of the same integral

template<class ...Args>
inline auto integrate(const double a, const double b, Args&&... args) const -> std::vector<double>

Adaptively integrate f_ over [a, b].

Starts from a single quadSingle() call over the whole [a, b] interval; if that single interval’s own largest elementwise error (max_k |err_[k]|) is already below absTol_, or its own largest elementwise relative error (max_k |err_[k]| / |quad_[k]|) is already below relTol_, returns immediately – no bisection needed.

Otherwise, repeatedly bisects whichever interval currently contributes the largest relative error to the running total (quadSum, tracked incrementally: each iteration undoes one interval’s own old contribution to quadSum/errSum, replaces that interval in place with a fresh quadSingle() over its own left half, appends a new interval for its own right half, and adds both halves’ own contributions back in), stopping – and returning quadSum – as soon as either the largest element of errSum drops below absTol_, or the largest elementwise relative error (|errSum[k]| / |quadSum[k]|) drops below relTol_, or the number of bisection iterations exceeds maxIter_ (0 = unlimited) – not a count of raw evaluations of f_ itself, which each iteration costs two further quadSingle() calls’ worth of (2 * (2 * ngk - 1) further evaluations).

Which interval to bisect next is decided by each interval’s own relative error against the global running total (|intervals[i].err_[k]| / |quadSum[k]|, maximized over k), not against that interval’s own local quad_ – so bisection effort concentrates on whichever interval’s own uncertainty matters most to the overall answer, mirroring slug2’s own identical convention in slug_imf_integrator<T>:: integrate_range (see quadSingle()’s own comment for that function’s sibling, integrate_gk, which this class’s quadSingle() itself translates).

Parameters:

• a – Lower limit of integration

• b – Upper limit of integration

• args – Any additional arguments f_ requires; see quadSingle()’s own args parameter

Returns:

The integral estimate, one value per quantity f_ returns (length nInt_)