# ClassicalOrthogonalPolynomials.jl

A Julia package for classical orthogonal polynomials and expansions

## Definitions

We follow the Digital Library of Mathematical Functions, which defines the following classical orthogonal polynomials:

- Legendre:
`P_n(x)`

- Chebyshev (1st kind, 2nd kind):
`T_n(x)`

,`U_n(x)`

- Ultraspherical:
`C_n^{(λ)}(x)`

- Jacobi:
`P_n^{(a,b)}(x)`

- Laguerre:
`L_n^{(α)}(x)`

- Hermite:
`H_n(x)`

## Evaluation

The simplest usage of this package is to evaluate classical orthogonal polynomials:

`julia> using ClassicalOrthogonalPolynomials`

`julia> n, x = 5, 0.1;`

`julia> legendrep(n, x) # P_n(x)`

`0.17882875000000004`

`julia> chebyshevt(n, x) # T_n(x) == cos(n*acos(x))`

`0.48016000000000003`

`julia> chebyshevu(n, x) # U_n(x) == sin((n+1)*acos(x))/sin(acos(x))`

`0.56832`

`julia> λ = 0.3; ultrasphericalc(n, λ, x) # C_n^(λ)(x)`

`0.08578714248`

`julia> a,b = 0.1,0.2; jacobip(n, a, b, x) # P_n^(a,b)(x)`

`0.17459116797117197`

`julia> laguerrel(n, x) # L_n(x)`

`0.5483540833333335`

`julia> α = 0.1; laguerrel(n, α, x) # L_n^(α)(x)`

`0.7329166666666659`

`julia> hermiteh(n, x) # H_n(x)`

`11.840320000000002`

## Continuum arrays

For expansions, recurrence relationships, and other operations linked with linear equations, it is useful to treat the families of orthogonal polynomials as *continuum arrays*, as implemented in ContinuumArrays.jl. The continuum arrays implementation is accessed as follows:

`julia> T = ChebyshevT() # Or just Chebyshev()`

`ChebyshevT()`

`julia> axes(T) # [-1,1] by 1:∞`

`(Inclusion(-1.0 .. 1.0 (Chebyshev)), OneToInf())`

`julia> T[x, n+1] # T_n(x) = cos(n*acos(x))`

`0.48016000000000003`

We can thereby access many points and indices efficiently using array-like language:

`julia> x = range(-1, 1; length=1000);`

`julia> T[x, 1:1000] # [T_j(x[k]) for k=1:1000, j=0:999]`

`1000×1000 Matrix{Float64}: 1.0 -1.0 1.0 -1.0 … -1.0 1.0 -1.0 1.0 -0.997998 0.992 -0.98203 -0.964818 0.946258 -0.92391 1.0 -0.995996 0.984016 -0.964156 -0.282676 0.195792 -0.107341 1.0 -0.993994 0.976048 -0.946378 0.807761 -0.867423 0.916664 1.0 -0.991992 0.968096 -0.928695 -0.827934 0.750471 -0.660989 1.0 -0.98999 0.96016 -0.911108 … 0.982736 -0.999011 0.995286 1.0 -0.987988 0.952241 -0.893616 0.73242 -0.618409 0.489542 1.0 -0.985986 0.944337 -0.87622 0.822718 -0.716355 0.589914 1.0 -0.983984 0.936449 -0.858918 0.923481 -0.977078 0.999377 1.0 -0.981982 0.928577 -0.84171 -0.495939 0.322905 -0.138236 ⋮ ⋱ 1.0 0.983984 0.936449 0.858918 -0.923481 -0.977078 -0.999377 1.0 0.985986 0.944337 0.87622 -0.822718 -0.716355 -0.589914 1.0 0.987988 0.952241 0.893616 -0.73242 -0.618409 -0.489542 1.0 0.98999 0.96016 0.911108 -0.982736 -0.999011 -0.995286 1.0 0.991992 0.968096 0.928695 … 0.827934 0.750471 0.660989 1.0 0.993994 0.976048 0.946378 -0.807761 -0.867423 -0.916664 1.0 0.995996 0.984016 0.964156 0.282676 0.195792 0.107341 1.0 0.997998 0.992 0.98203 0.964818 0.946258 0.92391 1.0 1.0 1.0 1.0 1.0 1.0 1.0`

## Expansions

We view a function expansion in say Chebyshev polynomials in terms of continuum arrays as follows:

\[f(x) = \sum_{k=0}^∞ c_k T_k(x) = \begin{bmatrix}T_0(x) | T_1(x) | … \end{bmatrix} \begin{bmatrix}c_0\\ c_1 \\ \vdots \end{bmatrix} = T[x,:] * 𝐜\]

To be more precise, we think of functions as continuum-vectors. Here is a simple example:

`julia> f = T * [1; 2; 3; zeros(∞)]; # T_0(x) + T_1(x) + T_2(x)`

`julia> f[0.1]`

`-1.74`

To find the coefficients for a given function we consider this as the problem of finding `𝐜`

such that `T*𝐜 == f`

, that is:

`julia> T \ f`

`ℵ₀-element LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf(): 1.0 2.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋮`

For a function given only pointwise we broadcast over `x`

, e.g., the following are the coefficients of `\exp(x)`

:

`julia> x = axes(T, 1);`

`julia> c = T \ exp.(x)`

`vcat(14-element Vector{Float64}, ℵ₀-element FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}} with indices OneToInf()) with indices OneToInf(): 1.2660658777520084 1.1303182079849703 0.27149533953407656 0.04433684984866379 0.0054742404420936785 0.0005429263119139232 4.497732295427654e-5 3.19843646253781e-6 1.992124804817033e-7 1.1036771869970875e-8 ⋮`

`julia> f = T*c; f[0.1] # ≈ exp(0.1)`

`1.1051709180756477`

With a little cheeky usage of Julia's order-of-operations this can be written succicently as:

`julia> f = T / T \ exp.(x);`

`julia> f[0.1]`

`1.1051709180756477`

(Or for more clarity just write `T * (T \ exp.(x))`

.)

## Jacobi matrices

Orthogonal polynomials satisfy well-known three-term recurrences:

\[x p_n(x) = c_{n-1} p_{n-1}(x) + a_n p_n(x) + b_n p_{n+1}(x).\]

In continuum-array language this has the form of a comuting relationship:

\[x \begin{bmatrix} p_0 | p_1 | \cdots \end{bmatrix} = \begin{bmatrix} p_0 | p_1 | \cdots \end{bmatrix} \begin{bmatrix} a_0 & c_0 \\ b_0 & a_1 & c_1 \\ & b_1 & a_2 & \ddots \\ &&\ddots & \ddots \end{bmatrix}\]

We can therefore find the Jacobi matrix naturally as follows:

`julia> T \ (x .* T)`

`ℵ₀×ℵ₀ LazyBandedMatrices.Tridiagonal{Float64, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, FillArrays.Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf()×OneToInf(): 0.0 0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … 1.0 0.0 0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 0.5 ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋮ ⋮ ⋱`

Alternatively, just call `jacobimatrix(T)`

(noting its the transpose of the more traditional convention).

## Derivatives

The derivatives of classical orthogonal polynomials are also classical OPs, and this can be seen as follows:

`julia> U = ChebyshevU();`

`julia> D = Derivative(x);`

`julia> U\D*T`

`ℵ₀×ℵ₀ BandedMatrix{Float64} with bandwidths (-1, 1) with data 1×ℵ₀ adjoint(::InfiniteArrays.InfStepRange{Float64, Float64}) with eltype Float64 with indices Base.OneTo(1)×OneToInf() with indices OneToInf()×OneToInf(): ⋅ 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ 2.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 3.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 4.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 5.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 6.0 ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 7.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋮ ⋮ ⋱`

Similarly, the derivative of *weighted* classical OPs are weighted classical OPs:

`julia> Weighted(T)\D*Weighted(U)`

`ℵ₀×ℵ₀ BandedMatrix{Float64} with bandwidths (1, -1) with data 1×ℵ₀ adjoint(::InfiniteArrays.InfStepRange{Float64, Float64}) with eltype Float64 with indices Base.OneTo(1)×OneToInf() with indices OneToInf()×OneToInf(): ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -2.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -3.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -4.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -5.0 ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ -6.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋮ ⋮ ⋱`

## Other recurrence relationships

Many other sparse recurrence relationships are implemented. Here's one:

`julia> U\T`

`ℵ₀×ℵ₀ BandedMatrix{Float64} with bandwidths (0, 2) with data vcat(1×ℵ₀ FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf(), 1×ℵ₀ FillArrays.Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf(), hcat(1×1 Ones{Float64}, 1×ℵ₀ FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf()) with indices Base.OneTo(1)×OneToInf()) with indices Base.OneTo(3)×OneToInf() with indices OneToInf()×OneToInf(): 1.0 0.0 -0.5 ⋅ ⋅ ⋅ ⋅ … ⋅ 0.5 0.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 -0.5 ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 0.0 … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋮ ⋮ ⋱`

(Probably best to ignore the type signature 😅)

## Index

### Polynomials

`ClassicalOrthogonalPolynomials.Chebyshev`

— TypeChebyshev{kind,T}()

is a quasi-matrix representing Chebyshev polynomials of the specified kind (1, 2, 3, or 4) on -1..1.

`ClassicalOrthogonalPolynomials.chebyshevt`

— Function` chebyshevt(n, z)`

computes the `n`

-th Chebyshev polynomial of the first kind at `z`

.

`ClassicalOrthogonalPolynomials.chebyshevu`

— Function` chebyshevt(n, z)`

computes the `n`

-th Chebyshev polynomial of the second kind at `z`

.

`ClassicalOrthogonalPolynomials.legendrep`

— Function` legendrep(n, z)`

computes the `n`

-th Legendre polynomial at `z`

.

`ClassicalOrthogonalPolynomials.jacobip`

— Function` jacobip(n, a, b, z)`

computes the `n`

-th Jacobi polynomial, orthogonal with respec to `(1-x)^a*(1+x)^b`

, at `z`

.

`ClassicalOrthogonalPolynomials.laguerrel`

— Function` laguerrel(n, α, z)`

computes the `n`

-th generalized Laguerre polynomial, orthogonal with respec to `x^α * exp(-x)`

, at `z`

.

` laguerrel(n, z)`

computes the `n`

-th Laguerre polynomial, orthogonal with respec to `exp(-x)`

, at `z`

.

`ClassicalOrthogonalPolynomials.hermiteh`

— Function` hermiteh(n, z)`

computes the `n`

-th Hermite polynomial, orthogonal with respec to `exp(-x^2)`

, at `z`

.

### Weights

`ClassicalOrthogonalPolynomials.OrthonormalWeighted`

— Type`OrthonormalWeighted(P)`

is the orthonormal with respect to L^2 basis given by `sqrt.(orthogonalityweight(P)) .* Normalized(P)`

.

`ClassicalOrthogonalPolynomials.HermiteWeight`

— TypeHermiteWeight()

is a quasi-vector representing `exp(-x^2)`

on ℝ.

`ClassicalOrthogonalPolynomials.Weighted`

— Type`Weighted(P)`

is equivalent to `orthogonalityweight(P) .* P`

`ClassicalOrthogonalPolynomials.ChebyshevWeight`

— TypeChebyshevWeight{kind,T}()

is a quasi-vector representing the Chebyshev weight of the specified kind on -1..1. That is, `ChebyshevWeight{1}()`

represents `1/sqrt(1-x^2)`

, and `ChebyshevWeight{2}()`

represents `sqrt(1-x^2)`

.

`ClassicalOrthogonalPolynomials.LaguerreWeight`

— TypeLaguerreWeight(α)

is a quasi-vector representing `x^α * exp(-x)`

on `0..Inf`

.

`ClassicalOrthogonalPolynomials.HalfWeighted`

— Type`HalfWeighted{lr}(Jacobi(a,b))`

is equivalent to `JacobiWeight(a,0) .* Jacobi(a,b)`

(`lr = :a`

) or `JacobiWeight(0,b) .* Jacobi(a,b)`

(`lr = :b`

)

### Recurrences

`ClassicalOrthogonalPolynomials.normalizationconstant`

— Functionnormalizationconstant

gives the normalization constants so that the jacobi matrix is symmetric, that is, so we have orthonormal OPs:

`Q == P*normalizationconstant(P)`

`ClassicalOrthogonalPolynomials.OrthogonalPolynomialRatio`

— Type`OrthogonalPolynomialRatio(P,x)`

is a is equivalent to the vector `P[x,:] ./ P[x,2:end]`

`but built from the recurrence coefficients of`

P`.

`ClassicalOrthogonalPolynomials.Clenshaw`

— Type`Clenshaw(a, X)`

represents the operator `a(X)`

where a is a polynomial. Here `a`

is to stored as a quasi-vector.

`ClassicalOrthogonalPolynomials.singularities`

— Function`singularities(f)`

gives the singularity structure of an expansion, e.g., `JacobiWeight`

.

`ClassicalOrthogonalPolynomials.jacobimatrix`

— Function`jacobimatrix(P)`

returns the Jacobi matrix `X`

associated to a quasi-matrix of orthogonal polynomials satisfying

```
x = axes(P,1)
x*P == P*X
```

Note that `X`

is the transpose of the usual definition of the Jacobi matrix.

`ClassicalOrthogonalPolynomials.recurrencecoefficients`

— Function`recurrencecoefficients(P)`

returns a `(A,B,C)`

associated with the Orthogonal Polynomials P, satisfying for `x in axes(P,1)`

```
P[x,2] == (A[1]*x + B[1])*P[x,1]
P[x,n+1] == (A[n]*x + B[n])*P[x,n] - C[n]*P[x,n-1]
```

Note that `C[1]`

` is unused.

The relationship with the Jacobi matrix is:

```
1/A[n] == X[n+1,n]
-B[n]/A[n] == X[n,n]
C[n+1]/A[n+1] == X[n,n+1]
```

### Internal

`ClassicalOrthogonalPolynomials.ShuffledIR2HC`

— TypeGives a shuffled version of the real IFFT, with order 1,sin(θ),cos(θ),sin(2θ)…

`ClassicalOrthogonalPolynomials.ShuffledR2HC`

— TypeGives a shuffled version of the real FFT, with order 1,sin(θ),cos(θ),sin(2θ)…

`ClassicalOrthogonalPolynomials.ShuffledIFFT`

— TypeGives a shuffled version of the IFFT, with order 1,sin(θ),cos(θ),sin(2θ)…

`ClassicalOrthogonalPolynomials.qr_jacobimatrix`

— Functionqr_jacobimatrix(w, P)

returns the Jacobi matrix `X`

associated to a quasi-matrix of polynomials orthogonal with respect to `w(x)`

by computing a QR decomposition of the square root weight modification.

The resulting polynomials are orthonormal on the same domain as `P`

. The supplied `P`

must be normalized. Accepted inputs for `w`

are the target weight as a function or `sqrtW`

, representing the multiplication operator of square root weight modification on the basis `P`

.

The underlying QR approach allows two methods, one which uses the Q matrix and one which uses the R matrix. To change between methods, an optional argument :Q or :R may be supplied. The default is to use the Q method.

`ClassicalOrthogonalPolynomials.MappedOPLayout`

— Type`MappedOPLayout`

represents an OP that is (usually affine) mapped OP

`ClassicalOrthogonalPolynomials.cholesky_jacobimatrix`

— Functioncholesky_jacobimatrix(w, P)

returns the Jacobi matrix `X`

associated to a quasi-matrix of polynomials orthogonal with respect to `w(x)`

by computing a Cholesky decomposition of the weight modification.

The resulting polynomials are orthonormal on the same domain as `P`

. The supplied `P`

must be normalized. Accepted inputs are `w`

as a function or `W`

as an infinite matrix representing the weight modifier multiplication by the function `w / w_P`

on `P`

where `w_P`

is the orthogonality weight of `P`

.

`ClassicalOrthogonalPolynomials.AbstractNormalizedOPLayout`

— Type`AbstractNormalizedOPLayout`

represents OPs that are of the form P * R where P is another family of OPs and R is upper-triangular.

`ClassicalOrthogonalPolynomials.ShuffledFFT`

— TypeGives a shuffled version of the FFT, with order 1,sin(θ),cos(θ),sin(2θ)…

`ClassicalOrthogonalPolynomials.legendre_grammatrix`

— Function`legendre_grammatrix`

computes the grammatrix by first re-expanding in Legendre

`ClassicalOrthogonalPolynomials.WeightedOPLayout`

— Type`WeightedOPLayout`

represents an OP multiplied by its orthogonality weight.

`ClassicalOrthogonalPolynomials.interlace!`

— FunctionThis function implements the algorithm described in:

`P. Jain, "A simple in-place algorithm for in-shuffle," arXiv:0805.1598, 2008.`

`ClassicalOrthogonalPolynomials._tritrunc`

— Function`_tritrunc(X,n)`

does a square truncation of a tridiagonal matrix.

`ClassicalOrthogonalPolynomials.SetindexInterlace`

— Type`SetindexInterlace(z, args...)`

is an analogue of `Basis`

for vector that replaces the `i`

th index of `z`

, takes the union of the first axis, and the second axis is a blocked interlace of args.

`ClassicalOrthogonalPolynomials.ConvertedOrthogonalPolynomial`

— TypeRepresent an Orthogonal polynomial which has a conversion operator from P, that is, Q = P * inv(U).

`ClassicalOrthogonalPolynomials.PiecewiseInterlace`

— Type`PiecewiseInterlace(args...)`

is an analogue of `Basis`

that takes the union of the first axis, and the second axis is a blocked interlace of args. If there is overlap, it uses the first in order.