`BSplineKit.BandedTensors.BandedTensor3D`

— Type`BandedTensor3D{T,b}`

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

**Storage**

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)
end
```

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

.

`BSplineKit.BandedTensors.SubMatrix`

— Type`SubMatrix{T} <: AbstractMatrix{T}`

Represents the submatrix `A[:, :, k]`

of a `BandedTensor3D`

`A`

.

Wraps the `SMatrix`

holding the submatrix.

`BSplineKit.BandedTensors.band_indices`

— Method`band_indices(A::BandedTensor3D, k)`

Return the range of indices `a:b`

for subarray `A[:, :, k]`

where values may be non-zero.

`BSplineKit.BandedTensors.bandshift`

— Method`bandshift(A::BandedTensor3D) -> (δi, δj, δk)`

Return tuple with band shifts along each dimension.

`BSplineKit.BandedTensors.muladd!`

— Method`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`

.

`BandedMatrices.bandwidth`

— Method`bandwidth(A::BandedTensor3D)`

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`

.

`Base.setindex!`

— Method`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`

.

`LinearAlgebra.dot`

— Method`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`

.

`BSplineKit.Recombinations`

— Module`BSplineKit.Recombinations.NoUniqueSolutionError`

— Type`NoUniqueSolutionError <: Exception`

Exception thrown when solving linear system using `RecombineMatrix`

, when the system has no unique solution.

`BSplineKit.Recombinations.RecombineMatrix`

— Type`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.

`BSplineKit.Recombinations.RecombinedBSplineBasis`

— Type`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 `Derivative`

s 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 `Derivative`

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

**Examples**

```
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:

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}$.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))`

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

**Examples**

```
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},)
```

`BSplineKit.Recombinations.constraints`

— Method```
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.

`BSplineKit.Recombinations.num_constraints`

— Method```
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)`

.

`BSplineKit.Recombinations.num_recombined`

— Method```
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.

`BSplineKit.Recombinations.nzrows`

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

Returns the range of row indices `i`

such that `A[i, col]`

is non-zero.

`BSplineKit.Recombinations.parent_coefficients`

— Method`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.

`BSplineKit.Recombinations.recombination_matrix`

— Method`recombination_matrix(R::AbstractBSplineBasis)`

Get `RecombineMatrix`

associated to the recombined basis.

For non-recombined bases such as `BSplineBasis`

, this returns the identity matrix (`LinearAlgebra.I`

).

`Base.length`

— Method`length(R::RecombinedBSplineBasis)`

Returns the number of functions in the recombined basis.

`Base.parent`

— Method`parent(R::RecombinedBSplineBasis)`

Get original B-spline basis.

`BSplineKit.Galerkin.galerkin_matrix`

— Function```
galerkin_matrix(
[f::Function],
B::AbstractBSplineBasis,
[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.

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.

`BSplineKit.Galerkin.galerkin_matrix!`

— Function```
galerkin_matrix!(
[f::Function], A::AbstractMatrix, B::AbstractBSplineBasis, [deriv = (Derivative(0), Derivative(0))];
[quadrature],
)
```

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.

`BSplineKit.Galerkin.galerkin_projection!`

— Function```
galerkin_projection!(
f, φ::AbstractVector, B::AbstractBSplineBasis, [deriv = Derivative(0)],
)
```

Compute Galerkin projection $φ_i = ⟨ b_i, f ⟩$.

See `galerkin_projection`

for details.

`BSplineKit.Galerkin.galerkin_projection`

— Method```
galerkin_projection(
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.

`BSplineKit.Galerkin.galerkin_tensor`

— Function```
galerkin_tensor(
B::AbstractBSplineBasis,
(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.

`BSplineKit.Galerkin.galerkin_tensor!`

— Function```
galerkin_tensor!(
A::BandedTensor3D,
B::AbstractBSplineBasis,
(D₁::Derivative, D₂::Derivative, D₃::Derivative),
)
```

Compute 3D Galerkin tensor in-place.

See `galerkin_tensor`

for details.

`BSplineKit.SplineExtrapolations.AbstractExtrapolationMethod`

— Type`AbstractExtrapolationMethod`

Abstract type representing an extrapolation method.

`BSplineKit.SplineExtrapolations.Flat`

— Type`Flat <: AbstractExtrapolationMethod`

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

`BSplineKit.SplineExtrapolations.Linear`

— Type`Linear <: AbstractExtrapolationMethod`

Represents a linear extrapolation: splines values extend linearly beyond the left and right boundaries.

`BSplineKit.SplineExtrapolations.Smooth`

— Type`Smooth <: AbstractExtrapolationMethod`

Represents a smooth extrapolation: derivatives up to order $k - 2$ are continuous at the boundaries.

`BSplineKit.SplineExtrapolations.SplineExtrapolation`

— Type`SplineExtrapolation`

Represents a spline which can be evaluated outside of its limits according to a given extrapolation method.

`BSplineKit.SplineExtrapolations.extrapolate`

— Method`extrapolate(S::Union{Spline, SplineWrapper}, method::AbstractExtrapolationMethod)`

Construct a `SplineExtrapolation`

from the given spline `S`

(which can also be the result of an interpolation).

`BSplineKit.SplineInterpolations.SplineInterpolation`

— Type`SplineInterpolation`

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.

`BSplineKit.SplineInterpolations.interpolate`

— Function`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 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 (`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).

**Examples**

```
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)
-1.0
julia> (Derivative(1) * S)(-1)
-0.01663433622896893
julia> (Derivative(2) * S)(-1)
10.527273287554928
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)
-1.0
julia> (Derivative(1) * Snat)(-1)
0.2872618670889516
julia> (Derivative(2) * Snat)(-1)
-3.33066907387547e-14
```

**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, 1.51788e-17, -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.9971071640321145, -0.9996420091470221)
julia> x = 0.998; cospi(x), Sper(x), Snat(x), S(x)
(-0.9999802608561371, -0.9999801044078943, -0.9994253145274461, -1.0000122303614758)
```

`BSplineKit.SplineInterpolations.interpolate!`

— Method`interpolate!(I::SplineInterpolation, y::AbstractVector)`

Update spline interpolation with new data.

This function allows to reuse a `SplineInterpolation`

returned by a previous call to `interpolate`

, using new data on the same locations `x`

.

See `interpolate`

for details.

`BSplineKit.SplineApproximations`

— Module`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.

`BSplineKit.SplineApproximations.AbstractApproxMethod`

— Type`AbstractApproxMethod`

Abstract type describing a type of approximation method.

`BSplineKit.SplineApproximations.ApproxByInterpolation`

— Type`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`

.

`BSplineKit.SplineApproximations.MinimiseL2Error`

— Type`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$.

`BSplineKit.SplineApproximations.VariationDiminishing`

— Type`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`

).

`BSplineKit.SplineApproximations.approximate!`

— Method`approximate!(f, A::SplineApproximation)`

Approximate function `f`

reusing a previous Spline approximation in a given basis.

See `approximate`

for details.

`BSplineKit.SplineApproximations.approximate`

— Method`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.

**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)
```

`BSplineKit.Splines.PeriodicVector`

— Type`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]`

.

`PeriodicVector(cs::AbstractVector)`

Wraps coefficient vector `cs`

such that it can be indexed in a periodic manner.

`BSplineKit.Splines.Spline`

— Type`Spline{T} <: Function`

Represents a spline function.

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

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

**Examples**

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

`(S::Spline)(x)`

Evaluate spline at coordinate `x`

.

`BSplineKit.Splines.SplineWrapper`

— Type`SplineWrapper{S <: Spline} <: Function`

Abstract type representing a type that wraps a `Spline`

.

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

`BSplineKit.BSplines.basis`

— Method`basis(S::Spline) -> AbstractBSplineBasis`

Returns the associated B-spline basis.

`BSplineKit.Splines.coefficients`

— Method`coefficients(S::Spline)`

Get B-spline coefficients of the spline.

`BSplineKit.Splines.integral`

— Method`integral(S::Spline)`

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

.

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`

).

`BSplineKit.Splines.spline`

— Method`spline(w::SplineWrapper) -> Spline`

Returns the `Spline`

wrapped by the object.

`Base.:*`

— Method`*(op::Derivative, S::Spline) -> Spline`

Returns `N`

-th derivative of spline `S`

as a new spline.

See also `diff`

.

**Examples**

```
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
true
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]
```

`Base.diff`

— Function`diff(S::Spline, [op::Derivative = Derivative(1)]) -> Spline`

Same as `op * S`

.

Returns `N`

-th derivative of spline `S`

as a new spline.

`Base.eltype`

— Method```
eltype(::Type{<:Spline})
eltype(S::Spline)
```

Returns type of element returned when evaluating the `Spline`

.

`Base.length`

— Method`length(S::Spline)`

Returns the number of coefficients in the spline.

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

.

`BSplineKit.BSplines`

— Module`BSplines`

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

`BSplineKit.BSplines.AbstractBSplineBasis`

— Type`AbstractBSplineBasis{k,T}`

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`

.

```
(B::AbstractBSplineBasis)(
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.

**Examples**

```
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)
1.0
julia> bs[1] - B[i](0.42)
0.0
julia> bs[2] - B[i - 1](0.42)
-5.551115123125783e-17
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))
```

`BSplineKit.BSplines.AugmentedKnots`

— Type`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.

`BSplineKit.BSplines.BSplineBasis`

— Type`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.

**Examples**

```
julia> breaks = range(-1, 1; length = 21)
-1.0:0.1:1.0
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):
0.0
0.0
0.0
0.1
0.2
0.6
1.0
1.0
1.0
```

`BSplineKit.BSplines.BSplineOrder`

— Type`BSplineOrder(k::Integer)`

Specifies the B-spline order `k`

.

`BSplineKit.BSplines.BasisFunction`

— Type`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.

`BSplineKit.BSplines.PeriodicBSplineBasis`

— Type`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]`

.

**Examples**

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

```
julia> ts = range(-1, 1; length = 11)
-1.0:0.2:1.0
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)
2.0
julia> length(B)
10
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}}:
-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
```

`BSplineKit.BSplines.PeriodicKnots`

— Type`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`

.

`BSplineKit.BSplines.augment_knots!`

— Method`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.

`BSplineKit.BSplines.basis_to_array_index`

— Method`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).

`BSplineKit.BSplines.boundaries`

— Function`boundaries(B::AbstractBSplineBasis)`

Returns `(xmin, xmax)`

tuple with the boundaries of the domain supported by the basis.

`BSplineKit.BSplines.common_support`

— Method`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.

`BSplineKit.BSplines.evaluate`

— Function```
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!`

.

`BSplineKit.BSplines.evaluate!`

— Method```
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`

.

`BSplineKit.BSplines.evaluate_all`

— Function```
evaluate_all(
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.

**Examples**

See `AbstractBSplineBasis`

for some examples using the alternative evaluation syntax `B(x, [op], [T]; [ileft])`

, which calls this function.

`BSplineKit.BSplines.find_knot_interval`

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

`BSplineKit.BSplines.has_parent_basis`

— Function```
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.

`BSplineKit.BSplines.isindomain`

— Method`isindomain(B::AbstractBSplineBasis, x::Real)`

Check if the coordinate `x`

is within the boundaries of the domain.

`BSplineKit.BSplines.knots`

— Function```
knots(B::AbstractBSplineBasis)
knots(g::Spline)
```

Returns the knots of the B-spline basis.

`BSplineKit.BSplines.multiplicity`

— Method`multiplicity(knots, i)`

Determine multiplicity of knot `knots[i]`

.

`BSplineKit.BSplines.nonzero_in_segment`

— Method`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.

`BSplineKit.BSplines.order`

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

Returns order of B-splines as an integer.

`BSplineKit.BSplines.static_length`

— Function```
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.

`BSplineKit.BSplines.support`

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

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

`BSplineKit.BSplines.support`

— Method`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)$.

`Base.getindex`

— Method`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.

**Examples**

```
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)
0.16666666666666666
julia> B[6, Float32](-0.5)
0.16666667f0
julia> B[6](-0.5, Derivative(1))
-5.000000000000001
```

`Base.length`

— Method`length(g::BSplineBasis) -> Int`

Returns the number of B-splines composing a spline.

`BSplineKit.Collocation.AvgKnots`

— Type`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).

`BSplineKit.Collocation.CollocationMatrix`

— Type`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`

.

`BSplineKit.Collocation.CyclicTridiagonalMatrix`

— Type`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.

`BSplineKit.Collocation.SameAsKnots`

— Type`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`

).

`BSplineKit.Collocation.SelectionMethod`

— Type`SelectionMethod`

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

`BSplineKit.Collocation.collocation_matrix!`

— Method```
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.

`BSplineKit.Collocation.collocation_matrix`

— Method```
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.

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

.

`BSplineKit.Collocation.collocation_points`

— Function```
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!`

.

`BSplineKit.Collocation.collocation_points!`

— Function```
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.

`LinearAlgebra.ldiv!`

— Method`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.

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

`LinearAlgebra.lu`

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

Returns LU factorisation of collocation matrix.

See also `lu!`

.

`BSplineKit.DifferentialOps`

— Module`DifferentialOps`

Module defining types describing differential operators and compositions thereof.

`BSplineKit.DifferentialOps.AbstractDifferentialOp`

— Type`AbstractDifferentialOp`

Represents a general differential operator.

`BSplineKit.DifferentialOps.AbstractNormalDirection`

— Type`AbstractNormalDirection`

Represents the normal direction on a given domain boundary.

`BSplineKit.DifferentialOps.Derivative`

— Type`Derivative{n} <: AbstractDifferentialOp`

Specifies the `n`

-th derivative of a function.

`BSplineKit.DifferentialOps.DerivativeUnitRange`

— Type`DerivativeUnitRange{m, n} <: AbstractDifferentialOp`

Specifies a range of derivatives.

**Examples**

Two ways of constructing derivative ranges:

```
julia> Derivative(2):Derivative(4)
Derivative(2:4)
julia> Derivative(2:4)
Derivative(2:4)
julia> Tuple(Derivative(2:4))
(D{2}, D{3}, D{4})
```

`BSplineKit.DifferentialOps.DifferentialOpSum`

— Type`DifferentialOpSum <: AbstractDifferentialOp`

Sum of two differential operators.

`BSplineKit.DifferentialOps.LeftNormal`

— Type`LeftNormal <: AbstractNormalDirection`

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

The left normal direction goes opposite to the coordinate axis.

`BSplineKit.DifferentialOps.RightNormal`

— Type`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.

`BSplineKit.DifferentialOps.ScaledDerivative`

— Type`ScaledDerivative{n} <: AbstractDifferentialOp`

`n`

-th derivative of a function scaled by a constant coefficient.

`BSplineKit.DifferentialOps.max_order`

— Function```
max_order(op::AbstractDifferentialOp)
max_order(ops...)
```

Get maximum derivative order of one or more differential operators.

`LinearAlgebra.dot`

— Method`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.

`BSplineKit.BoundaryConditions`

— Module`BSplineKit.BoundaryConditions`

Contains some boundary condition definitions.

`BSplineKit.BoundaryConditions.Natural`

— Type`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.

`BSplineKit.BoundaryConditions.Periodic`

— Type`Periodic <: BoundaryCondition`

Represents periodic boundary conditions with a given period `L`

.

`Periodic(L::Real)`

Constructs periodic boundary conditions with period `L`

.

`BSplineKit.BoundaryConditions.period`

— Method```
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.