Three-dimensional banded tensor with element type T.

Extended help

Band structure

The band structure is assumed to be symmetric, and is defined in terms of the band width $b$. For a cubic banded tensor of dimensions $N × N × N$, the element $A_{ijk}$ may be non-zero only if $|i - j| ≤ b$, $|i - k| ≤ b$ and $|j - k| ≤ b$.


The data is stored as a Vector of small matrices, each with size $r × r$, where $r = 2b + 1$ is the total number of bands. Each submatrix holds the non-zero values of a slice of the form A[:, :, k].

For $b = 2$, one of these matrices looks like the following, where dots indicate out-of-bands values (equal to zero):

| x  x  x  ⋅  ⋅ |
| x  x  x  x  ⋅ |
| x  x  x  x  x |
| ⋅  x  x  x  x |
| ⋅  ⋅  x  x  x |

These submatrices are stored as static matrices (SMatrix).

Setting elements

To define the elements of the tensor, each slice A[:, :, k] must be set at once. For instance:

A = BandedTensor3D(undef, 20, Val(2))  # tensor of size 20×20×20 and band width b = 2
for k in axes(A, 3)
    A[:, :, k] = rand(5, 5)

See setindex! for more details.

Non-cubic tensors

A slight departure from cubic tensors is currently supported, with dimensions of the form $N × N × M$. Moreover, bands may be shifted along the third dimension by an offset $δ$. In this case, the bands are given by $|i - j| ≤ b$, $|i - (k + δ)| ≤ b$ and $|j - (k + δ)| ≤ b$.

BandedTensor3D{T}(undef, (Ni, Nj, Nk), Val(b); [bandshift = (0, 0, 0)])
BandedTensor3D{T}(undef, N, Val(b); ...)

Construct 3D banded tensor with band widths b.

Right now, the first two dimension sizes Ni and Nj of the tensor must be equal. In the second variant, the tensor dimensions are N × N × N.

The tensor is constructed uninitialised. Each submatrix A[:, :, k] of size (2b + 1, 2b + 1), for k ∈ 1:Nk, should be initialised as in the following example:

A[:, :, k] = rand(2b + 1, 2b + 1)

The optional bandshift argument should be a tuple of the form (δi, δj, δk) describing a band shift. Right now, band shifts are limited to δi = δj = 0, so this argument should actually look like (0, 0, δk).

muladd!(Y::AbstractMatrix, A::BandedTensor3D, b::AbstractVector)

Perform contraction Y[i, j] += ∑ₖ A[i, j, k] * b[k].

Note that the result is added to previously existent values of Y.

As an (allocating) alternative, one can use Y = A * b, which returns Y as a BandedMatrix.


Get band width b of BandedTensor3D.

The band width is defined here such that the element A[i, j, k] may be non-zero only if $|i - j| ≤ b$, $|i - k| ≤ b$ and $|j - k| ≤ b$. This definition is consistent with the specification of the upper and lower band widths in BandedMatrices.

setindex!(A::BandedTensor3D, Ak::AbstractMatrix, :, :, k)

Set submatrix A[:, :, k] to the matrix Ak.

The Ak matrix must have dimensions (r, r), where r = 2b + 1 is the total number of bands of A.

dot(x, Asub::SubMatrix, y)

Efficient implementation of the generalised dot product dot(x, Asub * y).

To be used with a submatrix Asub = A[:, :, k] of a BandedTensor3D A.

RecombineMatrix{T} <: AbstractMatrix{T}

Matrix for transformation from coefficients of the recombined basis, to the corresponding B-spline basis coefficients.

Extended help

The transformation matrix $M$ is defined by

\[ϕ_j = M_{ij} b_i,\]

where $b_j(x)$ and $ϕ_i(x)$ are elements of the B-spline and recombined bases, respectively.

This matrix allows to pass from known coefficients $u_j$ in the recombined basis $ϕ_j$, to the respective coefficients $v_i$ in the B-spline basis $b_i$:

\[\bm{v} = \mathbf{M} \bm{u}.\]

Note that the matrix is not square: it has dimensions $N × M$, where $N$ is the length of the B-spline basis, and $M = N - δ$ is that of the recombined basis (see RecombinedBSplineBasis for details).

Due to the local support of B-splines, basis recombination can be performed by combining just a small set of B-splines near the boundaries (as discussed in RecombinedBSplineBasis). This leads to a recombination matrix which is almost a diagonal of ones, plus a few extra super- and sub-diagonal elements in the upper left and lower right corners, respectively. The matrix is stored in a memory-efficient way that also allows fast access to its elements.

Efficient implementations of matrix-vector products (using the * operator or LinearAlgebra.mul!) and of left division of vectors (using \ or LinearAlgebra.ldiv!) are included. These two operations can be used to transform between coefficients in the original and recombined bases.

Note that, since the recombined basis forms a subspace of the original basis (which is related to the rectangular shape of the matrix), it is generally not possible to obtain recombined coefficients from original coefficients, unless the latter already satisfy the constraints encoded in the recombined basis. The left division operation will throw a NoUniqueSolutionError if that is not the case.

RecombineMatrix(ops::Tuple{Vararg{AbstractDifferentialOp}}, B::BSplineBasis, [T])
RecombineMatrix(ops_left, ops_right, B::BSplineBasis, [T])

Construct recombination matrix describing a B-spline basis recombination.

In the first case, ops is the boundary condition (BC) to be applied on both boundaries. The second case allows to set different BCs on each boundary.

The default element type T is generally Float64, except for specific differential operators which yield a matrix of zeroes and ones, for which Bool is the default.

See the RecombinedBSplineBasis constructor for details on the ops argument.

RecombinedBSplineBasis{k, T}

Functional basis defined from the recombination of a BSplineBasis in order to satisfy certain homogeneous boundary conditions (BCs).

Extended help

The basis recombination technique is a common way of applying BCs in Galerkin methods. It is described for instance in Boyd 2000 (ch. 6), in the context of a Chebyshev basis. In this approach, the original basis, $\{b_i(x), 1 ≤ i ≤ N\}$, is "recombined" into a new basis, $\{ϕ_j(x), 1 ≤ j ≤ M\}$, so that each basis function $ϕ_j$ individually satisfies the chosen BCs.

The length $M$ of the recombined basis is always smaller than the length $N$ of the original basis. The difference, $δ = N - M$, is equal to the number of boundary conditions. In the simplest (and most common) case, a single BC is applied on each boundary, leading to $δ = 2$. More generally, as described further below, it is possible to simultaneously impose different BCs, which further decreases the number of degrees of freedom (increasing $δ$).

Thanks to the local support of B-splines, basis recombination involves only a little portion of the original B-spline basis. For instance, since there is only one B-spline that is non-zero at each boundary, removing that function from the basis is enough to apply homogeneous Dirichlet BCs. Imposing BCs for derivatives is slightly more complex, but still possible.

Note that, when combining basis recombination with collocation methods, there must be no collocation points at the boundaries, or the resulting collocation matrices may not be invertible.

Order of the boundary condition

In this section, we consider the simplest case where a single homogeneous boundary condition, $\mathrm{d}^n u / \mathrm{d}x^n = 0$, is to be satisfied by the basis.

The recombined basis requires the specification of a Derivative object determining the order of the homogeneous BCs to be applied at the two boundaries. Linear combinations of Derivatives are also supported. The order of the B-spline needs to be $k ≥ n + 1$, since a B-spline of order $k$ is a $C^{k - 1}$-continuous function (except on the knots where it is $C^{k - 1 - p}$, with $p$ the knot multiplicity).

Some usual choices are:

  • Derivative(0) sets homogeneous Dirichlet BCs ($u = 0$ at the boundaries) by removing the first and last B-splines, i.e. $ϕ_1 = b_2$;

  • Derivative(1) sets homogeneous Neumann BCs ($u' = 0$ at the boundaries) by adding the two first (and two last) B-splines, i.e. $ϕ_1 = b_1 + b_2$.

  • more generally, α Derivative(0) + β Derivative(1) sets homogeneous Robin BCs by defining $ϕ_1$ as a linear combination of $b_1$ and $b_2$. Here it's important to note that Derivative(1) denotes the normal derivative at the boundary, $\frac{∂u}{∂n}$, which is equal to $-\frac{∂u}{∂x}$ on the left boundary.

Higher order BCs are also possible. For instance, Derivative(2) recombines the first three B-splines into two basis functions that satisfy $ϕ_1'' = ϕ_2'' = 0$ at the left boundary, while ensuring that lower and higher-order derivatives keep degrees of freedom at the boundary. Note that simply adding the first three B-splines, as in $ϕ_1 = b_1 + b_2 + b_3$, makes the first derivative vanish as well as the second one, which is unwanted.

For Derivative(2), the chosen solution is to set $ϕ_i = α_i b_i + β_i b_{i + 1}$ for $i ∈ \{1, 2\}$. The $α_i$ and $β_i$ coefficients are chosen such that $ϕ_i'' = 0$ at the boundary. Moreover, they satisfy the (somewhat arbitrary) constraint $α_i + β_i = 2$ for each $i$, for consistency with the Neumann case described above. This generalises to higher order BCs. Note that, since each boundary function $ϕ_i$ is defined from only two neighbouring B-splines, its local support stays minimal, hence preserving the small bandwidth of the resulting matrices.

Finally, note that in the current implementation, it is not possible to impose different boundary conditions on both boundaries.

Multiple boundary conditions

As an option, the recombined basis may simultaneously satisfy homogeneous BCs of different orders. In this case, a tuple of Derivatives must be passed.

For more details on the supported combinations of BCs, see the different RecombinedBSplineBasis constructors documented further below.

RecombinedBSplineBasis(B::BSplineBasis, op::AbstractDifferentialOp)
RecombinedBSplineBasis(B::BSplineBasis, op_left, op_right)

Construct RecombinedBSplineBasis from B-spline basis B, satisfying homogeneous boundary conditions (BCs) associated to the given differential operator.

The second case allows setting different BCs on each boundary.

For instance, op = Derivative(0) and op = Derivative(1) correspond to homogeneous Dirichlet and Neumann BCs, respectively.

Linear combinations of differential operators are also supported. For instance, op = Derivative(0) + λ Derivative(1) corresponds to homogeneous Robin BCs.

Higher-order derivatives are also allowed, being only limited by the order of the B-spline basis.


julia> B = BSplineBasis(BSplineOrder(4), -1:0.2:1)
13-element BSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]

julia> R_neumann = RecombinedBSplineBasis(B, Derivative(1))
11-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 BCs left:  (D{1},)
 BCs right: (D{1},)

julia> R_robin = RecombinedBSplineBasis(B, Derivative(0) + 3Derivative(1))
11-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 BCs left:  (D{0} + 3 * D{1},)
 BCs right: (D{0} + 3 * D{1},)

RecombinedBSplineBasis(B::BSplineBasis, ops)
RecombinedBSplineBasis(B::BSplineBasis, ops_left, ops_right)

Construct RecombinedBSplineBasis simultaneously satisfying homogeneous BCs associated to multiple differential operators.

Currently, the following cases are supported:

  1. all derivatives up to order m:

     ops = (Derivative(0), ..., Derivative(m))

    This boundary condition is obtained by removing the first m + 1 B-splines from the original basis.

    For instance, if (Derivative(0), Derivative(1)) is passed, then the basis simultaneously satisfies homogeneous Dirichlet and Neumann BCs at the two boundaries. The resulting basis is $ϕ_1 = b_3, ϕ_2 = b_4, …, ϕ_{N - 4} = b_{N - 2}$.

  2. an extension of the above, with an additional differential operator of order n at the end:

     ops = (Derivative(0), ..., Derivative(m), D(n))

    The operator D(n) may be a Derivative, or a linear combination of derivatives. The only restriction is that its maximum degree must satisfy n ≥ m + 1.

    One example is the combination of homogeneous Dirichlet BCs, $u = 0$ on the boundaries, with Robin BCs for the derivative, $u' + λ u'' = 0$, which corresponds to ops = (Derivative(0), Derivative(1) + λ Derivative(2)).

  3. generalised natural boundary conditions:

     ops = Natural()

    This is equivalent to ops = (Derivative(2), Derivative(3), ..., Derivative(k ÷ 2)) where k is the spline order (which must be even). See Natural for details.

In the first two cases, the degrees of the differential operators must be in increasing order. For instance, ops = (Derivative(1), Derivative(0)) fails with an error.


julia> ops = (Derivative(0), Derivative(1));

julia> R1 = RecombinedBSplineBasis(B, ops)
9-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 BCs left:  (D{0}, D{1})
 BCs right: (D{0}, D{1})

julia> ops = (Derivative(0), Derivative(1) - 4Derivative(2));

julia> R2 = RecombinedBSplineBasis(B, ops)
9-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 BCs left:  (D{0}, D{1} + -4 * D{2})
 BCs right: (D{0}, D{1} + -4 * D{2})

julia> R3 = RecombinedBSplineBasis(B, Natural())
11-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 BCs left:  (D{2},)
 BCs right: (D{2},)
constraints(R::AbstractBSplineBasis) -> (left, right)
constraints(A::RecombineMatrix) -> (left, right)

Return the constraints (homogeneous boundary conditions) that the basis satisfies on each boundary.

Constraints are returned as a tuple (left, right) indicating the BCs that are satisfied on each boundary. Each element is a tuple of differential operators specifying the BCs.

For example, if both Dirichlet and Neumann BCs are satisfied on the left boundary, then left = (Derivative(0), Derivative(1)).

For non-recombined bases such as BSplineBasis, this returns a tuple of empty tuples: ((), ()), since no BCs are satisfied on either boundary.

num_constraints(R::AbstractBSplineBasis) -> (Int, Int)
num_constraints(A::RecombineMatrix) -> (Int, Int)

Returns the number of constraints (number of BCs to satisfy) on each boundary.

For instance, if R simultaneously satisfies Dirichlet and Neumann boundary conditions on each boundary, this returns (2, 2).

Note that for non-recombined bases such as BSplineBasis, the number of constraints is zero, and this returns (0, 0).

num_recombined(R::AbstractBSplineBasis) -> (Int, Int)
num_recombined(A::RecombineMatrix) -> (Int, Int)

Returns the number of recombined functions in the recombined basis for each boundary.

For instance, if R satisfies Neumann boundary conditions on both boundaries, then only the first and last basis functions are different from the original B-spline basis, e.g. $ϕ_1 = b_1 + b_2$, and this returns (1, 1).

For non-recombined bases such as BSplineBasis, this returns zero.

nzrows(A::RecombineMatrix, col::Integer) -> UnitRange{Int}

Returns the range of row indices i such that A[i, col] is non-zero.

parent_coefficients(R::RecombinedBSplineBasis, coefs::AbstractVector)

Returns the coefficients associated to the parent B-spline basis, from the coefficients coefs in the recombined basis.

Note that this function doesn't allocate, since it returns a lazy concatenation of two StaticArrays and a view of the coefs vector.


Returns the number of functions in the recombined basis.


Get original B-spline basis.

    [derivatives = (Derivative(0), Derivative(0))],
    [MatrixType = BandedMatrix{Float64}];
    [quadrature = default_quadrature((B, B))],

Compute Galerkin mass or stiffness matrix, as well as more general variants of these.

Extended help

The Galerkin mass matrix is defined as

\[M_{ij} = ⟨ ϕ_i, ϕ_j ⟩ \quad \text{for} \quad i ∈ [1, N] \text{ and } j ∈ [1, N],\]

where $ϕ_i(x)$ is the $i$-th basis function and N = length(B) is the number of functions in the basis B. Here, $⟨f, g⟩$ is the $L^2$ inner product between functions $f$ and $g$.

Since products of B-splines are themselves piecewise polynomials, integrals can be computed exactly using Gaussian quadrature rules. To do this, we use Gauss–Legendre quadratures via the FastGaussQuadrature package.

Matrix layout and types

The mass matrix is banded with $2k - 1$ bands. Moreover, the matrix is symmetric and positive definite, and only $k$ bands are needed to fully describe the matrix. Hence, a Hermitian view of an underlying matrix is returned.

By default, the underlying matrix holding the data is a BandedMatrix that defines the upper part of the symmetric matrix. Other types of container are also supported, including regular sparse matrices (SparseMatrixCSC) and dense arrays (Matrix). See collocation_matrix for a discussion on matrix types.

Periodic B-spline bases

The default matrix type is BandedMatrix, except for periodic bases (PeriodicBSplineBasis), in which case the Galerkin matrix has a few out-of-bands entries due to periodicity. For periodic bases, SparseMatrixCSC is the default. Note that this may change in the future.

Derivatives of basis functions

Galerkin matrices associated to the derivatives of basis functions may be constructed using the optional derivatives parameter. For instance, if derivatives = (Derivative(0), Derivative(2)), the matrix $⟨ ϕ_i, ϕ_j'' ⟩$ is constructed, where primes denote derivatives. Note that, if the derivative orders are different, the resulting matrix is not symmetric, and a Hermitian view is not returned in those cases.

Combining different bases

More generally, it is possible to compute matrices of the form $⟨ ψ_i^{(n)}, ϕ_j^{(m)} ⟩$, where n and m are derivative orders, and $ψ_i$ and $ϕ_j$ belong to two different (but related) bases B₁ and B₂. For this, instead of the B parameter, one must pass a tuple of bases (B₁, B₂). The restriction is that the bases must have the same parent B-spline basis. That is, they must share the same set of B-spline knots and be of equal polynomial order.

Note that, if both bases are different, the matrix will not be symmetric, and will not even be square if the bases have different lengths.

In practice, this feature may be used to combine a B-spline basis B, with a recombined basis R generated from B (see Basis recombination).

Integrating more general functions

Say we wanted to integrate the more general term

\[L_{ij} = ⟨ f(x) ϕ_i(x), ϕ_j(x) ⟩ = ∫_{a}^{b} f(x) ϕ_i(x) ϕ_j(x) \, \mathrm{d}x\]

To obtain an approximation of this matrix, one would pass the $f(x)$ function as a first positional argument to galerkin_matrix (or galerkin_matrix!). This can also be combined with the derivatives argument if one wants to consider derivatives of $ϕ_i$ or $ϕ_j$.

Note that, in this case, the computation of the integrals is not guaranteed to be exact, since Gauss–Legendre quadratures are only "exact" when the integrand is a polynomial (the product of two B-splines is a piecewise polynomial). To improve accuracy, one may want to increase the number n of quadrature nodes. For this, pass quadrature = Galerkin.gausslegendre(Val(n)) as a keyword argument.

    [f::Function], A::AbstractMatrix, B::AbstractBSplineBasis, [deriv = (Derivative(0), Derivative(0))];

Fill preallocated Galerkin matrix.

The matrix may be a Hermitian view, in which case only one half of the matrix will be filled. Note that, for the matrix to be symmetric, both derivative orders in deriv must be the same.

More generally, it is possible to combine different functional bases by passing a tuple of AbstractBSplineBasis as B.

See galerkin_matrix for details on arguments.

    f, B::AbstractBSplineBasis,
    [deriv = Derivative(0)], [VectorType = Vector{Float64}],

Perform Galerkin projection of a function f onto the given basis.

By default, returns a vector with values

\[φ_i = ⟨ b_i, f ⟩ = ∫_a^b b_i(x) \, f(x) \, \mathrm{d}x,\]

where $a$ and $b$ are the boundaries of the B-spline basis $\{ b_i \}_{i = 1}^N$.

The integrations are performed using Gauss–Legendre quadrature. The number of quadrature nodes is chosen so that the result is exact when $f$ is a polynomial of degree $k - 1$ (or, more generally, a spline belonging to the space spanned by the basis B). Here $k$ is the order of the B-spline basis. In the more general case, this function returns a quadrature approximation of the projection.

See also galerkin_projection! for the in-place operation, and galerkin_matrix for more details.

    (D₁::Derivative, D₂::Derivative, D₃::Derivative),
    [T = Float64],

Compute 3D banded tensor appearing from quadratic terms in Galerkin method.

As with galerkin_matrix, it is also possible to combine different functional bases by passing, instead of B, a tuple (B₁, B₂, B₃) of three AbstractBSplineBasis. For now, the first two bases, B₁ and B₂, must have the same length.

The tensor is efficiently stored in a BandedTensor3D object.

Flat <: AbstractExtrapolationMethod

Represents a flat extrapolation: spline values at domain limits are extended to the left and to the right.


Spline interpolation.

This is the type returned by interpolate.

A SplineInterpolation I can be evaluated at any point x using the I(x) syntax.

It can also be updated with new data on the same data points using interpolate!.

SplineInterpolation(undef, B::AbstractBSplineBasis, x::AbstractVector, [T = eltype(x)])

Initialise a SplineInterpolation from B-spline basis and a set of interpolation (or collocation) points x.

Note that the length of x must be equal to the number of B-splines.

Use interpolate! to actually interpolate data known on the x locations.

interpolate(x, y, BSplineOrder(k), [bc = nothing])

Interpolate values y at locations x using B-splines of order k.

Grid points x must be real-valued and are assumed to be in increasing order.

Returns a SplineInterpolation which can be evaluated at any intermediate point.

Optionally, one may pass one of the boundary conditions listed in the Boundary conditions section. Currently, the Natural and Periodic boundary conditions are available.

See also interpolate!.

Periodic boundary conditions

Periodic boundary conditions should be used if the interpolated data is supposed to represent a periodic signal. In this case, pass bc = Period(L), where L is the period of the x-axis. Note that the endpoint x[begin] + L should not be included in the x vector.

Cubic periodic splines

Cubic periodic splines (BSplineOrder(4)) are particularly well optimised compared to periodic splines of other orders. Just note that interpolations using cubic periodic splines modify their input (including x and y values).


julia> xs = -1:0.1:1;

julia> ys = cospi.(xs);

julia> S = interpolate(xs, ys, BSplineOrder(4))
SplineInterpolation containing the 21-element Spline{Float64}:
 basis: 21-element BSplineBasis of order 4, domain [-1.0, 1.0]
 order: 4
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3  …  0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0, 1.0, 1.0, 1.0]
 coefficients: [-1.0, -1.00111, -0.8975, -0.597515, -0.314147, 1.3265e-6, 0.314142, 0.597534, 0.822435, 0.96683  …  0.96683, 0.822435, 0.597534, 0.314142, 1.3265e-6, -0.314147, -0.597515, -0.8975, -1.00111, -1.0]
 interpolation points: -1.0:0.1:1.0

julia> S(-1)

julia> (Derivative(1) * S)(-1)

julia> (Derivative(2) * S)(-1)

julia> Snat = interpolate(xs, ys, BSplineOrder(4), Natural())
SplineInterpolation containing the 21-element Spline{Float64}:
 basis: 21-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0], BCs {left => (D{2},), right => (D{2},)}
 order: 4
 knots: [-1.0, -1.0, -1.0, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4  …  0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.0, 1.0, 1.0]
 coefficients: [-0.833333, -0.647516, -0.821244, -0.597853, -0.314057, -2.29076e-5, 0.314148, 0.597532, 0.822435, 0.96683  …  0.96683, 0.822435, 0.597532, 0.314148, -2.29076e-5, -0.314057, -0.597853, -0.821244, -0.647516, -0.833333]
 interpolation points: -1.0:0.1:1.0

julia> Snat(-1)

julia> (Derivative(1) * Snat)(-1)

julia> (Derivative(2) * Snat)(-1)

Periodic boundary conditions

Interpolate $f(x) = \cos(πx)$ for $x ∈ [-1, 1)$. Note that the period is $L = 2$ and that the endpoint ($x = 1$) must not be included in the data points.

julia> xp = -1:0.1:0.9;

julia> yp = cospi.(xp);

julia> Sper = interpolate(xp, yp, BSplineOrder(4), Periodic(2))
SplineInterpolation containing the 20-element Spline{Float64}:
 basis: 20-element PeriodicBSplineBasis of order 4, domain [-1.0, 1.0), period 2.0
 order: 4
 knots: [..., -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3  …  0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, ...]
 coefficients: [..., -1.01659, -0.96683, -0.822435, -0.597534, -0.314142, 1.10589e-17, 0.314142, 0.597534, 0.822435, 0.96683, 1.01659, 0.96683, 0.822435, 0.597534, 0.314142, 3.90313e-18, -0.314142, -0.597534, -0.822435, -0.96683, ...]
 interpolation points: -1.0:0.1:0.9

As expected, the periodic spline does a better job at approximating the periodic function $f(x) = \cos(πx)$ near the boundaries than the other interpolations:

julia> x = -0.99; cospi(x), Sper(x), Snat(x), S(x)
(-0.9995065603657316, -0.9995032595823043, -0.9971071640321146, -0.9996420091470221)

julia> x = 0.998; cospi(x), Sper(x), Snat(x), S(x)
(-0.9999802608561371, -0.9999801044078943, -0.9994253145274461, -1.0000122303614758)

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.

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.


Specifies an approximation by interpolation at the given points xs.


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>,\]


\[\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$.

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.


Currently, this method is not guaranteed to work well with recombined B-spline bases (of type RecombinedBSplineBasis).

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:

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.


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)
PeriodicVector{T} <: AbstractVector{T}

Describes a periodic (or "circular") vector wrapping a regular vector.

Used to store the coefficients of periodic splines.

The vector has an effective length N associated to a single period, but it is possible to index it outside of this "main" interval.

This is similar to BSplines.PeriodicKnots. It is simpler though, since here there is no notion of coordinates or of a period L. Periodicity is only manifest in the indexation of the vector, e.g. a PeriodicVector vs satisfies vs[i + N] == vs[i].


Wraps coefficient vector cs such that it can be indexed in a periodic manner.


Represents a spline function.

Spline(B::AbstractBSplineBasis, coefs::AbstractVector)

Construct a spline from a B-spline basis and a vector of B-spline coefficients.


julia> B = BSplineBasis(BSplineOrder(4), -1:0.2:1);

julia> coefs = rand(length(B));

julia> S = Spline(B, coefs)
13-element Spline{Float64}:
 basis: 13-element BSplineBasis of order 4, domain [-1.0, 1.0]
 order: 4
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 coefficients: [0.173575, 0.321662, 0.258585, 0.166439, 0.527015, 0.483022, 0.390663, 0.802763, 0.721983, 0.372347, 0.0301856, 0.0793339, 0.663758]

Spline{T = Float64}(undef, B::AbstractBSplineBasis)

Construct a spline with uninitialised vector of coefficients.


Evaluate spline at coordinate x.

SplineWrapper{S <: Spline}

Abstract type representing a type that wraps a Spline.

Such a type implements all common operations on splines, including evaluation, differentiation, etc…


Returns an antiderivative of the given spline as a new spline.

The algorithm is described in de Boor 2001, p. 127.

Note that the integral spline I returned by this function is defined up to a constant. By convention, here the returned spline I is zero at the left boundary of the domain. One usually cares about the integral of S from point a to point b, which can be obtained as I(b) - I(a).

Periodic splines

Note that the integral of a periodic function is in general not periodic. For periodic splines (backed by a PeriodicBSplineBasis), this function returns a non-periodic spline (backed by a regular BSplineBasis).

*(op::Derivative, S::Spline) -> Spline

Returns N-th derivative of spline S as a new spline.

See also diff.


julia> B = BSplineBasis(BSplineOrder(4), -1:0.2:1);

julia> S = Spline(B, rand(length(B)))
13-element Spline{Float64}:
 basis: 13-element BSplineBasis of order 4, domain [-1.0, 1.0]
 order: 4
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 coefficients: [0.173575, 0.321662, 0.258585, 0.166439, 0.527015, 0.483022, 0.390663, 0.802763, 0.721983, 0.372347, 0.0301856, 0.0793339, 0.663758]

julia> Derivative(0) * S === S

julia> Derivative(1) * S
12-element Spline{Float64}:
 basis: 12-element BSplineBasis of order 3, domain [-1.0, 1.0]
 order: 3
 knots: [-1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0]
 coefficients: [2.22131, -0.473071, -0.460734, 1.80288, -0.219964, -0.461794, 2.0605, -0.403899, -1.74818, -1.71081, 0.368613, 8.76636]

julia> Derivative(2) * S
11-element Spline{Float64}:
 basis: 11-element BSplineBasis of order 2, domain [-1.0, 1.0]
 order: 2
 knots: [-1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0]
 coefficients: [-26.9438, 0.0616849, 11.3181, -10.1142, -1.20915, 12.6114, -12.322, -6.72141, 0.186876, 10.3971, 83.9775]
diff(S::Spline, [op::Derivative = Derivative(1)]) -> Spline

Same as op * S.

Returns N-th derivative of spline S as a new spline.


Returns type of element returned when evaluating the Spline.


Returns the number of coefficients in the spline.

Note that this is equal to the number of basis functions, length(basis(S)).


Module defining B-spline bases and B-spline related functions.


Abstract type defining a B-spline basis, or more generally, a functional basis defined from B-splines.

The basis is represented by a B-spline order k and a knot element type T.

    x::Real, [op = Derivative(0)], [T = float(typeof(x))];
    [ileft = nothing],
) -> (i, bs)

Evaluates all basis functions which are non-zero at x.

This is a convenience alias for evaluate_all. See evaluate_all for details on optional arguments and on the returned values.


julia> B = BSplineBasis(BSplineOrder(4), -1:0.1:1)
23-element BSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4  …  0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.0, 1.0, 1.0]

julia> i, bs = B(0.42)
(18, (0.0013333333333333268, 0.28266666666666657, 0.6306666666666667, 0.08533333333333339))

julia> sum(bs)

julia> bs[1] - B[i](0.42)

julia> bs[2] - B[i - 1](0.42)

julia> B(0.44; ileft = i)
(18, (0.01066666666666666, 0.4146666666666667, 0.5386666666666665, 0.03599999999999999))

julia> B(0.42, Float32)
(18, (0.0013333336f0, 0.28266668f0, 0.6306667f0, 0.085333325f0))

julia> B(0.42, Derivative(1))
(18, (0.19999999999999937, 6.4, -3.3999999999999977, -3.200000000000001))
AugmentedKnots{T,k} <: AbstractVector{T}

Pads from both sides a vector of B-spline breakpoints, making sure that the first and last values are repeated k times.

BSplineBasis{k, T}

B-spline basis for splines of order k and knot element type T <: Real.

The basis is defined by a set of knots and by the B-spline order.

BSplineBasis(order::BSplineOrder{k}, ts::AbstractVector; augment = Val(true))
BSplineBasis(order::BSplineOrder{k}, ts::NTuple; augment = Val(true))

Create B-spline basis of order k with breakpoints ts.

If augment = Val(true), breakpoints will be "augmented" so that boundary knots have multiplicity k. Note that, if they are passed as a regular Vector, the input may be modified. See augment_knots! for details.


julia> breaks = range(-1, 1; length = 21)

julia> B = BSplineBasis(BSplineOrder(5), breaks)
24-element BSplineBasis of order 5, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5  …  0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.0, 1.0, 1.0, 1.0]

Note that first and last knots are repeated $k = 5$ times.

If augment = Val(false), input breakpoints are taken without modification as the knots $t_i$ of the B-spline basis. Note that the valid domain is reduced to $[-0.6, 0.6]$. The domain is always defined as the range $[t_k, t_{N + 1}]$, where $N$ is the length of the basis (below, $N = 16$).

julia> Bn = BSplineBasis(5, breaks, augment = Val(false))
16-element BSplineBasis of order 5, domain [-0.6, 0.6]
 knots: -1.0:0.1:1.0

Statically-sized bases

To define a basis with static size (i.e. size known at compile time), the breakpoints ts should be passed as a tuple or as an SVector (from the StaticArrays package):

julia> breaks = (0.0, 0.1, 0.2, 0.6, 1.0);

julia> B = BSplineBasis(BSplineOrder(3), breaks)
6-element BSplineBasis of order 3, domain [0.0, 1.0]
 knots: [0.0, 0.0, 0.0, 0.1, 0.2, 0.6, 1.0, 1.0, 1.0]

julia> knots(B)
9-element StaticArraysCore.SVector{9, Float64} with indices SOneTo(9):
BasisFunction{B <: AbstractBSplineBasis, T}

Describes a single basis function.

The basis function may belong to a BSplineBasis (in which case it's effectively a B-spline), or to a basis derived from a B-spline basis (such as a RecombinedBSplineBasis).

BasisFunction(basis::AbstractBSplineBasis, i::Int, [T = Float64])

Construct i-th basis function of the given basis.

The constructed function can be evaluated as b(x), returning a value of type T.

(b::BasisFunction)(x, [op::AbstractDifferentialOp])

Evaluate basis function at coordinate x.

To evaluate a derivative, pass Derivative(n) as the op argument, with n the derivative order.

To evaluate multiple derivatives, pass a derivative range Derivative(m:n). In particular, Derivative(m:n) evaluates the basis function itself and its first n derivatives.

More general differential operators, such as Derivative(n) + λ Derivative(m), are also supported.

PeriodicBSplineBasis{k, T}

B-spline basis for splines of order k and knot element type T <: Real.

The basis is defined by a set of knots and by the B-spline order.

PeriodicBSplineBasis(order::BSplineOrder{k}, ts::AbstractVector)

Create periodic B-spline basis of order k with knots ts.

The knot period is taken to be L = ts[end] - ts[begin].


Create B-spline basis on periodic domain with period $L = 2$.

julia> ts = range(-1, 1; length = 11)

julia> B = PeriodicBSplineBasis(BSplineOrder(4), ts)
10-element PeriodicBSplineBasis of order 4, domain [-1.0, 1.0), period 2.0
 knots: [..., -1.4, -1.2, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, ...]

julia> period(B)

julia> length(B)

julia> boundaries(B)
(-1.0, 1.0)

julia> B(-0.42)
(5, (0.12150000000000002, 0.6571666666666667, 0.22116666666666668, 0.00016666666666666563))

julia> B(-0.42 + 2)
(15, (0.12150000000000015, 0.6571666666666667, 0.22116666666666657, 0.00016666666666666674))

julia> knots(B)
14-element PeriodicKnots{Float64, 4, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}:
PeriodicKnots{T} <: AbstractVector{T}

Describes an infinite vector of knots with periodicity L.

PeriodicKnots(ξs::AbstractVector{T}, ::BSplineOrder{k})

Construct a periodic knot sequence from breakpoints ξs.

The knot period is taken to be L = ξs[end] - ξs[begin]. In other words, the breakpoints ξs are expected to include the endpoint ξs[begin] + L.

The breakpoints should be given in non-decreasing order.

Note that the indices of the returned knots ts are offset with respect to the input ξs according to ts[i] = ξs[i + offset] where offset = k ÷ 2.

augment_knots!(breaks::AbstractVector, k::Union{Integer,BSplineOrder})

Modifies the input breakpoints to make sure that the first and last knot have multiplicity k for splines of order k.

To prevent allocations, this function will modify the input when this is a standard Vector. Otherwise, the input will be wrapped inside an AugmentedKnots object.

It is assumed that the input breakpoints have multiplicity 1 at the borders. That is, border coordinates should not be repeated in the input.

basis_to_array_index(B::AbstractBSplineBasis, axs, i::Int) -> Int

Converts from a basis index i (for basis function $b_i$) to a valid index in an array A with indices axs, where axs is typically the result of axes(A, dim) for a given dimension dim.

This is mainly relevant for periodic bases, in which case i may be outside of the "main" period of the basis.

It is assumed that length(B) == length(axs) (this is not checked by this function).


Returns (xmin, xmax) tuple with the boundaries of the domain supported by the basis.

common_support(b1::BasisFunction, b2::BasisFunction, ...) -> UnitRange{Int}

Get range of knots commonly supported by different basis functions.

If the supports don't intersect, an empty range is returned (e.g. 6:5), following the behaviour of intersect. The lack of intersection can be checked using isempty, which returns true for such a range.

evaluate(B::AbstractBSplineBasis, i::Integer, x,
         [op::AbstractDifferentialOp], [T=Float64])

Evaluate $i$-th basis function in the given basis at x (can be a coordinate or a vector of coordinates).

To evaluate a derivative, pass Derivative(n) as the op argument, with n the derivative order.

More general differential operators, such as Derivative(n) + λ Derivative(m), are also supported.

See also evaluate!.

evaluate!(b::AbstractVector, B::BSplineBasis, i::Integer,
          x::AbstractVector, args...)

Evaluate i-th basis function at positions x and write result to b.

See also evaluate.

    B::AbstractBSplineBasis, x::Real,
    [op = Derivative(0)], [T = float(typeof(x))];
    [ileft = nothing],
) -> i, bs

Evaluate all B-splines which are non-zero at coordinate x.

Returns a tuple (i, bs), where i is an index identifying the basis functions that were computed, and bs is a tuple with the actual values.

More precisely:

  • i is the index of the first B-spline knot $t_{i}$ when going from $x$ towards the left. In other words, it is such that $t_{i} ≤ x < t_{i + 1}$.

    It can be effectively computed as i = searchsortedlast(knots(B), x). If the correct value of i is already known, one can avoid this computation by manually passing this index via the optional ileft keyword argument.

  • bs is a tuple of B-splines evaluated at $x$:

    \[(b_i(x), b_{i - 1}(x), …, b_{i - k + 1}(x)).\]

    It contains $k$ values, where $k$ is the order of the B-spline basis. Note that values are returned in backwards order starting from the $i$-th B-spline.

Computing derivatives

One can pass the optional op argument to compute B-spline derivatives instead of the actual B-spline values.


See AbstractBSplineBasis for some examples using the alternative evaluation syntax B(x, [op], [T]; [ileft]), which calls this function.

find_knot_interval(ts::AbstractVector, x::Real, [ileft = nothing]) -> (i, zone)

Finds the index $i$ corresponding to the knot interval $[t_i, t_{i + 1}]$ that should be used to evaluate B-splines at location $x$.

The knot vector is assumed to be sorted in non-decreasing order.

It also returns a zone integer, which is:

  • 0 if x is within the knot domain (ts[begin] ≤ x ≤ ts[end]),
  • -1 if x < ts[begin],
  • 1 if x > ts[end].

This function is functionally equivalent to de Boor's INTERV routine (de Boor 2001, p. 74).

If one already knows the location i associated to the knot interval, then one can pass it as the optional ileft argument, in which case only the zone needs to be computed.

has_parent_basis(::Type{<:AbstractBSplineBasis}) -> Bool
has_parent_basis(::AbstractBSplineBasis) -> Bool

Trait determining whether a basis has a parent B-spline basis.

This is notably the case for RecombinedBSplineBasis, which are derived from regular B-spline bases.

isindomain(B::AbstractBSplineBasis, x::Real)

Check if the coordinate x is within the boundaries of the domain.


Returns the knots of the B-spline basis.

nonzero_in_segment(B::AbstractBSplineBasis, n::Int) -> UnitRange{Int}

Returns the range of basis functions that are non-zero in a given knot segment.

The $n$-th knot segment is defined by $Ω_n = [t_n, t_{n + 1}]$.

For BSplineBasis and RecombinedBSplineBasis, the number of non-zero functions in any given segment is generally equal to the B-spline order $k$. This number decreases near the borders, but this is not significant when B-spline knots have multiplicity $k$ there (as is the default).

For B-spline bases, and excepting the borders, the non-zero B-splines are $\left\{ b_i \right\}_{i = n - k + 1}^{n}$. This function thus returns (n - k + 1):N when B is a BSplineBasis.

See also support for the inverse operation.

order(::Type{AbstractBSplineBasis}) -> Int
order(::AbstractBSplineBasis) -> Int
order(::Type{Spline}) -> Int
order(::BSplineOrder) -> Int

Returns order of B-splines as an integer.

static_length(::Type{<:AbstractBSplineBasis}) -> Union{Int, Nothing}
static_length(::AbstractBSplineBasis) -> Union{Int, Nothing}

Return the basis' length if it is statically known (i.e. at compile time); return nothing otherwise.

Typically, bases with statically-known length are those constructed using an SVector (from the StaticArrays package) to describe the basis breakpoints.

support(B::AbstractBSplineBasis, i::Integer) -> UnitRange{Int}

Get range of knots supported by the $i$-th basis function.

support(b::BasisFunction) -> UnitRange{Int}

Get range of knots supported by the basis function.

Returns the knot range i:j such that the basis function support is $t ∈ [t_i, t_j)$.

getindex(B::AbstractBSplineBasis, i, [T = Float64])

Get $i$-th basis function.

This is an alias for BasisFunction(B, i, T) (see BasisFunction for details).

The returned object can be evaluated at any point within the boundaries defined by the basis.


julia> B = BSplineBasis(BSplineOrder(4), -1:0.1:1)
23-element BSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4  …  0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.0, 1.0, 1.0]

julia> B[6]
Basis function i = 6
  from 23-element BSplineBasis of order 4, domain [-1.0, 1.0]
  support: [-0.8, -0.4) (knots 6:10)

julia> B[6](-0.5)

julia> B[6, Float32](-0.5)

julia> B[6](-0.5, Derivative(1))
length(g::BSplineBasis) -> Int

Returns the number of B-splines composing a spline.

AvgKnots <: SelectionMethod

Each collocation point is chosen as a sliding average over k - 1 knots:

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

The resulting collocation points are sometimes called Greville sites (de Boor 2001).

CollocationMatrix{T} <: AbstractBandedMatrix{T}

B-spline collocation matrix, defined by

\[C_{ij} = b_j(x_i),\]

where $\bm{x}$ is a set of collocation points.

Provides an efficient LU factorisation without pivoting adapted from de Boor (1978). The factorisation takes advantage of the total positivity of spline collocation matrices (de Boor 2001, p. 175).


CollocationMatrix supports in-place LU factorisation using lu!, as well as out-of-place factorisation using lu. LU decomposition may also be performed via factorize.

SameAsKnots <: SelectionMethod

Collocation points are chosen to match B-spline knots.

Note that this only makes sense when the number of degrees of freedom of the B-spline basis (i.e. the length of the basis) matches the number of knots, which is generally not the case.

Some examples of bases that satisfy this condition are:

    [deriv::Derivative = Derivative(0)],
    [MatrixType = CollocationMatrix{Float64}];
    clip_threshold = eps(eltype(MatrixType)),

Return collocation matrix mapping B-spline coefficients to spline values at the collocation points x.

Extended help

The matrix elements are given by the B-splines evaluated at the collocation points:

\[C_{ij} = b_j(x_i) \quad \text{for} \quad i ∈ [1, N_x] \text{ and } j ∈ [1, N_b],\]

where Nx = length(x) is the number of collocation points, and Nb = length(B) is the number of B-splines in B.

To obtain a matrix associated to the B-spline derivatives, set the deriv argument to the order of the derivative.

Given the B-spline coefficients $\{u_j, 1 ≤ j ≤ N_b\}$, one can recover the values (or derivatives) of the spline at the collocation points as v = C * u. Conversely, if one knows the values $v_i$ at the collocation points, the coefficients $u$ of the spline passing by the collocation points may be obtained by inversion of the linear system u = C \ v.

The clip_threshold argument allows one to ignore spurious, negligible values obtained when evaluating B-splines. These values are typically unwanted, as they artificially increase the number of elements (and sometimes the bandwidth) of the matrix. They may appear when a collocation point is located on a knot. By default, clip_threshold is set to the machine epsilon associated to the matrix element type (see eps in the Julia documentation). Set it to zero to disable this behaviour.

Matrix types

The MatrixType optional argument allows to select the type of returned matrix.

Due to the compact support of B-splines, the collocation matrix is banded if the collocation points are properly distributed.

Supported matrix types

  • CollocationMatrix{T}: similar to a BandedMatrix{T}, with efficient LU factorisations without pivoting (see CollocationMatrix for details). This option performs much better than sparse matrices for inversion of linear systems. On the other hand, for matrix-vector or matrix-matrix multiplications, SparseMatrixCSC may perform better, especially when using OpenBLAS (see BandedMatrices issue). May fail with an error for non-square matrix shapes, or if the distribution of collocation points is not adapted. In these cases, the effective bandwidth of the matrix may be larger than the expected bandwidth.

  • SparseMatrixCSC{T}: regular sparse matrix; correctly handles all matrix shapes.

  • Matrix{T}: a regular dense matrix. Generally performs worse than the alternatives, especially for large problems.

Periodic B-spline bases

The default matrix type is CollocationMatrix{T}, except for periodic bases (PeriodicBSplineBasis), in which case the collocation matrix has a few out-of-bands entries due to periodicity. For cubic periodic bases, the Collocation.CyclicTridiagonalMatrix matrix type is used, which implements efficient solution of linear problems. For non-cubic periodic bases, SparseMatrixCSC is the default. Note that this may change in the future.

See also collocation_matrix!.

    method::SelectionMethod = default_method(B),

Define and return adapted collocation points for evaluation of splines.

The number of returned collocation points is equal to the number of functions in the basis.

Note that if B is a RecombinedBSplineBasis (adapted for boundary value problems), collocation points are not included at the boundaries, since the boundary conditions are implicitly satisfied by the basis.

In principle, the choice of collocation points is not unique. The selection method can be chosen via the method argument. For now, the following methods are accepted:

The former is the default, except for periodic B-spline bases (PeriodicBSplineBasis) of even order $k$, for which SameAsKnots is the default. (Note that for odd-order B-splines, this can lead to non-invertible collocation matrices.)

See also collocation_points!.

LinearAlgebra.ldiv!(x::AbstractVector, A::CyclicTridiagonalMatrix, y::AbstractVector)

Solve cyclic tridiagonal linear system Ax = y.

Note that all three arrays are modified by this function, with the result being stored in x.

Uses the algorithm described in Wikipedia based on the Sherman–Morrison formula.!Function!(C::CollocationMatrix, pivot = Val(false); check = true)

Perform in-place LU factorisation of collocation matrix without pivoting.

Takes advantage of the totally positive property of collocation matrices appearing in spline calculations (de Boor 1978).

The code is ported from Carl de Boor's BANFAC routine in FORTRAN77, via its FORTRAN90 version by John Burkardt.

LinearAlgebra.luMethod, pivot = NoPivot(); check = true)

Returns LU factorisation of collocation matrix.

See also lu!.


Module defining types describing differential operators and compositions thereof.

DerivativeUnitRange{m, n} <: AbstractDifferentialOp

Specifies a range of derivatives.


Two ways of constructing derivative ranges:

julia> Derivative(2):Derivative(4)

julia> Derivative(2:4)

julia> Tuple(Derivative(2:4))
(D{2}, D{3}, D{4})
LeftNormal <: AbstractNormalDirection

Specifies the normal direction on the left boundary of a 1D domain.

The left normal direction goes opposite to the coordinate axis.

RightNormal <: AbstractNormalDirection

Specifies the normal direction on the right boundary of a 1D domain.

The right normal direction is equal to that of the coordinate axis.

dot(op::AbstractDifferentialOp, dir::AbstractNormalDirection) -> AbstractDifferentialOp

Project derivative along a normal direction.

This should be used to convert from a normal derivative at the boundaries, to a derivative along the coordinate axes of the domain.

In practice, this returns op for RightNormal. For LeftNormal, it multiplies the odd-order derivatives by -1.

Natural <: BoundaryCondition

Generalised natural boundary condition.

This boundary condition is convenient for spline interpolations, as it provides extra constraints enabling to equate the number of unique B-spline knots to the number of data points.

For cubic splines (order $k = 4$), this corresponds to natural cubic splines, imposing the second derivatives to be zero at the boundaries ($S''(a) = S''(b) = 0$).

For higher-order splines, this boundary condition generalises the standard natural cubic splines, by setting derivatives of order $2, 3, …, k/2$ to be zero at the boundaries. For instance, for $k = 6$ (quintic splines), this imposes $S'' = S''' = 0$. In practice, BSplineKit.jl achieves this by using basis recombination.

Note that, for symmetry reasons, only even-order splines are supported by this BC.

Periodic <: BoundaryCondition

Represents periodic boundary conditions with a given period L.


Constructs periodic boundary conditions with period L.

period(bc::Periodic) -> Real
period(B::PeriodicBSplineBasis) -> Real
period(ts::PeriodicKnots) -> Real

Returns the period L associated to a periodic boundary condition or B-spline basis.