SimpleSplines.jl

CirculantMass(M, n)

The mass matrix of a uniform periodic mesh, represented by the eigenvalues of its circulant structure and a pair of planned real transforms.

On a uniform mesh every basis function is a translate of one cardinal spline, so $\mathbb{M}_{kl}$ depends only on $k - l \bmod N$ and

\[\mathbb{M} = F^{*} \operatorname{diag}(\hat{c}) F , \qquad \hat{c} = \mathcal{F}(\mathbb{M}_{:,1}) ,\]

with $F$ the discrete Fourier transform. A solve is therefore a forward transform, a pointwise division and an inverse transform.

The plans are created once, at construction, and the spectral buffer is preallocated, so a solve allocates nothing beyond its result. \ allocates the result; mass_solve! does not.

The first column is read from the assembled matrix rather than recomputed, and the construction checks that the matrix really is circulant — a silent mismatch here would give wrong answers on every mesh that is uniform by accident rather than by construction.

FactorizedMass(M)

The sparse Cholesky factorisation of the mass matrix, for meshes on which it is not circulant.

$\mathbb{M}$ is symmetric positive definite and banded modulo $N$ — a basis function overlaps only the $2p+1$ others whose supports meet its own — so the factorisation is cheap and the solve is $O(Np)$.

GradedMesh(n, L; amplitude = 0.12)
GradedMesh{T}(n, L; amplitude = 0.12)

The image of a uniform mesh under the fixed smooth map

\[s \mapsto L \left( s + \frac{a}{2\pi} \sin 2 \pi s \right) , \qquad s_i = \frac{i-1}{n} ,\]

with amplitude $= a$. The map does not depend on n, so refining n gives a genuine family of meshes and the rates measured across it are meaningful — which is exactly what RandomMesh does not give.

The map is a bijection of $[0,1]$ with positive derivative for $|a| < 1$, and the cell widths vary by a factor of $(1+a)/(1-a)$ across the domain.

julia> m = GradedMesh(8, 1.0);

julia> issorted(breakpoints(m)) && breakpoints(m)[1] == 0
true

MassOperator{T}

The mass matrix in the form the solves actually want: something that can be applied and, above all, inverted, without forming $\mathbb{M}^{-1}$.

Two representations are provided, and which one is built is decided by the mesh:

• CirculantMass on a UniformMesh, where the basis functions are translates of a single cardinal spline and $\mathbb{M}$ is circulant, hence diagonalised by the discrete Fourier transform. A solve is two transforms and a pointwise division, $O(N \log N)$, with the transforms planned once.
• FactorizedMass otherwise. On a GradedMesh or a RandomMesh the matrix is banded modulo $N$ but not circulant — the basis functions are no longer translates of one another — so there is nothing for a Fourier transform to diagonalise, and a sparse Cholesky factorisation is what is left.

Both answer \, ldiv! and Matrix.

Mesh{T}

A subdivision of the periodic domain $\Omega = [0,L)$ into n cells by the breakpoints

\[0 = y_1 < y_2 < \dots < y_n < L ,\]

which are continued periodically, $y_{i+n} = y_i + L$, so that the last cell is $[y_n, y_1 + L)$. There are n breakpoints and n cells, not n+1: on a torus the right endpoint of the last cell is the left endpoint of the first.

The breakpoints are to be distinguished from the knot vector of PeriodicBSplineBasis, whose entries may repeat.

Three families are provided, and the difference between them matters for testing rather than for use:

• UniformMesh — equally spaced. The assembled matrices are circulant, which makes several quantities vanish identically that do not vanish in general.
• GradedMesh — the image of a uniform mesh under a fixed smooth map. Refining n therefore gives a genuine family of meshes, and rates of convergence measured across it are meaningful.
• RandomMesh — randomly perturbed cell widths. A different mesh for every n, so convergence rates across it mean nothing; its use is to check the properties that must hold exactly on any mesh.

PeriodicBSplineBasis(mesh, p)
PeriodicBSplineBasis{T}(mesh, p)
PeriodicBSplineBasis(n, p; L = 2π)

The periodic B-spline basis of degree p on the Mesh mesh, spanning the spline space $\mathcal{S}^p_n$ of $\mathcal{C}^{p-1}$ piecewise polynomials of degree p on the torus $\Omega = [0,L)$.

julia> b = PeriodicBSplineBasis(UniformMesh(8, 1.0), 3);

julia> nbasis(b), degree(b), order(b)
(8, 3, 4)

julia> sum(b[0.3, j] for j in eachindex(b)) ≈ 1      # partition of unity
true

Dimension

The number of degrees of freedom is $N = n$, the number of cells — not $n + p$, as it would be on a bounded interval. The construction differs from the bounded case only in how the knot vector is closed up: instead of repeating the end knots to clamp the basis at the two ends, the breakpoints are continued periodically, $y_{i+n} = y_i + L$, and the recursion is applied to the resulting bi-infinite sequence. The splines it produces satisfy $\phi_{j+n}^p (x) = \phi_j^p (x - L)$, so only n of them are distinct on $\Omega$, and the periodic basis is obtained by wrapping these onto the torus,

\[\phi_j^p \big\vert_\Omega (x) = \sum_{r \in \mathbb{Z}} \phi_j^p (x + rL) , \qquad 1 \le j \le n ,\]

a finite sum, since $\phi_j^p$ is supported on the p+1 cells $[y_j, y_{j+p+1}]$. The p extra functions of the bounded case are precisely those the clamping introduces at the two ends, and the wrapping identifies them in pairs.

Why this matters

There are no boundary functions at all: every basis function spans p+1 cells, and the basis is $\mathcal{C}^{p-1}$ across the seam $a \equiv b$ as well as inside. Integration by parts over $\Omega$ therefore leaves no boundary terms, which is what makes the discrete brackets assembled from this basis exactly antisymmetric.

On a UniformMesh the basis functions are in addition translates of a single cardinal B-spline, and the mass, stiffness and derivative matrices are circulant.

Requirements

p ≥ 0 and n > p. The second is what makes the wrapping well defined: a basis function spans p+1 cells, so with n ≤ p it would wrap onto itself and the sum above would not terminate.

See also SplineQuadrature for the assembly built on this basis, and evaluate for derivatives of arbitrary order.

PeriodicBSplineDerivative

The type of Derivative(axes(b,1)) * b for a PeriodicBSplineBasis b, equivalently of b'.

A lazy product: it stores the basis and runs the derivative recursion on indexing. For derivatives of order higher than one use evaluate with an explicit d.

RandomMesh(n, L; seed = 1, spread = 0.6)
RandomMesh{T}(n, L; seed = 1, spread = 0.6)

A mesh of n cells whose widths are drawn as $1 + \sigma \, r_i$ with $r_i$ uniform on $[0,1)$ and spread $= \sigma$, then rescaled to sum to L. The stream is seeded, so the mesh is reproducible.

The widths vary by up to a factor of $1 + \sigma$, and a different n gives an unrelated mesh rather than a refinement of this one. Convergence rates measured across a family of RandomMeshes are therefore meaningless; use GradedMesh for those. What this mesh is for is the properties that must hold exactly on any mesh — the antisymmetry of the discrete brackets, for instance, which a uniform mesh would confirm for the wrong reason, its assemblies being circulant.

SplineQuadrature(basis; nq = quadrature_order(degree(basis)), dmax = 3)

The assembly table of a PeriodicBSplineBasis: its basis functions and their derivatives tabulated at the global Gauß-Legendre quadrature points, together with the quadrature weights and the mass matrix.

This is the one data structure every assembly here is built from. With Φ[d+1][i,q] the d-th derivative of $\phi_i$ at the q-th quadrature point and w the weight vector, every matrix of the form $\int_\Omega f(x) \, D^a \phi_k \, D^b \phi_l \, dx$ is a single weighted contraction,

\[A = \Phi_a \, \mathrm{diag}(f \odot w) \, \Phi_b^T ,\]

which is what mixed_matrix and weighted_matrix do. Writing the assembly this way rather than element-by-element is what makes a variable coefficient — the $\int u_h \phi_k \phi_l'$ of the second KdV bracket, say — cost no more than a constant one.

Arguments

• nq: quadrature points per cell. The default integrates degree $3p-1$ exactly; see quadrature_order.
• dmax: the highest derivative order tabulated. Three is what a third-order operator needs after one integration by parts.

julia> q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(16, 2π), 3));

julia> S = derivative_matrix(q);

julia> maximum(abs, S + S') < 1e-14        # antisymmetric on a periodic mesh
true

Storage

Φ is stored densely, N by n * nq. Only the entries inside each basis function's support of p+1 cells are ever nonzero and only those are computed, but keeping the array dense lets the contractions above run as one BLAS call, which is the faster arrangement at the sizes these discretisations are used at.

UniformMesh(n, L)
UniformMesh{T}(n, L)

The uniform mesh of n cells on $[0,L)$, $y_i = (i-1) L / n$.

julia> breakpoints(UniformMesh(4, 1.0))
4-element Vector{Float64}:
 0.0
 0.25
 0.5
 0.75

On a uniform mesh the B-spline basis functions are translates of a single cardinal spline, so the mass, stiffness and derivative matrices assembled from it are circulant. Several identities of the discrete brackets hold on a uniform mesh and nowhere else; a test that means to check a property valid on any mesh should use RandomMesh instead.

nodes(b::PeriodicBSplineBasis)

The Greville abscissae of b, reduced onto $[0,L)$,

\[\xi_j = \frac{1}{p} \sum_{i=1}^{p} x_{j+i} ,\]

the averages of the p interior knots of each basis function.

These are the points a spline basis is naturally interpolated at: $\xi_j$ lies in the support of $\phi_j$ and the collocation matrix $\phi_j(\xi_i)$ is invertible (Schoenberg-Whitney). For p = 1 they are the breakpoints themselves.

order(b::PeriodicBSplineBasis)

The order $k = p + 1$ of the spline basis b.

Note that this is the spline meaning of the word — a B-spline of order k is piecewise of degree k-1, and the knot vector of a basis of N functions has N + k entries. It is not the meaning order carries for the bases of CompactBasisFunctions, where it is the number of basis functions. For a spline basis the two differ: order(b) == degree(b) + 1, while nbasis(b) is the number of cells.

_bspline(kv, k, p, x, d)

The d-th derivative at x of the B-spline of degree p starting at knot index k of the knot vector kv, i.e. the one supported on kv[k] .. kv[k+p+1].

This is the Cox-de Boor recursion written out directly,

\[\phi_j^p = w_j^p \, \phi_j^{p-1} + (1 - w_{j+1}^p) \, \phi_{j+1}^{p-1} , \qquad w_j^p = \frac{x - x_j}{x_{j+p} - x_j} ,\]

with $\phi_j^0$ the indicator of $[x_j, x_{j+1})$, together with the derivative recursion

\[\frac{d}{dx} \phi_j^p = p \left( \frac{\phi_j^{p-1}}{x_{j+p} - x_j} - \frac{\phi_{j+1}^{p-1}}{x_{j+p+1} - x_{j+1}} \right) .\]

Both are evaluated as written rather than through a faster equivalent, so that the assembly built on them doubles as a check of the formulae themselves.

A knot span of zero length — a repeated knot — contributes nothing, which is what makes the recursion well defined at a knot of multiplicity greater than one. The guard is on the span being positive rather than on the numerator, because it is the division that fails.

basis_integrals(q::SplineQuadrature)

The vector $\int_\Omega \phi_i \, dx$.

This is the gradient of the total mass $C_0 = \int_\Omega u \, dx$ with respect to the degrees of freedom, and it equals $\mathbb{M} \mathbf{1}$ because the basis is a partition of unity. It spans the kernel of the first discrete bracket, which is why the mass is a Casimir there.

Assembled once when the SplineQuadrature is built and returned by reference, as mass_matrix is — do not mutate the result.

basis_values(q::SplineQuadrature, d = 0)

The table Φ[i,r] of the d-th derivative of $\phi_i$ at the r-th quadrature point.

breakpoints(m::Mesh)

The n breakpoints of m in $[0,L)$, in increasing order and starting at zero.

The cell boundaries, which repeat the first breakpoint at $L$, are cellbounds.

cellbounds(m::Mesh)

The n+1 cell boundaries [y_1, …, y_n, y_1 + L], i.e. breakpoints with the periodic image of the first breakpoint appended, so that cell k is cellbounds(m)[k] .. cellbounds(m)[k+1].

derivative_matrix(q::SplineQuadrature)

The matrix $S_{kl} = \int_\Omega \phi_k \phi_l' \, dx$.

On a periodic mesh this is already antisymmetric, with no skew-symmetrisation needed: $\int_\Omega \partial_x (\phi_k \phi_l) \, dx = 0$ because the basis is $\mathcal{C}^{p-1}$ across the seam and there are no boundary terms. Its skew part is therefore S itself and not S/2, which is where the factor of one half in the first discrete bracket comes from.

The proviso is that the quadrature integrate that total derivative, of degree $2p-1$, exactly — which means nq ≥ p, one point per cell more than quadrature_order would suggest is needed at low degree. Below it the identity fails outright rather than gracefully: at p = 3 and nq = 2 the defect is $9 \times 10^{-3}$ on a RandomMesh. On a UniformMesh it holds at any nq, the assemblies being circulant, so a check run only there confirms it for the wrong reason.

domainlength(m::Mesh)

The period $L$ of the domain $\Omega = [0,L)$.

evaluate(b::PeriodicBSplineBasis, j, x, d = 0)
evaluate(b::PeriodicBSplineBasis, û::AbstractVector, x, d = 0)

The d-th derivative of the j-th basis function of b at x, or of the spline $u_h = \sum_j \hat{u}_j \phi_j$ with coefficients û.

x may be any real number: it is reduced onto $[0,L)$ first, so the result is the periodic extension. x may also be a vector, in which case a vector is returned.

The explicit derivative order is what the lazy b' products cannot give: the discrete brackets need d up to 3, and stacking Derivative that deep does not compose.

julia> b = PeriodicBSplineBasis(UniformMesh(8, 1.0), 3);

julia> evaluate(b, 1, 0.1) ≈ b[0.1, 1]
true

julia> abs(evaluate(b, 1, 1.1) - evaluate(b, 1, 0.1)) < 1e-14     # periodic
true

knotvector(b::PeriodicBSplineBasis)

The periodically extended knot sequence the recursion runs on.

Five images of the breakpoints, [y .- 2L; y .- L; y; y .+ L; y .+ 2L]. Basis function j is supported on knotvector(b)[2n+j] .. knotvector(b)[2n+j+p+1] before wrapping.

l2_projection!(û, q::SplineQuadrature, f)

In-place l2_projection, writing the coefficients into û.

On a UniformMesh this allocates nothing: the $f \odot w$ product goes into a buffer held by the quadrature, the load vector is formed with mul! straight into û, and the CirculantMass solve is itself allocation-free. On a non-uniform mesh the CHOLMOD solve still allocates a temporary, as mass_solve! notes.

A sample whose element type is wider than the quadrature's — a complex f — gets its own product instead of being narrowed into that buffer, so this method accepts exactly what l2_projection accepts and merely stops being allocation-free there.

The shared buffer is one of the two reasons a SplineQuadrature may not be used from two threads at once; see the warning there.

l2_projection(q::SplineQuadrature, f)

The coefficients of the $L^2$ projection of f onto the spline space, $\hat{u} = \mathbb{M}^{-1} \int_\Omega f \phi_i \, dx$.

f may be a function or a vector of values at quadrature_nodes.

julia> q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(32, 2π), 3));

julia> û = l2_projection(q, sin);

julia> abs(evaluate(basis(q), û, 1.0) - sin(1.0)) < 1e-5
true

mass_factorization(q::SplineQuadrature)

The Cholesky factorization of mass_matrix, for solving with the mass matrix without refactorizing.

mass_matrix(op::MassOperator)

The assembled mass matrix behind the operator.

mass_matrix(q::SplineQuadrature)

The mass matrix $\mathbb{M}_{kl} = \int_\Omega \phi_k \phi_l \, dx$, symmetric positive definite.

Assembled once when the SplineQuadrature is built and returned by reference; see mass_factorization for the Cholesky factor that goes with it.

mass_operator(M, mesh)

Build the MassOperator appropriate to mesh: a CirculantMass on a uniform mesh, a FactorizedMass otherwise.

mass_operator(q::SplineQuadrature)

The MassOperator of the quadrature — a CirculantMass on a uniform mesh, a FactorizedMass otherwise.

mass_solve!(y, op::MassOperator, x)

Solve $\mathbb{M} y = x$ in place. y and x may alias.

Attached to the function rather than to either method, so that the two representations of MassOperator share one piece of documentation and a cross-reference to the name resolves.

meshwidth(m::Mesh)

The largest cell width, the $h$ that convergence rates are measured against.

mixed_matrix(q::SplineQuadrature, a, b)

The matrix $\int_\Omega D^a \phi_k \, D^b \phi_l \, dx$.

julia> q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(16, 2π), 3));

julia> K0 = mixed_matrix(q, 1, 2);      # ∫ φ_k' φ_l''  =  -∫ φ_k φ_l'''

julia> maximum(abs, K0 + mixed_matrix(q, 0, 3)) < 1e-10
true

ncells(m::Mesh)

The number of cells n, which on a periodic mesh equals the number of breakpoints.

quadrature_nodes(q::SplineQuadrature)

The global quadrature points, nq Gauß-Legendre points in each of the n cells, concatenated in cell order.

quadrature_order(p)

The number of Gauß-Legendre points per cell that integrates degree $3p-1$ exactly, $n_q = \lceil 3p/2 \rceil$.

An nq-point Gauß-Legendre rule is exact to degree $2 n_q - 1$, so this is the smallest nq with $2 n_q - 1 \ge 3p - 1$. Degree $3p-1$ is what the consistency of a Galerkin discretisation of a quadratic nonlinearity needs — an identity such as $\int 6 u_h u_{h,x} \phi_i = -3 \int u_h^2 \phi_i'$ pairs three basis functions and one derivative. The mass matrix alone would need only $2p-1$.

What this rule is not needed for is antisymmetry. A bracket assembled in explicitly skew-symmetrised form is antisymmetric to the last bit at any nq and on any mesh; it is accuracy, not structure, that the quadrature buys.

julia> quadrature_order.(1:4)
4-element Vector{Int64}:
 2
 3
 5
 6

quadrature_weights(q::SplineQuadrature)

The global quadrature weights, scaled by the cell widths, so that sum(quadrature_weights(q)) == L.

stiffness_matrix(q::SplineQuadrature)

The stiffness matrix $\mathbb{K}^1_{kl} = \int_\Omega \phi_k' \phi_l' \, dx$, symmetric positive semi-definite with the constants in its kernel.

weighted_matrix(q::SplineQuadrature, f, a, b)

The matrix $\int_\Omega f(x) \, D^a \phi_k \, D^b \phi_l \, dx$.

f may be a function of the coordinate, or a vector already sampled at quadrature_nodes — the second form is what a variable coefficient given as a spline expansion becomes, and it avoids resampling a field that is already in hand.

julia> q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(16, 2π), 3));

julia> A = weighted_matrix(q, sin, 0, 1);   # ∫ sin(x) φ_k φ_l'

julia> size(A)
(16, 16)