# 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 polynomials`

`Chebyshev polynomials (1st kind), EndpointGrid(), dimension: 5`

`julia> is_function_basis(basis) # ie we support the interface below`

`true`

`julia> dimension(basis) # number of basis functions`

`5`

`julia> domain(basis) # domain`

`(-1, 1)`

`julia> grid(basis) # Gauss-Lobatto grid`

`SpectralKit.ChebyshevGridIterator{Float64, Chebyshev{EndpointGrid}}(Chebyshev polynomials (1st kind), EndpointGrid(), dimension: 5)`

`julia> collect(basis_at(basis, 0.41)) # iterator for basis functions at 0.41`

`5-element Vector{Float64}: 1.0 0.41 -0.6638000000000001 -0.9543159999999999 -0.11873911999999986`

`julia> collect(basis_at(basis, derivatives(0.41))) # values and 1st derivatives`

`5-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 coefficients`

`5-element Vector{Float64}: 1.0 0.5 0.2 0.3 0.001`

`julia> linear_combination(basis, θ, 0.41) # combination at some value`

`0.7858264608800001`

`julia> linear_combination(basis, θ)(0.41) # also as a callable`

`0.7858264608800001`

`julia> basis2 = Chebyshev(EndpointGrid(), 8) # 8 Chebyshev polynomials`

`Chebyshev 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 it`

`8-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)) end`

`f2 (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 coefficients`

`49-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).

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

`SpectralKit.InteriorGrid`

— Type`struct InteriorGrid <: SpectralKit.AbstractGrid`

Grid with interior points (eg Gauss-Chebyshev).

`SpectralKit.InteriorGrid2`

— Type`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.

`SpectralKit.to_pm1`

— Function`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

`SpectralKit.from_pm1`

— Function`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

`SpectralKit.BoundedLinear`

— Type`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.

`SpectralKit.InfRational`

— Type```
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$

`SpectralKit.SemiInfRational`

— Type```
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$

`SpectralKit.coordinate_transformations`

— Function```
coordinate_transformations(transformations)
```

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.

`SpectralKit.Chebyshev`

— Type`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).

`SpectralKit.SmolyakParameters`

— Type```
SmolyakParameters(B)
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.

`SpectralKit.smolyak_basis`

— Function```
smolyak_basis(
univariate_family,
grid_kind,
smolyak_parameters,
_
)
```

Create a sparse Smolyak basis.

**Arguments**

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

**Example**

```
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)
81
julia> domain(basis)
((-1, 1), (-1, 1))
```

**Properties**

*Grids nest*: increasing arguments of `SmolyakParameters`

result in a refined grid that contains points of the cruder grid.

## Using bases

### Introspection

`SpectralKit.is_function_basis`

— Function`is_function_basis(F)`

`is_function_basis(f::F)`

Test if the argument is a *function basis*, supporting the following interface:

`domain`

for querying the domain,`dimension`

for the dimension,`basis_at`

for function evaluation,`grid`

to obtain collocation points.

`linear_combination`

and `collocation_matrix`

are also supported, building on the above.

Can be used on both types (preferred) and values (for convenience).

`SpectralKit.dimension`

— Function`dimension(basis)`

Return the dimension of `basis`

, a positive `Int`

.

`SpectralKit.domain`

— Function`domain(basis)`

The domain of a function basis. A tuple of numbers (of arbitrary type, but usually `Float64`

), or a tuple of domains by coordinate.

### Evaluation

`SpectralKit.basis_at`

— Function`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, `Tuple`

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

`SpectralKit.linear_combination`

— Function```
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

`SpectralKit.grid`

— Function`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.

`SpectralKit.collocation_matrix`

— Function```
collocation_matrix(basis)
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

`SpectralKit.augment_coefficients`

— Function`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`

.

`SpectralKit.is_subset_basis`

— Function```
is_subset_basis(basis1, basis2)
```

Return a `Bool`

indicating whether coefficients in `basis1`

can be augmented to `basis2`

with `augment_coefficients`

.

`true`

does not mean that coefficients from `basis1`

can just be padded with zeros, since they may be in different positions. Always use `augment_coefficients`

.

## Derivatives

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.

`SpectralKit.derivatives`

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

Obtain `N`

derivatives (and the function value) at a scalar `x`

. The `i`

th 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
0.0
```

**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
4.0
```

## Internals

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

`SpectralKit.UnivariateTransformation`

— Type`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

`SpectralKit.gridpoint`

— Function```
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.