This is a very simple package for building blocks of spectral methods. Its intended audience is users who are familiar with the theory and practice of these methods, and prefer to assemble their code from modular building blocks. If you need an introduction, a book like Boyd (2001): Chebyshev and Fourier spectral methods is a good place to start.

The package is optimized for solving functional equations, as usually encountered in economics when solving discrete and continuous-time problems. It uses static arrays extensively to avoid allocation and unroll some loops. Key functionality includes evaluating a set of basis functions, their linear combination at arbitrary points in a fast manner, for use in threaded code. These should work seamlessly with automatic differentiation frameworks, but also has its own primitives for obtaining derivatives of basis functions.


In this package,

  1. A basis is a finite family of functions for approximating other functions. The dimension of a basis tells you how many functions are in there, while domain can be used to query its domain.

  2. A grid is vector of suggested gridpoints for evaluating the function to be approximated that has useful theoretical properties. You can contruct a collocation_matrix using this grid (or any other set of points). Grids are associated with bases at the time of their construction: a basis with the same set of functions can have different grids.

  3. basis_at returns an iterator for evaluating basis functions at an arbitrary point inside their domain. This iterator is meant to be heavily optimized and non-allocating. linear_combination is a convenience wrapper for obtaining a linear combination of basis functions at a given point.

A basis is constructed using

  1. a family on a fixed domain, eg Chebyshev,

  2. a grid specification like InteriorGrid,

  3. a number of parameters that determine the number of basis functions.

A set of coordinates for a particular basis can be augmented for a wider basis with augment_coefficients.

Currenly, all bases have the domain $[-1,1]$ or $[-1,1]^n$. Facilities are provided for coordinatewise transformations to other domains.


Univariate family on [-1,1]

julia> using SpectralKit
julia> basis = Chebyshev(EndpointGrid(), 5) # 5 Chebyshev polynomialsChebyshev polynomials (1st kind), EndpointGrid(), dimension: 5
julia> is_function_basis(basis) # ie we support the interface belowtrue
julia> dimension(basis) # number of basis functions5
julia> domain(basis) # domain(-1, 1)
julia> grid(basis) # Gauss-Lobatto gridSpectralKit.ChebyshevGridIterator{Float64, Chebyshev{EndpointGrid}}(Chebyshev polynomials (1st kind), EndpointGrid(), dimension: 5)
julia> collect(basis_at(basis, 0.41)) # iterator for basis functions at 0.415-element Vector{Float64}: 1.0 0.41 -0.6638000000000001 -0.9543159999999999 -0.11873911999999986
julia> collect(basis_at(basis, derivatives(0.41))) # values and 1st derivatives5-element Vector{SpectralKit.Derivatives{0, 2, Float64}}: SpectralKit.Derivatives{0, 2, Float64}((1.0, 0.0)) SpectralKit.Derivatives{0, 2, Float64}((0.41, 1.0)) SpectralKit.Derivatives{0, 2, Float64}((-0.6638000000000001, 1.64)) SpectralKit.Derivatives{0, 2, Float64}((-0.9543159999999999, -0.9828000000000003)) SpectralKit.Derivatives{0, 2, Float64}((-0.11873911999999986, -4.354528))
julia> θ = [1, 0.5, 0.2, 0.3, 0.001] # a vector of coefficients5-element Vector{Float64}: 1.0 0.5 0.2 0.3 0.001
julia> linear_combination(basis, θ, 0.41) # combination at some value0.7858264608800001
julia> linear_combination(basis, θ)(0.41) # also as a callable0.7858264608800001
julia> basis2 = Chebyshev(EndpointGrid(), 8) # 8 Chebyshev polynomialsChebyshev polynomials (1st kind), EndpointGrid(), dimension: 8
julia> is_subset_basis(basis, basis2) # we could augment θ …true
julia> augment_coefficients(basis, basis2, θ) # … so let's do it8-element Vector{Float64}: 1.0 0.5 0.2 0.3 0.001 0.0 0.0 0.0

Smolyak approximation on a transformed domain

julia> using SpectralKit, StaticArrays
julia> function f2(x) # bivariate function we approximate x1, x2 = x # takes vectors exp(x1) + exp(-abs2(x2)) endf2 (generic function with 1 method)
julia> ct = coordinate_transformations(BoundedLinear(-1, 2.0), SemiInfRational(-3.0, 3.0))coordinate transformations (-1.0,2.0) ↔ (-1, 1) [linear transformation] (-3.0,∞) ↔ (-1, 1) [rational transformation with scale 3.0]
julia> basis = smolyak_basis(Chebyshev, InteriorGrid2(), SmolyakParameters(3), 2)Sparse multivariate basis on ℝ² Smolyak indexing, ∑bᵢ ≤ 3, all bᵢ ≤ 3, dimension 49 using Chebyshev polynomials (1st kind), InteriorGrid2(), dimension: 15
julia> x = grid(basis)SpectralKit.SmolyakGridIterator{Tuple{Float64, Float64}, SpectralKit.SmolyakIndices{2, 15, 3, 3, 4}, StaticArraysCore.SVector{15, Float64}}(Smolyak indexing, ∑bᵢ ≤ 3, all bᵢ ≤ 3, dimension 49, [0.0, -0.7071067811865476, 0.7071067811865476, -0.9238795325112867, -0.3826834323650898, 0.3826834323650898, 0.9238795325112867, -0.9807852804032304, -0.8314696123025452, -0.5555702330196022, -0.19509032201612828, 0.19509032201612828, 0.5555702330196022, 0.8314696123025452, 0.9807852804032304])
julia> θ = collocation_matrix(basis) \ f2.(from_pm1.(ct, x)) # find the coefficients49-element Vector{Float64}: 2.7793220217324697 3.236988643057215 1.11399024245779 0.2663479965031001 0.048598256445250274 0.00715729546175104 0.0008829533666790555 9.366852829757329e-5 8.71376891398461e-6 7.21659851465705e-7 ⋮ -8.652321510390804e-17 0.030107750183198955 0.048411759340962354 -0.03222472955157697 -0.16530952363028042 0.0163966331843244 0.012259562724398284 -0.015910872405879215 -0.1408521904743876
julia> z = (0.5, 0.7) # evaluate at this point(0.5, 0.7)
julia> isapprox(f2(z), linear_combination(basis, θ, to_pm1(ct, z)), rtol = 0.005)true

Constructing bases

Grid specifications

struct EndpointGrid <: SpectralKit.AbstractGrid

Grid that includes endpoints (eg Gauss-Lobatto).


For small dimensions may fall back to a grid that does not contain endpoints.

struct InteriorGrid2 <: SpectralKit.AbstractGrid

Grid with interior points that results in smaller grids than InteriorGrid when nested. Equivalent to an EndpointGrid with endpoints dropped.


Univariate and multivariate transformations

Bases are defined on the domain $[-1, 1]$ or $[-1, 1]^n$. Transformations map other uni- and multivariate sets into these domains.


to_pm1(transformation, x)

Transform x to $[-1, 1]$ using transformation.

Supports partial application as to_pm1(transformation).

!!! FIXME document, especially differentiability requirements at infinite endpoints


from_pm1(transformation, x)

Transform x from $[-1, 1]$ using transformation.

Supports partial application as from_pm1(transformation).

!!! FIXME document, especially differentiability requirements at infinite endpoints

struct BoundedLinear{T<:Real} <: SpectralKit.UnivariateTransformation

Transform x ∈ (-1,1) to y ∈ (a, b), using $y = x ⋅ s + m$.

m and s are calculated and checked by the constructor; a < b is enforced.

InfRational(A, L)

Chebyshev polynomials transformed to the domain (-Inf, Inf) using $y = A + L ⋅ x / √(1 - x^2)$, with L > 0.

Example mappings

  • $0 ↦ A$
  • $±0.5 ↦ A ± L / √3$
SemiInfRational(A, L)

[-1,1] transformed to the domain [A, Inf) (when L > 0) or (-Inf,A] (when L < 0) using $y = A + L ⋅ (1 + x) / (1 - x)$.

When used with Chebyshev polynomials, also known as a “rational Chebyshev” basis.

Example mappings

  • $-1/2 ↦ A + L / 3$
  • $0 ↦ A + L$
  • $1/2 ↦ A + 3 ⋅ L$

Wrapper for coordinate-wise transformations.

julia> using StaticArrays

julia> ct = coordinate_transformations(BoundedLinear(0, 2), SemiInfRational(2, 3))
coordinate transformations
  (0.0,2.0) ↔ (-1, 1) [linear transformation]
  (2,∞) ↔ (-1, 1) [rational transformation with scale 3]

julia> x = from_pm1(ct, (0.4, 0.5))
(1.4, 11.0)

julia> y = to_pm1(ct, x)
(0.3999999999999999, 0.5)

Univariate bases

Currently, only Chebyshev polynomials are implemented. Univariate bases operate on real numbers.

struct Chebyshev{K<:SpectralKit.AbstractGrid} <: SpectralKit.FunctionBasis

The first N Chebyhev polynomials of the first kind, defined on [-1,1].


Multivariate bases

Multivariate bases operate on tuples or vectors (StaticArrays.SVector is preferred for performance, but all <:AbstractVector types should work).

SmolyakParameters(B, M)

Parameters for Smolyak grids that are independent of the dimension of the domain.

Polynomials are organized into blocks of 1, 2, 2, 4, 8, 16, … polynomials (and corresponding gridpoints), indexed with a block index b that starts at 0. B ≥ ∑ bᵢ and 0 ≤ bᵢ ≤ M constrain the number of blocks along each dimension i.

M > B is not an error, but will be normalized to M = B with a warning.


Create a sparse Smolyak basis.


  • univariate_family: should be a callable that takes a grid_kind and a dimension parameter, eg Chebyshev.

  • grid_kind: the grid kind, eg InteriorGrid() etc.

  • smolyak_parameters: the Smolyak grid specification parameters, see SmolyakParameters.

  • N: the dimension. wrapped in a Val for type stability, a convenience constructor also takes integers.


julia> basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(3), 2)
Sparse multivariate basis on ℝ²
  Smolyak indexing, ∑bᵢ ≤ 3, all bᵢ ≤ 3, dimension 81
  using Chebyshev polynomials (1st kind), InteriorGrid(), dimension: 27

julia> dimension(basis)

julia> domain(basis)
((-1, 1), (-1, 1))


Grids nest: increasing arguments of SmolyakParameters result in a refined grid that contains points of the cruder grid.


Using bases




The domain of a function basis. A tuple of numbers (of arbitrary type, but usually Float64), or a tuple of domains by coordinate.




basis_at(basis, x)

Return an iterable with known element type and length (Base.HasEltype(), Base.HasLength()) of basis functions in basis evaluated at x.

Univariate bases operate on real numbers, while for multivariate bases, Tuples or StaticArrays.SVector are preferred for performance, though all <:AbstractVector types should work.

Methods are type stable.


Consequences are undefined when evaluating outside the domain.

linear_combination(basis, θ, x)

Evaluate the linear combination of $∑ θₖ⋅fₖ(x)$ of function basis $f₁, …$ at x, for the given order.

The length of θ should equal dimension(θ).

linear_combination(basis, θ)

Return a callable that calculates linear_combination(basis, θ, x) when called with x.


Grids and collocation


grid([T], basis)

Return an iterator for the grid recommended for collocation, with dimension(basis) elements.

T for the element type of grid coordinates, and defaults to Float64. Methods are type stable.

collocation_matrix(basis, x)

Convenience function to obtain a “collocation matrix” at points x, which is assumed to have a concrete eltype. The default is x = grid(basis), specialized methods may exist for this when it makes sense.

The collocation matrix may not be an AbstractMatrix, all it needs to support is C \ y for compatible vectors y = f.(x).

Methods are type stable.


Augment coordinates for a wider basis


augment_coefficients(basis1, basis2, θ1)

Return a set of coefficients θ2 for basis2 such that

linear_combination(basis1, θ1, x) == linear_combination(basis2, θ2, x)

for any x in the domain. In practice this means padding with zeros.

Throw a ArgumentError if the bases are incompatible with each other or x, or this is not possible. Methods may not be defined for incompatible bases, compatibility between bases can be checked with is_subset_basis.




API for derivatives is still experimental and subject to change.

If derivatives along a coordinate are needed, use derivatives. For multiple coordinates, the result will be nested in the order of increasing tags. When unspecified, tags are assigned automatically from left to right.

derivatives(::Val(I) = Val(0), x, ::Val(N) = Val(1))

Obtain N derivatives (and the function value) at a scalar x. The ith derivative can be accessed with [i] from results, with [0] for the function value.

I is an integer “tag” for determining the nesting order. Lower I always end up outside when nested. When the defaults are used in multiple coordinates, increasing numbers replace zeros from left to right, starting after the highest explicitly assigned tag.

Consequently, for most applications, you only need to specify tags if you want a different nesting than left-to-right.

Univariate example

julia> basis = Chebyshev(InteriorGrid(), 3)
Chebyshev polynomials (1st kind), InteriorGrid(), dimension: 3

julia> C = collect(basis_at(basis, derivatives(0.1)))
3-element Vector{SpectralKit.Derivatives{0, 2, Float64}}:
 SpectralKit.Derivatives{0, 2, Float64}((1.0, 0.0))
 SpectralKit.Derivatives{0, 2, Float64}((0.1, 1.0))
 SpectralKit.Derivatives{0, 2, Float64}((-0.98, 0.4))

julia> C[1][1]                         # 1st derivative of the linear term is 1

Multivariate example

julia> basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(2), 2)
Sparse multivariate basis on ℝ²
  Smolyak indexing, ∑bᵢ ≤ 2, all bᵢ ≤ 2, dimension 21
  using Chebyshev polynomials (1st kind), InteriorGrid(), dimension: 9

julia> C = collect(basis_at(basis, (derivatives(0.1), derivatives(0.2, Val(2)))));

julia> C[14][1][2]                  # ∂/∂x₁ ∂/∂x₂² of the 14th basis function at x


This section of the documentation is probably only relevant to contributors and others who want to understand the internals.

Simplified API for adding custom transformations

abstract type UnivariateTransformation

An abstract type for univariate transformations. Transformations are not required to be subtypes, this just documents the interface they need to support:


Abstract type used for code organization, not exported.


Grid internals

gridpoint(_, basis, i)

Return a gridpoint for collocation, with 1 ≤ i ≤ dimension(basis).

T is used as a hint for the element type of grid coordinates, and defaults to Float64. The actual type can be broadened as required. Methods are type stable.


Not all grids have this method defined, especially if it is impractical. See grid, which is part of the API, this function isn't.