Numerical Integration

AppleAccelerate wraps Apple's Quadrature library (a C port of QUADPACK) for adaptive numerical integration of real scalar functions.

Namespace

integrate is not exported. Access it via the AppleAccelerate. prefix.

`integrate`

AppleAccelerate.integrate(f, a, b; integrator=:qags, abstol=1e-8, reltol=1e-8,
                          max_intervals=200, qag_points=0)

Integrate f over the interval (a, b). f is an ordinary Julia function f(x::Float64) -> Float64; Accelerate evaluates it in batches, which is handled internally. The call returns a named tuple (value, abserr, status) where value is the integral estimate, abserr is the estimated absolute error, and status is the quadrature_status enum (AppleAccelerate.LibAccelerate.QUADRATURE_SUCCESS on success).

r = AppleAccelerate.integrate(x -> x^2, 0, 1)
@assert isapprox(r.value, 1/3; atol = 1e-8)
r.value, r.abserr

(0.33333333333333337, 3.700743415417189e-15)

AppleAccelerate.integrate(sin, 0, π).value   # ≈ 2

2.0

Integrators

The integrator keyword selects the QUADPACK routine:

`integrator`	Routine	Use case
`:qng`	Non-adaptive Gauss-Kronrod	Smooth integrands; fastest, fixed rule
`:qag`	Adaptive Gauss-Kronrod	General-purpose adaptive on a finite interval
`:qags`	Adaptive with extrapolation (default)	Endpoint singularities; supports infinite bounds

For :qag, the qag_points keyword selects the Gauss-Kronrod rule and must be one of 15, 21, 31, 41, 51, 61 (or 0 for the library default). Any other value throws an ArgumentError (the underlying library would otherwise silently return value = 0.0).

AppleAccelerate.integrate(x -> exp(-x), 0, 5; integrator = :qag, qag_points = 21).value

0.9932620530009145

Infinite and reversed bounds

For integrator = :qags, one or both bounds may be infinite (±Inf):

r = AppleAccelerate.integrate(x -> exp(-x^2), -Inf, Inf)
@assert isapprox(r.value, sqrt(π); atol = 1e-6)
r.value   # ≈ √π

1.7724538509055159

Reversed bounds (a > b) follow the usual sign convention, ∫ₐᵇ f = -∫ᵦᵃ f:

fwd = AppleAccelerate.integrate(sin, 0, π).value
rev = AppleAccelerate.integrate(sin, π, 0).value
@assert isapprox(fwd, -rev; atol = 1e-8)

Tolerances and subdivision

abstol and reltol set the requested absolute and relative error; max_intervals bounds the number of subintervals the adaptive routines may create. If the integrand throws, that exception is captured and rethrown from integrate (rather than returning a meaningless value).

AppleAccelerate.integrate — Function

integrate(f, a, b; integrator=:qags, abstol=1e-8, reltol=1e-8,
          max_intervals=200, qag_points=0)

Numerically integrate the real scalar function f over the interval (a, b) using Apple's Accelerate Quadrature library (a C port of QUADPACK).

f is an ordinary Julia function f(x::Float64) -> Float64; Accelerate calls it in batches, which is handled internally.

For integrator=:qags, one or both bounds may be infinite (±Inf).

For integrator=:qag, qag_points selects the Gauss-Kronrod rule and must be one of 15, 21, 31, 41, 51, 61 (or 0 for the library default). Any other value throws an ArgumentError (the underlying library would otherwise silently return value = 0.0).

Returns a named tuple (value, abserr, status) where value is the integral estimate, abserr is the estimated absolute error, and status is the quadrature_status enum (QUADRATURE_SUCCESS on success).

Examples

integrate(x -> x^2, 0, 1).value          # ≈ 0.3333333
integrate(sin, 0, π).value               # ≈ 2.0
integrate(x -> exp(-x^2), -Inf, Inf).value  # ≈ √π

source