Collocation tools

collocation_points(
    B::AbstractBSplineBasis,
    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:

• Collocation.AvgKnots();
• Collocation.SameAsKnots(), which requires the length of the basis to be equal to the number of knots.

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

collocation_points!(
    x::AbstractVector, B::AbstractBSplineBasis,
    method::SelectionMethod = default_method(B),
)

Fill vector with collocation points for evaluation of splines.

See collocation_points for details.

SelectionMethod

Abstract type defining a method for choosing collocation points from spline knots.

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

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:

• recombined B-spline bases (RecombinedBSplineBasis) with Natural boundary conditions;

• periodic B-spline bases (PeriodicBSplineBasis).

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

Factorisation

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.

collocation_matrix(
    B::AbstractBSplineBasis,
    x::AbstractVector,
    [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.

See also collocation_matrix!.

collocation_matrix!(
    C::AbstractMatrix{T}, B::AbstractBSplineBasis, x::AbstractVector,
    [deriv::Derivative = Derivative(0)]; clip_threshold = eps(T))

Fill preallocated collocation matrix.

See collocation_matrix for details.

LinearAlgebra.lu(C::CollocationMatrix, pivot = NoPivot(); check = true)

Returns LU factorisation of collocation matrix.

See also lu!.

LinearAlgebra.lu!(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.

CyclicTridiagonalMatrix{T} <: AbstractMatrix{T}

Represents an almost tridiagonal matrix with non-zero values at the lower-left and upper-right corners.

This kind of matrix appears when working with periodic splines.

Linear systems involving this matrix can be efficiently solved using a combination of the Thomas algorithm (for regular tridiagonal matrices) and of the Sherman–Morrison formula.

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.