Integrating pre-sampled data

In some cases, instead of a function that acts as integrand, one only possesses a list of data points y at a set of sampling locations x, that must be integrated. This package contains functionality for doing that.

Example

Say, by some means we have generated a dataset x and y:

f = x -> x^2
x = range(0, 1, length = 20)
y = f.(x)

20-element Vector{Float64}:
 0.0
 0.0027700831024930744
 0.011080332409972297
 0.02493074792243767
 0.04432132963988919
 0.06925207756232686
 0.09972299168975068
 0.13573407202216065
 0.17728531855955676
 0.22437673130193903
 0.27700831024930744
 0.33518005540166207
 0.39889196675900274
 0.46814404432132967
 0.5429362880886426
 0.6232686980609419
 0.709141274238227
 0.8005540166204986
 0.8975069252077561
 1.0

Now, we can integrate this data set as follows:

problem = SampledIntegralProblem(y, x)
method = TrapezoidalRule()
solve(problem, method)

retcode: Success
u: 0.33379501385041543

The exact answer is of course $1/3$.

Details

SampledIntegralProblem and the generic solve interface are provided by SciMLBase and reexported by Integrals.jl for convenience. The following sections describe how Integrals.jl algorithms use sampled data.

Non-equidistant grids

If the sampling points x are provided as an AbstractRange (constructed with the range function for example), faster methods are used that take advantage of the fact that the points are equidistantly spaced. Otherwise, general methods are used for non-uniform grids.

Example:

f = x -> x^7
x = [0.0; sort(rand(1000)); 1.0]
y = f.(x)
problem = SampledIntegralProblem(y, x)
method = TrapezoidalRule()
solve(problem, method)

retcode: Success
u: 0.12500383831565015

Evaluating multiple integrals at once

If the provided data set y is a multidimensional array, the integrals are evaluated across only one of its axes. For performance reasons, the last axis of the array y is chosen by default, but this can be modified with the dim keyword argument to the problem definition.

f1 = x -> x^2
f2 = x -> x^3
f3 = x -> x^4
x = range(0, 1, length = 20)
y = [f1.(x) f2.(x) f3.(x)]
problem = SampledIntegralProblem(y, x; dim = 1)
method = SimpsonsRule()
solve(problem, method)

retcode: Success
u: 3-element Vector{Float64}:
 0.33333333333333326
 0.24999999999999997
 0.20000122504525594

Supported methods

Right now, only the TrapezoidalRule and SimpsonsRule are supported.

TrapezoidalRule

using Integrals

f = x -> x^2
x = range(0, 1, length = 20)
y = f.(x)
prob = SampledIntegralProblem(y, x)
sol = solve(prob, TrapezoidalRule())

SimpsonsRule

using Integrals

f = x -> x^2
x = range(0, 1, length = 21)
y = f.(x)
prob = SampledIntegralProblem(y, x)
sol = solve(prob, SimpsonsRule())