SpectralKit

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, potentially reusing calculations. If you need an introduction, a book like Boyd (2001): Chebyshev and Fourier spectral methods is a good place to start.

This package was designed primarily for solving functional equations, as usually encountered in economics when solving discrete and continuous-time problems. Key features include

evaluation of univariate and multivatiate basis functions, including Smolyak combinations,
transformed to the relevant domains of interest, eg $[a,b] × [0,∞)$,
(partial) derivatives, with correct limits at endpoints,
allocation-free, thread safe linear combinations for the above with a given set of coefficients,
using static arrays extensively to avoid allocation and unroll some loops.

While there is some functionality in this package to fit approximations to existing functions, it does not use optimized algorithms (DCT) for that, as it was optimized for mapping a set of coefficients to residuals of functional equations at gridpoints.

Also, while the package should interoperate seamlessly with most AD frameworks, only the derivative API (explained below) is guaranteed to have correct derivatives of limits near infinity.

Concepts

In this package,

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

Univariate basis example

We construct a basis of Chebyshev polynomials on $[0, 4]$. This requires a transformation since their canonical domain is $[-1,1]$. Other transformations include SemiInfRational for $[A, \infty]$ intervals and InfRational for $[-\infty, \infty]$ intervals.

We display the domian and the dimension (number of basis functions).

julia> using SpectralKit
julia> basis = Chebyshev(InteriorGrid(), 5) ∘ BoundedLinear(0, 4)SpectralKit.TransformedBasis{Chebyshev{InteriorGrid}, BoundedLinear{Float64}}(Chebyshev polynomials (1st kind), InteriorGrid(), dimension: 5, (0.0,4.0) ↔ domain [linear transformation])
julia> domain(basis)[0.0,4.0]
julia> dimension(basis)5

We have chosen an interior grid, shown below. We collect the result for the purpose of this tutorial, since grid returns an iterable to avoid allocations.

julia> collect(grid(basis))5-element Vector{Float64}:
 3.9021130325903073
 3.1755705045849463
 2.0
 0.8244294954150537
 0.09788696740969272

We can show evaluate the basis functions at a given point. Again, it is an iterable, so we collect to show it here.

julia> collect(basis_at(basis, 0.41))5-element Vector{Float64}:
  1.0
 -0.795
  0.2640500000000001
  0.37516049999999984
 -0.8605551949999999

We can evaluate linear combination as directly, or via partial application.

julia> θ = [1, 0.5, 0.2, 0.3, 0.001];           # a vector of coefficients
julia> x = 0.410.41
julia> linear_combination(basis, θ, x)          # combination at some value0.7669975948050001
julia> linear_combination(basis, θ)(x)          # also as a callable0.7669975948050001

We can also evaluate derivatives of either the basis or the linear combination at a given point. Here we want the derivatives up to order 3.

julia> dx = (𝑑^2)(x)0.41 + 1.0⋅Δ + 0.0⋅Δ²
julia> collect(basis_at(basis, dx))5-element Vector{SpectralKit.𝑑Expansion{3, Float64}}:
 1.0 + 0.0⋅Δ + 0.0⋅Δ²
 -0.795 + 0.5⋅Δ + 0.0⋅Δ²
 0.2640500000000001 + -1.59⋅Δ + 1.0⋅Δ²
 0.37516049999999984 + 2.2921500000000004⋅Δ + -4.7700000000000005⋅Δ²
 -0.8605551949999999 + -1.6793580000000008⋅Δ + 11.168600000000001⋅Δ²
julia> fdx = linear_combination(basis, θ, dx)0.7669975948050001 + 0.617965642⋅Δ + -1.2198314000000001⋅Δ²
julia> fdx[0] # the value0.7669975948050001
julia> fdx[1] # the first derivative0.617965642

Having an approximation, we can embed it in a larger basis, extending the coefficients accordingly.

julia> basis2 = Chebyshev(EndpointGrid(), 8) ∘ transformation(basis)       # 8 Chebyshev polynomialsSpectralKit.TransformedBasis{Chebyshev{EndpointGrid}, BoundedLinear{Float64}}(Chebyshev polynomials (1st kind), EndpointGrid(), dimension: 8, (0.0,4.0) ↔ domain [linear transformation])
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

Multivariate (Smolyak) approximation example

We set up a Smolyak basis to approximate functions on $[-1,2] \times [-3, \infty]$, where the second dimension has a scaling of $3$.

using SpectralKit, StaticArrays
basis = smolyak_basis(Chebyshev, InteriorGrid2(), SmolyakParameters(3), 2)
ct = coordinate_transformations(BoundedLinear(-1, 2.0), SemiInfRational(-3.0, 3.0))
basis_t = basis ∘ ct

SpectralKit.TransformedBasis{SpectralKit.SmolyakBasis{SpectralKit.SmolyakIndices{2, 15, 3, 3, 4}, Chebyshev{InteriorGrid2}}, SpectralKit.CoordinateTransformations{Tuple{BoundedLinear{Float64}, SemiInfRational{Float64}}}}(Sparse multivariate basis on ℝ²
  Smolyak indexing, ∑bᵢ ≤ 3, all bᵢ ≤ 3, dimension 49
  using Chebyshev polynomials (1st kind), InteriorGrid2(), dimension: 15, coordinate transformations
  (-1.0,2.0) ↔ domain [linear transformation]
  (-3.0,∞) ↔ domain [rational transformation with scale 3.0])

Note how the basis can be combined with a transformation using ∘.

We will approximate the following function:

f2((x1, x2)) = exp(x1) + exp(-abs2(x2))

f2 (generic function with 1 method)

We find the coefficients by solving with the collocation matrix.

θ = collocation_matrix(basis_t) \ f2.(grid(basis_t))

49-element Vector{Float64}:
  2.7793220217324706
  3.2369886430572143
  1.1139902424577899
  0.2663479965030996
  0.04859825644525028
  0.007157295461750643
  0.0008829533666788979
  9.366852829718765e-5
  8.713768913578302e-6
  7.216598510789816e-7
  ⋮
 -2.8429038809505643e-16
  0.03010775018319918
  0.048411759340962084
 -0.03222472955157677
 -0.16530952363028062
  0.01639663318432438
  0.012259562724398244
 -0.01591087240587917
 -0.14085219047438763

Finally, we check the approximation at a point.

z = (0.5, 0.7)                            # evaluate at this point
isapprox(f2(z), linear_combination(basis_t, θ)(z), rtol = 0.005)

true

Infinite endpoints

Values and derivatives at $\pm\infty$ should provide the correct limits.

julia> using SpectralKit
julia> basis = Chebyshev(InteriorGrid(), 4) ∘ InfRational(0.0, 1.0)SpectralKit.TransformedBasis{Chebyshev{InteriorGrid}, InfRational{Float64}}(Chebyshev polynomials (1st kind), InteriorGrid(), dimension: 4, (-∞,∞) ↔ domain [rational transformation with center 0.0, scale 1.0])
julia> collect(basis_at(basis, 𝑑(Inf)))4-element Vector{SpectralKit.𝑑Expansion{2, Float64}}:
 1.0 + 0.0⋅Δ
 1.0 + 0.0⋅Δ
 1.0 + 0.0⋅Δ
 1.0 + 0.0⋅Δ
julia> collect(basis_at(basis, 𝑑(-Inf)))4-element Vector{SpectralKit.𝑑Expansion{2, Float64}}:
 1.0 + 0.0⋅Δ
 -1.0 + 0.0⋅Δ
 1.0 + -0.0⋅Δ
 -1.0 + 0.0⋅Δ

Constructing bases

Grid specifications

SpectralKit.EndpointGrid — Type

struct EndpointGrid <: SpectralKit.AbstractGrid

Grid that includes endpoints (eg Gauss-Lobatto).

Note

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