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)$, defined over $[-1, 1]$ with weight $w(x) = 1$.
- Chebyshev (1st kind, 2nd kind): $T_n(x)$ and $U_n(x)$, defined over $[-1, 1]$ with weights $w(x) = 1/\sqrt{1-x^2}$ and $w(x) = \sqrt{1-x^2}$, respectively.
- Ultraspherical: $C_n^{(\lambda)}(x)$, defined over $[-1, 1]$ with weight $w(x) = (1-x^2)^{\lambda-1/2}$.
- Jacobi: $P_n^{(a,b)}(x)$, defined over $[-1, 1]$ with weight $w(x) = (1-x)^a(1+x)^b$.
- Laguerre: $L_n^{(\alpha)}(x)$, defined over $[0, ∞)$ with weight $w(x) = x^\alpha \mathrm{e}^{-x}$.
- Hermite: $H_n(x)$, defined over $(-∞, ∞)$ with weight $w(x) = \mathrm{e}^{-x^2}$.
These special polynomials have many applications and can be used as a basis for any function given their domain conditions are met, however these polynomials have some advantages due to their formulation:
- Because of their relation to Laplace’s equation, Legendre polynomials can be useful as a basis for functions with spherical symmetry.
- Chebyshev polynomials are generally effective in reducing errors from numerical methods such as quadrature, interpolation, and approximation.
- Due to the flexibility of its parameters, Jacobi polynomials are capable of tailoring the behavior of an approximation around its endpoints, making these polynomials particularly useful in boundary value problems.
- Ultraspherical polynomials are advantageous in spectral methods for solving differential equations.
- Laguerre polynomials have a semi-infinite domain, therefore they are beneficial for problems involving exponential decay.
- Because of its weight function, Hermite polynomials can be useful in situations where functions display a Gaussian-like distribution.
These are just a few applications of these polynomials. They have many more uses across mathematics, physics, and engineering.
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.Legendre — TypeLegendre{T=Float64}(a,b)The quasi-matrix representing Legendre polynomials, where the first axes represents the interval and the second axes represents the polynomial index (starting from 1). See also legendre, legendrep, LegendreWeight and Jacobi.
Examples
julia> P = Legendre()
Legendre()
julia> typeof(P) # default eltype
Legendre{Float64}
julia> axes(P)
(Inclusion(-1.0 .. 1.0 (Chebyshev)), OneToInf())
julia> P[0,:] # Values of polynomials at x=0
ℵ₀-element view(::Legendre{Float64}, 0.0, :) with eltype Float64 with indices OneToInf():
1.0
0.0
-0.5
-0.0
0.375
0.0
-0.3125
-0.0
0.2734375
0.0
⋮
julia> P₀=P[:,1]; # P₀ is the first Legendre polynomial which is constant.
julia> P₀[0],P₀[0.5]
(1.0, 1.0)ClassicalOrthogonalPolynomials.legendrep — Function legendrep(n, z)computes the n-th Legendre polynomial at z.
ClassicalOrthogonalPolynomials.Jacobi — TypeJacobi{T}(a,b)
Jacobi(a,b)The quasi-matrix representing Jacobi polynomials, where the first axes represents the interval and the second axes represents the polynomial index (starting from 1). See also jacobi, jacobip and JacobiWeight.
The eltype, when not specified, will be converted to a floating point data type.
Examples
julia> J=Jacobi(0, 0) # The eltype will be converted to float
Jacobi(0, 0)
julia> axes(J)
(Inclusion(-1.0 .. 1.0 (Chebyshev)), OneToInf())
julia> J[0,:] # Values of polynomials at x=0
ℵ₀-element view(::Jacobi{Float64, Int64}, 0.0, :) with eltype Float64 with indices OneToInf():
1.0
0.0
-0.5
-0.0
0.375
0.0
-0.3125
-0.0
0.2734375
0.0
⋮
julia> J0=J[:,1]; # J0 is the first Jacobi polynomial which is constant.
julia> J0[0],J0[0.5]
(1.0, 1.0)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.LegendreWeight — TypeLegendreWeight{T}()
LegendreWeight()The quasi-vector representing the Legendre weight function (const 1) on [-1,1]. See also legendreweight and Legendre.
Examples
julia> w = LegendreWeight()
LegendreWeight{Float64}()
julia> w[0.5]
1.0
julia> axes(w)
(Inclusion(-1.0 .. 1.0 (Chebyshev)),)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.JacobiWeight — TypeJacobiWeight{T}(a,b)
JacobiWeight(a,b)The quasi-vector representing the Jacobi weight function $(1-x)^a (1+x)^b$ on $[-1,1]$. See also jacobiweight and Jacobi.
Examples
julia> J=JacobiWeight(1.0,1.0)
(1-x)^1.0 * (1+x)^1.0 on -1..1
julia> J[0.5]
0.75
julia> axes(J)
(Inclusion(-1.0 .. 1.0 (Chebyshev)),)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)
Affine-mapped
ClassicalOrthogonalPolynomials.legendre — Functionlegendre(d::AbstractInterval)The Legendre polynomials affine-mapped to interval d.
Examples
julia> P = legendre(0..1)
Legendre() affine mapped to 0 .. 1
julia> axes(P)
(Inclusion(0 .. 1), OneToInf())
julia> P[0.5,:]
ℵ₀-element view(::Legendre{Float64}, 0.0, :) with eltype Float64 with indices OneToInf():
1.0
0.0
-0.5
-0.0
0.375
0.0
-0.3125
-0.0
0.2734375
0.0
⋮ClassicalOrthogonalPolynomials.jacobi — Functionjacobi(a,b, d::AbstractInterval)The Jacobi polynomials affine-mapped to interval d.
Examples
julia> J = jacobi(1, 1, 0..1)
Jacobi(1, 1) affine mapped to 0 .. 1
julia> axes(J)
(Inclusion(0 .. 1), OneToInf())
julia> J[0,:]
ℵ₀-element view(::Jacobi{Float64, Int64}, -1.0, :) with eltype Float64 with indices OneToInf():
1.0
-2.0
3.0
-4.0
5.0
-6.0
7.0
-8.0
9.0
-10.0
⋮ClassicalOrthogonalPolynomials.legendreweight — Functionlegendreweight(d::AbstractInterval)The LegendreWeight affine-mapped to interval d.
Examples
julia> w = legendreweight(0..1)
LegendreWeight{Float64}() affine mapped to 0 .. 1
julia> axes(w)
(Inclusion(0 .. 1),)
julia> w[0]
1.0ClassicalOrthogonalPolynomials.jacobiweight — Functionjacobiweight(a,b, d::AbstractInterval)The JacobiWeight affine-mapped to interval d.
Examples
julia> J = jacobiweight(1, 1, 0..1)
(1-x)^1 * (1+x)^1 on -1..1 affine mapped to 0 .. 1
julia> axes(J)
(Inclusion(0 .. 1),)
julia> J[0.5]
1.0Recurrences
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.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.ShuffledFFT — TypeGives a shuffled version of the 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.ShuffledR2HC — TypeGives a shuffled version of the real FFT, with order 1,sin(θ),cos(θ),sin(2θ)…
ClassicalOrthogonalPolynomials.ShuffledIR2HC — TypeGives a shuffled version of the real 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.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.MappedOPLayout — TypeMappedOPLayoutrepresents an OP that is (usually affine) mapped OP
ClassicalOrthogonalPolynomials.WeightedOPLayout — TypeWeightedOPLayoutrepresents an OP multiplied by its orthogonality weight.
ClassicalOrthogonalPolynomials.legendre_grammatrix — Functionlegendre_grammatrixcomputes the grammatrix by first re-expanding in Legendre
ClassicalOrthogonalPolynomials.weightedgrammatrix — Functiongives the inner products of OPs with their weight, i.e., Weighted(P)'P.
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 — TypeRepresents orthonormal polynomials defined via 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.