DomainIntegrals is a package designed to numerically evaluate integrals on domains as defined by the DomainSets package.

The package does not include new methods for numerical integration. It relies on other Julia packages such as QuadGK and HCubature. The methods of those packages are leveraged to evaluate integrals on more general domains than intervals and boxes.

Evaluate the integral of cos on the interval [0,1] using integral or integrate. The integral function simply returns a value, while integrate returns both the value and an estimated accuracy (as returned by the underlying packages). Integrand and domain can be specified separately or in generator form:

julia> using DomainSets, DomainIntegrals

julia> integral(cos, 0..1.0)
0.8414709848078965

julia> integral(exp(x) for x in 2..3)
12.6964808242570

julia> integral(exp(x+y) for (x,y) in (0..1)^2)
2.9524924420120535

julia> integrate(cos(x) for x in UnionDomain(0..1, 2..3))
(0.07329356604208204, 1.1102230246251565e-16)

It is possible to specify singularities of the integrand. The integration domain is split such that the singularity lies on the boundary:

julia> integral( (sin(log(abs(t))) for t in  -1..1), LogSingPoint(0.0))
-1.0000000021051316

julia> using DomainSets: ×

julia> integral( ( exp(log(abs(x-y))) for (x,y) in (2..3) × (1..4) ), SingularDiagonal())
2.333333333333333

Weighted integrals are supported through the definition of measures. A few standard weight functions are included, in particular those associated with the classical orthogonal polynomials (Legendre, Chebyshev, Jacobi, Laguerre and Hermite):

julia> integral(cos, ChebyshevTMeasure())
2.403939430634413

julia> integral(cos(t)*1/sqrt(1-t^2) for t in  -1.0..1.0)
2.403939410869398

For the particular example of the ChebyshevT measure (associated with Chebyshev polynomials of the first kind), the typical cosine map is applied which removes the algebraic endpoint singularities of the weight function, before it is evaluated numerically.

Optionally, as a first argument to integral or quadrature the user can specify a quadrature strategy. The default is AdaptiveStrategy. Explicitly providing this argument allows setting optional parameters:

julia> I, E = quadrature(QuadAdaptive(atol=1e-3, rtol = 1e-3), t->cos(t^2), 0..10)
(0.6011251848111901, 0.0004364150560137517)

A few well-known quadrature rules are included, as provided by the GaussQuadrature and FastGaussQuadrature packages. They have corresponding strategies. For example, the application of a 10-point Gauss-Laguerre rule:

julia> integral(Q_GaussLaguerre(10), cos)
0.5000005097999486

julia> integral(cos(t)*exp(-t) for t in HalfLine())
0.5

The DomainIntegrals package is extensible. The quadrature routine invokes a series of functions (integrate_property, integrate_measure, integrate_domain) that allow to dispatch on the type of singularity, measure and domain respectively. The user can add methods to these functions to teach DomainIntegrals how to evaluate new kinds of integrals. As an example of a rule that is included, the integrate_domain function dispatches on the DomainUnion type and recursively evaluates the integrals on each of the composing parts separately (if they do not overlap). The cosine map of Chebyshev measures is implemented by specializing integrate_measure for the case of a ChebyshevTMeasure. See the file in src/processing for other examples.