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 ClassicalOrthogonalPolynomialsjulia> n, x = 5, 0.1;julia> legendrep(n, x) # P_n(x)0.17882875julia> chebyshevt(n, x) # T_n(x) == cos(n*acos(x))0.48016julia> 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.17459116797117194julia> laguerrel(n, x) # L_n(x)0.5483540833333331julia> α = 0.1; laguerrel(n, α, x) # L_n^(α)(x)0.732916666666666julia> hermiteh(n, x) # H_n(x)11.84032
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.48016
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 — TypeOrthonormalWeighted(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 — TypeWeighted(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 — TypeHalfWeighted{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 — TypeOrthogonalPolynomialRatio(P,x)is a is equivalent to the vector P[x,:] ./ P[x,2:end]but built from the recurrence coefficients ofP`.
ClassicalOrthogonalPolynomials.Clenshaw — TypeClenshaw(a, X)represents the operator a(X) where a is a polynomial. Here a is to stored as a quasi-vector.
ClassicalOrthogonalPolynomials.singularities — Functionsingularities(f)gives the singularity structure of an expansion, e.g., JacobiWeight.
ClassicalOrthogonalPolynomials.jacobimatrix — Functionjacobimatrix(P)returns the Jacobi matrix X associated to a quasi-matrix of orthogonal polynomials satisfying
x = axes(P,1)
x*P == P*XNote that X is the transpose of the usual definition of the Jacobi matrix.
ClassicalOrthogonalPolynomials.recurrencecoefficients — Functionrecurrencecoefficients(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 — TypeMappedOPLayoutrepresents 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 — TypeAbstractNormalizedOPLayoutrepresents 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 — Functionlegendre_grammatrixcomputes the grammatrix by first re-expanding in Legendre
ClassicalOrthogonalPolynomials.WeightedOPLayout — TypeWeightedOPLayoutrepresents 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 — TypeSetindexInterlace(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.
ClassicalOrthogonalPolynomials.ConvertedOrthogonalPolynomial — TypeRepresent an Orthogonal polynomial which has a conversion operator from P, that is, Q = P * inv(U).
ClassicalOrthogonalPolynomials.PiecewiseInterlace — TypePiecewiseInterlace(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.