Function approximation

SplineApproximations

Module for function approximation using splines.

The general idea is to find the spline $g(x)$ in a given spline space that best approximates a known function $f(x)$. In other words, given a predefined BSplineBasis, the objective is to find some appropriate B-spline coefficients such that the resulting spline appropriately approximates $f$. Different approximation approaches are implemented, trading accuracy and speed.

approximate(f, B::AbstractBSplineBasis, [method = ApproxByInterpolation(B)])

Approximate function f in the given basis, using the chosen method.

From lower to higher accuracy (and cost), the possible approximation methods are:

• VariationDiminishing,
• ApproxByInterpolation,
• MinimiseL2Error.

See their respective documentations for details.

Note that, once an approximation has been performed, it's possible to efficiently perform additional approximations (of other functions f) by calling the in-place interpolate!. This completely avoids allocations and strongly reduces computation time.

Examples

julia> B = BSplineBasis(BSplineOrder(3), -1:0.4:1);

julia> S_interp = approximate(sin, B)
SplineApproximation containing the 7-element Spline{Float64}:
 basis: 7-element BSplineBasis of order 3, domain [-1.0, 1.0]
 order: 3
 knots: [-1.0, -1.0, -1.0, -0.6, -0.2, 0.2, 0.6, 1.0, 1.0, 1.0]
 coefficients: [-0.841471, -0.731727, -0.39727, 2.85767e-17, 0.39727, 0.731727, 0.841471]
 approximation method: interpolation at [-1.0, -0.8, -0.4, 0.0, 0.4, 0.8, 1.0]

julia> sin(0.3), S_interp(0.3)
(0.29552020666133955, 0.2959895327282942)

julia> approximate!(exp, S_interp)
SplineApproximation containing the 7-element Spline{Float64}:
 basis: 7-element BSplineBasis of order 3, domain [-1.0, 1.0]
 order: 3
 knots: [-1.0, -1.0, -1.0, -0.6, -0.2, 0.2, 0.6, 1.0, 1.0, 1.0]
 coefficients: [0.367879, 0.440373, 0.65701, 0.980127, 1.46223, 2.18111, 2.71828]
 approximation method: interpolation at [-1.0, -0.8, -0.4, 0.0, 0.4, 0.8, 1.0]

julia> exp(0.3), S_interp(0.3)
(1.3498588075760032, 1.3491015490105396)

julia> S_fast = approximate(exp, B, VariationDiminishing())
SplineApproximation containing the 7-element Spline{Float64}:
 basis: 7-element BSplineBasis of order 3, domain [-1.0, 1.0]
 order: 3
 knots: [-1.0, -1.0, -1.0, -0.6, -0.2, 0.2, 0.6, 1.0, 1.0, 1.0]
 coefficients: [0.367879, 0.449329, 0.67032, 1.0, 1.49182, 2.22554, 2.71828]
 approximation method: VariationDiminishing()

julia> S_opt = approximate(exp, B, MinimiseL2Error())
SplineApproximation containing the 7-element Spline{Float64}:
 basis: 7-element BSplineBasis of order 3, domain [-1.0, 1.0]
 order: 3
 knots: [-1.0, -1.0, -1.0, -0.6, -0.2, 0.2, 0.6, 1.0, 1.0, 1.0]
 coefficients: [0.368074, 0.440342, 0.657077, 0.980279, 1.46216, 2.18201, 2.71669]
 approximation method: MinimiseL2Error()

julia> x = 0.34; exp(x), S_opt(x), S_interp(x), S_fast(x)
(1.4049475905635938, 1.4044530324752076, 1.4044149581073813, 1.4328668494041878)

approximate!(f, A::SplineApproximation)

Approximate function f reusing a previous Spline approximation in a given basis.

See approximate for details.

AbstractApproxMethod

Abstract type describing a type of approximation method.

VariationDiminishing <: AbstractApproxMethod

Schoenberg's variation diminishing spline approximation.

Approximates a function $f$ by setting the B-spline coefficients as $c_i = f(x_i)$, where the locations $x_i$ are chosen as the Greville sites:

\[x_i = \frac{1}{k - 1} ∑_{j = 1}^{k - 1} t_{i + j}.\]

This method is expected to preserve the shape of the function. However, it may be very inaccurate as an actual approximation of it.

For details, see e.g. de Boor 2001, chapter XI.

ApproxByInterpolation <: AbstractApproxMethod

Approximate function by a spline that passes through the given set of points.

The number of points must be equal to the length of the B-spline basis defining the approximation space.

See belows for different ways of specifying the interpolation points.

ApproxByInterpolation(xs::AbstractVector)

Specifies an approximation by interpolation at the given points xs.

ApproxByInterpolation(B::AbstractBSplineBasis)

Specifies an approximation by interpolation at an automatically-determined set of points, which are expected to be appropriate for the given basis.

The interpolation points are determined by calling collocation_points.

MinimiseL2Error <: AbstractApproxMethod

Approximate a given function $f(x)$ by minimisation of the $L^2$ distance between $f$ and its spline approximation $g(x)$.

Extended help

Minimises the $L^2$ distance between the two functions:

\[{\left\lVert f - g \right\rVert}^2 = \left< f - g, f - g \right>,\]

where

\[\left< u, v \right> = ∫_a^b u(x) \, v(x) \, \mathrm{d}x\]

is the inner product between two functions, and $a$ and $b$ are the boundaries of the prescribed B-spline basis. Here, $g$ is the spline $g(x) = ∑_{i = 1}^N c_i \, b_i(x)$, and $\{ b_i \}_{i = 1}^N$ is a prescribed B-spline basis.

One can show that the optimal coefficients $c_i$ minimising the $L^2$ error are the solution to the linear system $\bm{M} \bm{c} = \bm{φ}$, where $M_{ij} = \left< b_i, b_j \right>$ and $φ_i = \left< b_i, f \right>$. These two terms are respectively computed by galerkin_matrix and galerkin_projection.

The integrals associated to $\bm{M}$ and $\bm{φ}$ are computed via Gauss–Legendre quadrature. The number of quadrature nodes is chosen as a function of the order $k$ of the prescribed B-spline basis, ensuring that $\bm{M}$ is computed exactly (see also galerkin_matrix). In the particular case where $f$ is a polynomial of degree $k - 1$, this also results in an exact computation of $\bm{φ}$. In more general cases, as long as $f$ is smooth enough, this is still expected to yield a very good approximation of the integral, and thus of the optimal coefficients $c_i$.