# 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:

1. Legendre: P_n(x)
2. Chebyshev (1st kind, 2nd kind): T_n(x), U_n(x)
3. Ultraspherical: C_n^{(λ)}(x)
4. Jacobi: P_n^{(a,b)}(x)
5. Laguerre: L_n^{(α)}(x)
6. Hermite: H_n(x)

## Evaluation

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

julia> using ClassicalOrthogonalPolynomialsjulia> n, x = 5, 0.1;julia> legendrep(n, x) # P_n(x)0.17882875000000004julia> chebyshevt(n, x) # T_n(x) == cos(n*acos(x))0.48016000000000003julia> chebyshevu(n, x) # U_n(x) == sin((n+1)*acos(x))/sin(acos(x))0.56832julia> λ = 0.3; ultrasphericalc(n, λ, x) # C_n^(λ)(x)0.08578714248julia> a,b = 0.1,0.2; jacobip(n, a, b, x) # P_n^(a,b)(x)0.17459116797117197julia> laguerrel(n, x) # L_n(x)0.5483540833333335julia> α = 0.1; laguerrel(n, α, x) # L_n^(α)(x)0.7329166666666659julia> 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.laguerrelFunction
 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.

### Weights

ClassicalOrthogonalPolynomials.ChebyshevWeightType

ChebyshevWeight{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).

### Recurrences

ClassicalOrthogonalPolynomials.jacobimatrixFunction
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.recurrencecoefficientsFunction
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.qr_jacobimatrixFunction

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

cholesky_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.SetindexInterlaceType
SetindexInterlace(z, args...)

is an analogue of Basis for vector that replaces the ith index of z`, takes the union of the first axis, and the second axis is a blocked interlace of args.