Collocation tools

Collocation points

BSplineKit.Collocation.collocation_pointsFunction
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:

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

source

Point selection methods

BSplineKit.Collocation.AvgKnotsType
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).

source
BSplineKit.Collocation.SameAsKnotsType
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:

source

Matrices

BSplineKit.Collocation.CollocationMatrixType
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.

source
BSplineKit.Collocation.collocation_matrixFunction
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.

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

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

Returns LU factorisation of collocation matrix.

See also lu!.

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

source

Internals

LinearAlgebra.ldiv!Function
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.

source