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

Introduction

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.

A basis is constructed using

a family on a fixed domain, eg Chebyshev,
a grid specification like InteriorGrid,
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.

Examples

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

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.