Implement monic

Author: DanielVandH

src/ClassicalOrthogonalPolynomials.jl

@@ -285,6 +285,7 @@ end
 include("clenshaw.jl")
 include("ratios.jl")
 include("normalized.jl")
+include("monic.jl")
 include("lanczos.jl")
 include("choleskyQR.jl")
 
@@ -318,70 +319,50 @@ end
 
 golubwelsch(V::SubQuasiArray) = golubwelsch(parent(V), maximum(parentindices(V)[2]))
 
-function gaussradau(P, n::Integer, endpt)
-    ## See Thm 3.2 in Gautschi (2004)'s book
-    # n is the number of interior points
-    J = symtridiagonalize(jacobimatrix(P, n + 1))
-    α = diagonaldata(J) 
-    β = supdiagonaldata(J) 
-    T = eltype(P) 
+function gaussradau(P::Monic{T}, n::Integer, endpt) where {T}
+    ## See Thm 3.2 in Gautschi (2004)'s book 
+    # n is number of interior points 
+    J = jacobimatrix(P.P, n + 1)
+    α, β = diagonaldata(J), supdiagonaldata(J)
     endpt = T(endpt)
-    p0 = one(T) 
-    p1 = (endpt - α[1]) * p0
-    for i in 2:n 
-        # Evaluate the monic polynomials πₙ₋₁(endpt), πₙ(endpt)
-        _p1 = p0 
-        p0 = p1 
-        p1 = (endpt - α[i]) * p0 - β[i - 1]^2 * _p1
-    end 
-    a = endpt - β[end]^2 * p0 / p1 
+    p0 = P[endpt, n] # πₙ₋₁
+    p1 = P[endpt, n+1] # πₙ. Could be faster by computing p1 from p0 and πₙ₋₂, but the cost is tiny relative to eigen()
+    a = (endpt - β[end]^2 * p0 / p1)::T
     α′ = vcat(@view(α[begin:end-1]), a)
     J′ = SymTridiagonal(α′, β)
-    x, w = golubwelsch(J′)
+    x, w = golubwelsch(J′) # not inferred
     w .*= sum(orthogonalityweight(P))
     if endpt ≤ axes(P, 1)[begin]
         x[1] = endpt # avoid rounding errors
-    else 
-        x[end] = endpt 
+    else
+        x[end] = endpt
     end
     return x, w
 end
+gaussradau(P, n::Integer, endpt) = gaussradau(Monic(P), n, endpt)
 
-function gausslobatto(P, n::Integer)
-    ## See Thm 3.6 of Gautschi (2004)'s book
+function gausslobatto(P::Monic{T}, n::Integer) where {T}
+    ## See Thm 3.6 of Gautschi (2004)'s book 
     # n is the number of interior points 
     a, b = axes(P, 1)[begin], axes(P, 1)[end]
-    J = symtridiagonalize(jacobimatrix(P, n + 2))
-    α = diagonaldata(J) 
-    β = supdiagonaldata(J) 
-    T = eltype(P) 
-    a, b = T(a), T(b)
-    p0a = one(T) 
-    p0b = one(T) 
-    p1a = (a - α[1]) * p0a 
-    p1b = (b - α[1]) * p0b
-    for i in 2:(n+1)
-        # Evaluate the monic polynomials πₙ₋₁(a), πₙ₋₁(b), πₙ(a), πₙ(b)
-        _p1a = p0a 
-        _p1b = p0b 
-        p0a = p1a 
-        p0b = p1b
-        p1a = (a - α[i]) * p0a - β[i - 1]^2 * _p1a 
-        p1b = (b - α[i]) * p0b - β[i - 1]^2 * _p1b
-    end
+    J = jacobimatrix(P.P, n + 2)
+    α, β = diagonaldata(J), supdiagonaldata(J)
+    p0a, p0b = P[a, n+1], P[b, n+1]
+    p1a, p1b = P[a, n+2], P[b, n+2]
     # Solve Eq. 3.1.2.8 
     Δ = p1a * p0b - p1b * p0a # This could underflow/overflow for large n 
-    aext = (a * p1a * p0b - b * p1b * p0a) / Δ 
-    bext = (b - a) * p1a * p1b / Δ 
+    aext = (a * p1a * p0b - b * p1b * p0a) / Δ
+    bext = (b - a) * p1a * p1b / Δ
     α′ = vcat(@view(α[begin:end-1]), aext)
     β′ = vcat(@view(β[begin:end-1]), sqrt(bext)) # LazyBandedMatrices doesn't like when we use Vcat(array, scalar) apparently
     J′ = LazyBandedMatrices.SymTridiagonal(α′, β′)
     x, w = golubwelsch(J′)
     w .*= sum(orthogonalityweight(P))
-    x[1] = a 
+    x[1] = a
     x[end] = b # avoid rounding errors. Doesn't affect the actual result but avoids e.g. domain errors
-    return x, w 
+    return x, w
 end
+gausslobatto(P, n::Integer) = gausslobatto(Monic(P), n)
 
 # Default is Golub–Welsch expansion
 # note this computes the grid an extra time.

src/monic.jl

@@ -0,0 +1,42 @@
+struct Monic{T,OPs<:AbstractQuasiMatrix{T},NL} <: OrthogonalPolynomial{T}
+    P::Normalized{T,OPs,NL}
+    # X::AbstractMatrix{T} # Should X be stored? Probably not. 
+    # X would be the Jacobi matrix of the normalised polynomials, not for the monic polynomials.
+end
+# Will need to figure out what this should be exactly. 
+# Consider this internal for now until it stabilises.
+
+Monic(P::AbstractQuasiMatrix) = Monic(Normalized(P))
+Monic(P::Monic) = P
+
+Normalized(P::Monic) = P.P
+
+axes(P::Monic) = axes(P.P)
+
+orthogonalityweight(P::Monic) = orthogonalityweight(P.P)
+
+_p0(::Monic{T}) where {T} = one(T)
+
+show(io::IO, P::Monic) = print(io, "Monic($(P.P.P))")
+show(io::IO, ::MIME"text/plain", P::Monic) = show(io, P)
+
+function getindex(P::Monic{T}, x::Number, n::Int)::T where {T}
+    # TODO: Rewrite this to be more efficient using forwardrecurrence!
+    p0 = _p0(P)
+    n == 1 && return p0
+    t = convert(T, x)
+    J = jacobimatrix(P.P, n)
+    α = diagonaldata(J)
+    β = supdiagonaldata(J)
+    p1 = ((t - α[1]) * p0)::T
+    n == 2 && return p1
+    for i in 2:(n-1)
+        _p1 = p0::T
+        p0 = p1::T
+        p1 = ((t - α[i]) * p1 - β[i-1]^2 * _p1)::T
+    end
+    return p1
+end
+# Should a method be written that makes this more efficient when requesting multiple n?
+
+

test/runtests.jl

@@ -39,6 +39,7 @@ include("test_fourier.jl")
 include("test_odes.jl")
 include("test_ratios.jl")
 include("test_normalized.jl")
+include("test_monic.jl")
 include("test_lanczos.jl")
 include("test_interlace.jl")
 include("test_choleskyQR.jl")

test/test_monic.jl

@@ -0,0 +1,58 @@
+using ClassicalOrthogonalPolynomials
+using Test
+using ClassicalOrthogonalPolynomials: Monic, _p0, orthogonalityweight
+
+@testset "Basic definition" begin
+    P1 = Legendre()
+    P2 = Normalized(P1)
+    P3 = Monic(P1)
+    @test P3.P == P2
+    @test Monic(P3) === P3
+    @test axes(P3) == axes(Legendre())
+    @test Normalized(P3) === P3.P
+    @test _p0(P3) == 1
+    @test orthogonalityweight(P3) == orthogonalityweight(P1)
+    @test sprint(show, MIME"text/plain"(), P3) == "Monic(Legendre())"
+end
+
+@testset "getindex" begin
+    function _pochhammer(x, n)
+        y = one(x)
+        for i in 0:(n-1)
+            y *= (x + i)
+        end
+        return y
+    end
+    jacobi_kn = (α, β, n) -> _pochhammer(n + α + β + 1, n) / (2.0^n * factorial(n))
+    ultra_kn = (λ, n) -> 2^n * _pochhammer(λ, n) / factorial(n)
+    chebt_kn = n -> n == 0 ? 1.0 : 2.0 .^ (n - 1)
+    chebu_kn = n -> 2.0^n
+    leg_kn = n -> 2.0^n * _pochhammer(1 / 2, n) / factorial(n)
+    lag_kn = n -> (-1)^n / factorial(n)
+    herm_kn = n -> 2.0^n
+    _Jacobi(α, β, x, n) = Jacobi(α, β)[x, n+1] / jacobi_kn(α, β, n)
+    _Ultraspherical(λ, x, n) = Ultraspherical(λ)[x, n+1] / ultra_kn(λ, n)
+    _ChebyshevT(x, n) = ChebyshevT()[x, n+1] / chebt_kn(n)
+    _ChebyshevU(x, n) = ChebyshevU()[x, n+1] / chebu_kn(n)
+    _Legendre(x, n) = Legendre()[x, n+1] / leg_kn(n)
+    _Laguerre(α, x, n) = Laguerre(α)[x, n+1] / lag_kn(n)
+    _Hermite(x, n) = Hermite()[x, n+1] / herm_kn(n)
+    Ps = [
+        Jacobi(2.0, 5.0) (x, n)->_Jacobi(2.0, 5.0, x, n)
+        Ultraspherical(1.7) (x, n)->_Ultraspherical(1.7, x, n)
+        ChebyshevT() _ChebyshevT
+        ChebyshevU() _ChebyshevU
+        Legendre() _Legendre
+        Laguerre(1.5) (x, n)->_Laguerre(1.5, x, n)
+        Hermite() _Hermite
+    ]
+    for (P, _P) in eachrow(Ps)
+        Q = Monic(P)
+        @test Q[0.2, 1] == 1.0
+        @test Q[0.25, 2] ≈ _P(0.25, 1)
+        @test Q[0.17, 3] ≈ _P(0.17, 2)
+        @test Q[0.4, 17] ≈ _P(0.4, 16)
+        @test Q[0.9, 21] ≈ _P(0.9, 20)
+        @inferred Q[0.4, 5]
+    end
+end