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 ClassicalOrthogonalPolynomials
julia> n, x = 5, 0.1;
julia> legendrep(n, x) # P_n(x)0.17882875
julia> chebyshevt(n, x) # T_n(x) == cos(n*acos(x))0.48016
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.17459116797117194
julia> laguerrel(n, x) # L_n(x)0.5483540833333331
julia> α = 0.1; laguerrel(n, α, x) # L_n^(α)(x)0.732916666666666
julia> 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 — Type

Chebyshev{kind,T}()

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