Library

PolynomialBasis{T} <: Basis{T}

Supertype of the four bases of this package: a ContinuumArrays Basis of element type T, spanning the polynomials of degree $\le p$ on the reference interval $[0,1]$.

A subtype supplies one internal method, _eval, evaluating one basis function at one point, and the data that method needs. Everything a caller sees follows from it here — basis, nbasis, order, degree, eachindex, axes, the four indexing forms, equality and the derivative product — so that a family's own file holds only what is peculiar to it.

The hierarchy splits by whether a basis carries nodes, into NodalBasis and ModalBasis; see Nodal and modal bases.

NodalBasis{T} <: PolynomialBasis{T}

A PolynomialBasis built from a set of nodes, which it stores in a field x: Lagrange and Chebyshev. It has one basis function per node.

nodes, nnodes and grid are answered from that field. See Nodal and modal bases.

ModalBasis{T} <: PolynomialBasis{T}

A PolynomialBasis without nodes, built from a number of basis functions, which it stores in a field n: Bernstein and Legendre.

Its coefficients are those of an expansion rather than values at points, so nodes, nnodes and grid have no answer and throw. See Nodal and modal bases.

basis(b::PolynomialBasis)

Return the basis functions of b as callables, indexed as b itself is, so that basis(b)[j](x) == b[x,j].

The functions are built on each call. Prefer indexing b directly; this accessor exists for the rare case that the individual functions are needed as values.

Note that ContinuumArrays exports a different basis, which for a basis object returns the basis itself. Loading both packages with using therefore makes the name ambiguous; qualify it, or import the one that is wanted.

nbasis(b::PolynomialBasis)

Return the number of basis functions of b.

This equals nnodes for the nodal bases, where each basis function belongs to one node, and is defined for the modal bases too, where nnodes is not.

nodes(b::NodalBasis)

Return the nodes of b, the points $x_i \in [0,1]$ that its construction is based on.

Defined for the nodal bases Lagrange and Chebyshev. The modal bases Bernstein and Legendre have no nodes and throw an error; see Nodal and modal bases.

See also nnodes for their number and grid for the same vector under the name ContinuumArrays uses.

nnodes(b::NodalBasis)

Return the number of nodes of b, i.e. the length of nodes.

Defined only for the nodal bases; see nodes.

order(b::PolynomialBasis)

Return the order of b, the number of basis functions, hence the number of coefficients an expansion in b has.

order(b) == nbasis(b) and order(b) == degree(b) + 1 for every basis here. The name is the one the ecosystem uses for the accuracy of an approximation: a basis of order $p$ spans the polynomials of degree $\le p-1$ and so reproduces them exactly.

degree(b::PolynomialBasis)

Return the degree of b, the highest polynomial degree it spans.

degree(b) == order(b) - 1 == nbasis(b) - 1 for every basis here.

Lagrange(x)
Lagrange{T}(x)

The Lagrange basis on the nodes x, generally taken in the interval $[0,1]$,

\[\ell_j(x) = \prod_{i \neq j} \frac{x - x_i}{x_j - x_i} ,\]

indexed from 1, unlike the other three bases, which are indexed from 0.

This is the cardinal basis of the nodes: $\ell_j(x_i) = \delta_{ij}$, so the coefficients of an expansion are the values of the function at the nodes, and no linear system is needed to interpolate. The basis also forms a partition of unity, $\sum_j \ell_j(x) = 1$.

julia> l = Lagrange([0.0, 0.5, 1.0]);

julia> all(l[nodes(l)[i], j] ≈ (i == j) for i in 1:3, j in 1:3)   # ℓⱼ(xᵢ) = δᵢⱼ
true

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

The nodes must be distinct and finite, since the denominators are products of their differences; anything that leaves such a product zero or non-finite throws an ArgumentError. That covers a repeated node, the pair 0.0 and -0.0, whose difference is zero although the two are not isequal, and a NaN or Inf node. They need not be sorted, nor confined to $[0,1]$, although the declared domain is $[0,1]$.

Which nodes to use matters: equidistant nodes make high-degree interpolation diverge near the endpoints (the Runge phenomenon), so the Gauß-Legendre and Lobatto-Legendre node sets are provided as LagrangeGauß and LagrangeLobatto; see Choice of nodes.

The element type T is taken from the nodes, and all internal quantities are formed in it, so an arbitrary-precision basis carries its full precision:

julia> setprecision(BigFloat, 256) do
           l = Lagrange(BigFloat[0, 1//4, 1])
           abs(sum(l[BigFloat(1)/3, j] for j in eachindex(l)) - 1) < 1e-70
       end
true

See also Chebyshev for the other nodal basis, and Lagrange basis for the full discussion.

LagrangeGauß(n)

The Lagrange basis on the n Gauß-Legendre nodes of $[0,1]$.

These lie strictly inside the interval, which suits a basis whose expansion is integrated rather than matched at the boundary. Compare LagrangeLobatto, whose nodes include the endpoints.

julia> nodes(LagrangeGauß(2)) ≈ [(1 - 1/sqrt(3)) / 2, (1 + 1/sqrt(3)) / 2]
true

LagrangeLobatto(n)

The Lagrange basis on the n Lobatto-Legendre nodes of $[0,1]$.

These include both endpoints, so an expansion has coefficients that are the boundary values themselves — what a method needs when it has to impose or read off conditions there. Compare LagrangeGauß, whose nodes lie strictly inside.

julia> nodes(LagrangeLobatto(3)) == [0.0, 0.5, 1.0]
true

Chebyshev{kind}(n)
Chebyshev{kind}(T, n)
ChebyshevT(n)     # kind = 1
ChebyshevU(n)     # kind = 2

Chebyshev basis of the first (kind = 1) or second (kind = 2) kind on the interval $[0,1]$, of n functions, i.e. of degree $p = n-1$, indexed from 0.

The basis functions are the Chebyshev polynomials evaluated at 2x-1, and the nodes are the Chebyshev points shifted onto [0..1] and returned in ascending order.

\[T_j(\tilde{x}) = \cos(j \arccos \tilde{x}) , \qquad U_j(\tilde{x}) = \frac{\sin\big((j+1) \arccos \tilde{x}\big)}{\sin(\arccos \tilde{x})} , \qquad \tilde{x} = 2x-1 ,\]

both satisfying the recurrence $\phi_j = 2\tilde{x} \phi_{j-1} - \phi_{j-2}$, which is how they are evaluated; see The reference interval for the shift and its chain-rule factor 2 on derivatives.

Unlike the modal bases, a Chebyshev basis carries nodes: the Chebyshev points of the corresponding kind, from QuadratureRules, shifted onto $[0,1]$ and ascending. Those of the first kind are the roots of $T_n$ and lie strictly inside the interval; those of the second kind are the extrema of $T_{n-1}$ and include both endpoints, which is why kind = 2 requires n ≥ 2.

julia> t = ChebyshevT(3);

julia> t[0.0, 2], t[0.5, 2], t[1.0, 2]      # T₂(x̃) at x̃ = -1, 0, +1
(1.0, -1.0, 1.0)

julia> nodes(ChebyshevU(3))                 # includes both endpoints
3-element Vector{Float64}:
 0.0
 0.5
 1.0

julia> all(x -> x ∈ axes(t, 1), nodes(t))   # the nodes lie in the domain
true

T is the element type of the nodes, and must be able to represent them. This means a floating-point type in general; an integer-like T only works in the special cases where the nodes are exactly representable, such as ChebyshevU(Integer, 2), whose nodes are 0 and 1. Otherwise the constructor throws an InexactError.

See also Chebyshev basis for the full discussion, and Lagrange for the other nodal basis.

Legendre(n)
Legendre(T, n)

The Legendre basis of n functions, i.e. of degree $p = n-1$, on the interval $[0,1]$,

\[L_j(x) = \sqrt{2j+1} \, P_j(2x-1) , \qquad j = 0, \dots, p ,\]

indexed from 0. T is the element type and defaults to Float64.

The Legendre polynomials $P_j$ are defined on $[-1,+1]$, so the argument is shifted by $\tilde{x} = 2x-1$; see The reference interval.

The basis is modal: its coefficients are those of an expansion rather than values at points, so it has no nodes and no grid, and those accessors throw.

The factor $\sqrt{2j+1}$ makes the basis orthonormal on $[0,1]$: since $\int_0^1 P_i(2x-1) P_j(2x-1) \, dx = \delta_{ij} / (2j+1)$,

\[\int_0^1 L_i(x) \, L_j(x) \, dx = \delta_{ij} ,\]

so the mass matrix is the identity and the coefficients of a projection are just the inner products against the basis functions.

julia> l = Legendre(3);

julia> l[0.5, 0], l[0.5, 1], l[0.5, 2]      # the midpoint is x̃ = 0
(1.0, 0.0, -1.118033988749895)

julia> l[1.0, 1] ≈ sqrt(3)                  # L₁(1) = √3 P₁(1) = √3
true

julia> using QuadratureRules

julia> quad = GaussLegendreQuadrature(8);

julia> sum(weights(quad)[k] * l[nodes(quad)[k], 1]^2 for k in eachindex(nodes(quad))) ≈ 1
true

The derivative follows from differentiating Bonnet's recurrence, with the chain-rule factor 2 from the shift, and is obtained as Derivative(axes(l,1)) * l; see Derivatives.

See also Bernstein for the other modal basis, and Legendre basis for the full discussion.

Bernstein(n)
Bernstein(T, n)

The Bernstein basis of n functions, i.e. of degree $p = n-1$, on the interval $[0,1]$,

\[B_{j,p}(x) = \binom{p}{j} \, x^j \, (1-x)^{p-j} , \qquad j = 0, \dots, p ,\]

indexed from 0. T is the element type and defaults to Float64.

The basis is modal: its coefficients are not values at points, so it has no nodes and no grid, and those accessors throw. The coefficients are the control values of a Bézier curve.

On $[0,1]$ the basis functions are non-negative and form a partition of unity, $\sum_j B_{j,p}(x) = 1$, so an expansion is a convex combination of its coefficients and therefore lies in their convex hull. The first and last function interpolate the endpoints, $B_{0,p}(0) = B_{p,p}(1) = 1$, while every other function vanishes at both.

julia> b = Bernstein(3);

julia> b[0.0, 0], b[0.5, 0], b[1.0, 0]     # B₀ interpolates the left endpoint
(1.0, 0.25, 0.0)

julia> b[0.0, 2], b[0.5, 2], b[1.0, 2]     # B₂ interpolates the right endpoint
(0.0, 0.25, 1.0)

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

The derivative is again a Bernstein expansion, of one degree less,

\[B_{j,p}'(x) = p \, \big( B_{j-1,p-1}(x) - B_{j,p-1}(x) \big) ,\]

and is obtained as Derivative(axes(b,1)) * b; see Derivatives.

See also Legendre for the other modal basis, and Bernstein basis for the full discussion.

PolynomialBasisDerivative

The type of Derivative(axes(b,1)) * b for any PolynomialBasis b, equivalently of b', and the supertype of the four per-family derivative types.

A lazy product: it stores the basis and evaluates the derivative on indexing, in the same four forms as the basis itself. See Derivatives.

BernsteinDerivative

The type of Derivative(axes(b,1)) * b for a Bernstein basis b, equivalently of b'. See PolynomialBasisDerivative and Derivatives.

ChebyshevDerivative

The type of Derivative(axes(b,1)) * b for a Chebyshev basis b of either kind, equivalently of b'. The two kinds need different formulas, and ChebyshevTDerivative and ChebyshevUDerivative name them apart. See PolynomialBasisDerivative and Derivatives.

ChebyshevTDerivative

ChebyshevDerivative narrowed to the first kind, evaluated through $T_j' = j \, U_{j-1}$.

The alias names the type. Evaluation dispatches on ChebyshevDerivative, which picks the formula from the kind of the basis it was built from.

ChebyshevUDerivative

ChebyshevDerivative narrowed to the second kind, evaluated through the differentiated recurrence rather than the closed form, which is singular at the endpoints.

The alias names the type. Evaluation dispatches on ChebyshevDerivative, which picks the formula from the kind of the basis it was built from.

LagrangeDerivative

The type of Derivative(axes(l,1)) * l for a Lagrange basis l, equivalently of l'. See PolynomialBasisDerivative and Derivatives.

LegendreDerivative

The type of Derivative(axes(l,1)) * l for a Legendre basis l, equivalently of l'. See PolynomialBasisDerivative and Derivatives.

vandermonde_matrix(x)

The Vandermonde matrix of the nodes x, $V_{ij} = x_i^{j-1}$.

V is the matrix of the monomial basis evaluated at the nodes, so V * c are the values at the nodes of the polynomial with monomial coefficients c, and V \ y are the monomial coefficients of the polynomial interpolating the values y.

julia> vandermonde_matrix([0.0, 0.5, 1.0])
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 1.0  0.5  0.25
 1.0  1.0  1.0

See also vandermonde_matrix_inverse and The Vandermonde matrix.

vandermonde_matrix_inverse(x)

The inverse of vandermonde_matrix, computed in closed form.

V is factored as $V^{-1} = U L$ with L lower and U upper triangular, whose entries are known explicitly from the nodes: L holds the reciprocals of the products of node differences that also appear in the Lagrange cardinal functions, and U is built by a two-term recurrence. This avoids a general factorisation of a matrix that is notoriously ill-conditioned, though the result is still limited by that conditioning.

V⁻¹ * y are the monomial coefficients of the polynomial interpolating the values y at the nodes x.

julia> vandermonde_matrix_inverse([0.0, 0.5, 1.0])
3×3 Matrix{Float64}:
  1.0   0.0   0.0
 -3.0   4.0  -1.0
  2.0  -4.0   2.0

The nodes must be distinct, as for Lagrange; repeated nodes divide by zero.