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

SciMLBase.SampledIntegralProblem — Type

Defines a integral problem over pre-sampled data. Documentation Page: https://docs.sciml.ai/Integrals/stable/

Mathematical Specification of a data Integral Problem

Sampled integral problems are defined as:

\[\sum_i w_i y_i\]

where y_i are sampled values of the integrand, and w_i are weights assigned by a quadrature rule, which depend on sampling points x.

Problem Type

Constructors

SampledIntegralProblem(y::AbstractArray, x::AbstractVector; dim=ndims(y), kwargs...)

y: The sampled integrand, must be a subtype of AbstractArray. It is assumed that the values of y along dimension dim correspond to the integrand evaluated at sampling points x
x: Sampling points, must be a subtype of AbstractVector.
dim: Dimension along which to integrate. Defaults to the last dimension of y.
kwargs: Keyword arguments copied to the solvers.

Fields

The fields match the names of the constructor arguments.

source

CommonSolve.solve — Method

solve(prob::SampledIntegralProblem, alg::SciMLBase.AbstractIntegralAlgorithm; kwargs...)

Keyword Arguments

There are no keyword arguments used to solve SampledIntegralProblems

source

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.12500292787107686

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.

Integrals.TrapezoidalRule — Type

TrapezoidalRule

Struct for evaluating an integral via the trapezoidal rule.

Example with sampled data:

using Integrals
f = x -> x^2
x = range(0, 1, length = 20)
y = f.(x)
problem = SampledIntegralProblem(y, x)
method = TrapezoidalRule()
solve(problem, method)

source

Integrals.SimpsonsRule — Type

SimpsonsRule

Struct for evaluating an integral via the Simpson's composite 1/3-3/8 rule over AbstractRanges (evenly spaced points) and Simpson's composite 1/3 rule for non-equidistant grids.

Example with equidistant data:

using Integrals
f = x -> x^2
x = range(0, 1, length = 20)
y = f.(x)
problem = SampledIntegralProblem(y, x)
method = SimpsonsRule()
solve(problem, method)

source