IntegrationToolbox.jl

A Julia package to facilitate numerical integration in one or several dimensions, including the flexible use of finite and infinite limits across different dimensions by a change of variables, and the relative ease of switching between different integration routines.

The main functionality is provided by the callable objects Integration (for 1D integration) and MultivariateIntegration (for multidimensional integration). Additionally, the types IntegralTransformation and MultivariateIntegralTransformation allow for specifying a change of variables. In particular, this is useful to treat infinite and semi-infinite limits, as well as changing the region of integration to the standard interval of a fixed point Gaussian quadrature algorithm based on orthogonal polynomials. See also the QuadGK.jl documentation for further details. Several often needed transformations are provided, but it is also possible to write custom transformations.

To integrate a function f with the routine routine (which might take additional arguments args and keyword arguments kwargs), the general usage is as follows

using IntegrationToolbox
integrate_function = Integration(integrand=f, routine=routine, args=args, kwargs=kwargs)
result = integrate_function()

The minimum requirement for the integration routine method is that it takes the integrand function as its first argument, i.e. it can be called as routine(f), but may also take additional arguments and/or keyword arguments routine(f, args...; kwargs...).

Examples

To integrate the function f(x) =x^2 over the interval [0,1] with QuadGK.jl, we can write

using IntegrationToolbox, QuadGK
f(x) = x^2
integrate_f = Integration(integrand=f, routine=first ∘ quadgk, args=(0, 1))
result = integrate_f()

where the integration limits are given as mandatory positional arguments to quadgk.

For integrating the function f(x, c) = c x^2 parametrized by the value c=2 over the interval [-1,1] with a 5-point Gauss-Legendre rule, we first define the function

function quad(f, rule, N)
    x, w = rule(N)
    return dot(w, f.(x))
end

written as to accept f as its first argument. It also requires a Gaussian quadrature rule to calculate the corresponding nodes and weights, here provided by FastGaussQuadrature,jl, as well as the order N of the quadrature rule. For convenience, it is also exported by MultivariateQuadrature.jl We can then use it as follows

using IntegrationToolbox, FastGaussQuadrature
f(x, c) = c * x^2
integrate_f = Integration(integrand=f, routine=quad, args=(gausslegendre, 5))
result = integrate_f(2)

No explicit specification of limits was necessary, as we are integrating over the default interval [-1,1] for Gauss-Legendre quadrature.

As an example of a multidimensional integration, we integrate the function h(x_1, x_2, c) = c x_1^3 * exp(-x_2^2) with c=3 over the interval [2,3] x [-∞,∞] with the Monte-Carlo integrator Vegas from Cuba.jl. The latter requires the input to integrated over the square region [0, 1] x [0, 1], so we need to introduce the different integration limits via a change of variables. The call now looks like

using IntegrationToolbox, Cuba 
transformation_x2 = ABToZeroOne(-1,1) ∘ MinusInfInfToMinusOneOne()
transformation = MultivariateIntegralTransformation(ABToZeroOne(2,3), transformation_x2)
h(x, c) = c * x[1]^3 * exp(-x[2]^2)
integrate_h = MultivariateIntegration(integrand=h, dimension=2, transformation=transformation, routine= first ∘ first ∘ cubature_vegas, args=((2,)))
result = integrate_h(3)

This also demonstrates the use of a chained transform that combines several changes of variables into a single one via the composition operator ∘.

Finally, we use MultivariateQuadrature.jl to perform a two-dimensional fixed-point Gaussian quadrature to integrate g(x_1, x_2) = x_1 * exp(-x_2) over the interval [0,1] x [0,∞], using a Gauss-Laguerre and a Gauss-Legendre rule.

using IntegrationToolbox, MultivariateQuadrature 
transformation = MultivariateIntegralTransformation(ABToMinusOneOne(0, 1), IdentityTransformation())
g(x) = x[1]
integrate_g = MultivariateIntegration(integrand=g, dimension=2, routine=mvquad, transformation=transformation, args=((gausslegendre, gausslaguerre), (30, 30)))
result = integrate_g()
end

Since the exponential decay is implicitly contained in the Gauss-Laguerre rule, we use an identity transformation for the second argument that leaves everything unchanged.

Further examples can be found in the testfolder.

Related packages

Integrals.jl is a meta integration package providing a standardized interface for the flexible use of different integration routines fully compatible with the SciML universe. Parallelization is also supported.