Implement gausslobatto / gaussradau and Monic (#199)

* radau/lobatto

* Take the first iteration out of the loop to avoid if

* Implement monic

* typo

* Recurrence coefficients

* Doesn't really do much

* Typo

Author: DanielVandH

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.13.4"
+version = "0.13.5"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

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,6 +319,50 @@ end
 
 golubwelsch(V::SubQuasiArray) = golubwelsch(parent(V), maximum(parentindices(V)[2]))
 
+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, p1 = P[endpt,n:n+1]
+    a = endpt - β[end]^2 * p0 / p1
+    α′ = vcat(@view(α[begin:end-1]), a)
+    J′ = SymTridiagonal(α′, β)
+    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
+    end
+    return x, w
+end
+gaussradau(P, n::Integer, endpt) = gaussradau(Monic(P), n, endpt)
+
+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 = jacobimatrix(P.P, n + 2)
+    α, β = diagonaldata(J), supdiagonaldata(J)
+    p0a, p1a = P[a, (n+1):(n+2)]
+    p0b, p1b = P[b, (n+1):(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 / Δ
+    α′ = 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[end] = b # avoid rounding errors. Doesn't affect the actual result but avoids e.g. domain errors
+    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.
 function plan_transform_layout(::AbstractOPLayout, P, szs::NTuple{N,Int}, dims=ntuple(identity,Val(N))) where N

src/monic.jl

@@ -0,0 +1,38 @@
+struct Monic{T,OPs<:AbstractQuasiMatrix{T},NL} <: OrthogonalPolynomial{T}
+    P::Normalized{T,OPs,NL}
+    α::AbstractVector{T} # diagonal of Jacobi matrix of P
+    β::AbstractVector{T} # squared supdiagonal of Jacobi matrix of P
+end
+
+Monic(P::AbstractQuasiMatrix) = Monic(Normalized(P))
+function Monic(P::Normalized)
+    X = jacobimatrix(P)
+    α = diagonaldata(X)
+    β = supdiagonaldata(X)
+    return Monic(P, α, β.^2)
+end
+Monic(P::Monic) = Monic(P.P, P.α, 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 recurrencecoefficients(P::Monic{T}) where {T}
+    α = P.α
+    β = P.β
+    return _monicrecurrencecoefficients(α, β) # function barrier 
+end
+function _monicrecurrencecoefficients(α::AbstractVector{T}, β) where {T}
+    A = Ones{T}(∞)
+    B = -α
+    C = Vcat(zero(T), β)
+    return A, B, C
+end
+

test/runtests.jl

@@ -39,10 +39,12 @@ 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")
 include("test_roots.jl")
+include("test_lobattoradau.jl")
 
 @testset "Auto-diff" begin
     U = Ultraspherical(1)
@@ -83,4 +85,4 @@ end
 @testset "Issue #179" begin
     @test startswith(sprint(show, MIME"text/plain"(), Chebyshev()[0.3, :]; context=(:compact=>true, :limit=>true)), "ℵ₀-element view(::ChebyshevT{Float64}, 0.3, :)")
     @test startswith(sprint(show, MIME"text/plain"(), Jacobi(0.2, 0.5)[-0.7, :]; context=(:compact=>true, :limit=>true)), "ℵ₀-element view(::Jacobi{Float64}, -0.7, :)")
-end
+end

test/test_jacobi.jl

@@ -81,7 +81,7 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
 
             M = P[x,1:3]'Diagonal(w)*P[x,1:3]
             @test M ≈ Diagonal(M)
-            x,w = gaussradau(3,a,b)
+            x,w = FastGaussQuadrature.gaussradau(3,a,b)
             M = P[x,1:3]'Diagonal(w)*P[x,1:3]
             @test M ≈ Diagonal(M)
 

test/test_lobattoradau.jl

@@ -0,0 +1,91 @@
+using ClassicalOrthogonalPolynomials, FastGaussQuadrature
+const COP = ClassicalOrthogonalPolynomials
+const FGQ = FastGaussQuadrature
+using Test
+using ClassicalOrthogonalPolynomials: symtridiagonalize
+using LinearAlgebra
+
+@testset "gaussradau" begin
+    @testset "Compare with FastGaussQuadrature" begin
+        x1, w1 = COP.gaussradau(Legendre(), 5, -1.0)
+        x2, w2 = FGQ.gaussradau(6)
+        @test x1 ≈ x2 && w1 ≈ w2
+        @test x1[1] == -1
+
+        x1, w1 = COP.gaussradau(Jacobi(1.0, 3.5), 25, -1.0)
+        x2, w2 = FGQ.gaussradau(26, 1.0, 3.5)
+        @test x1 ≈ x2 && w1 ≈ w2
+        @test x1[1] == -1
+
+        I0, I1 = COP.ChebyshevInterval(), COP.UnitInterval()
+        P = Jacobi(2.0, 0.0)[COP.affine(I1, I0), :]
+        x1, w1 = COP.gaussradau(P, 18, 0.0)
+        x2, w2 = FGQ.gaussradau(19, 2.0, 0.0)
+        @test 2x1 .- 1 ≈ x2 && 2w1 ≈ w2
+
+        x1, w1 = COP.gaussradau(Jacobi(1 / 2, 1 / 2), 4, 1.0)
+        x2, w2 = FGQ.gaussradau(5, 1 / 2, 1 / 2)
+        @test sort(-x1) ≈ x2
+        @test_broken w1 ≈ w2 # What happens to the weights when inverting the interval?
+    end
+
+    @testset "Example 3.5 in Gautschi (2004)'s book" begin
+        P = Laguerre(3.0)
+        n = 5
+        J = symtridiagonalize(jacobimatrix(P))[1:(n-1), 1:(n-1)]
+        _J = zeros(n, n)
+        _J[1:n-1, 1:n-1] .= J
+        _J[n-1, n] = sqrt((n - 1) * (n - 1 + P.α))
+        _J[n, n-1] = _J[n-1, n]
+        _J[n, n] = n - 1
+        x, V = eigen(_J)
+        w = 6V[1, :] .^ 2
+        xx, ww = COP.gaussradau(P, n - 1, 0.0)
+        @test xx ≈ x && ww ≈ w
+    end
+
+    @testset "Some numerical integration" begin
+        f = x -> 2x + 7x^2 + 10x^3 + exp(-x)
+        x, w = COP.gaussradau(Chebyshev(), 10, -1.0)
+        @test dot(f.(x), w) ≈ 14.97303754807069897 # integral of (2x + 7x^2 + 10x^3 + exp(-x))/sqrt(1-x^2)
+        @test x[1] == -1
+        
+        f = x -> -1.0 + 5x^6
+        x, w = COP.gaussradau(Jacobi(-1/2, -1/2), 2, 1.0)
+        @test dot(f.(x), w) ≈ 9π/16
+        @test x[end] == 1 
+        @test length(x) == 3
+    end
+end
+
+@testset "gausslobatto" begin
+    @testset "Compare with FastGaussQuadrature" begin
+        x1, w1 = COP.gausslobatto(Legendre(), 5)
+        x2, w2 = FGQ.gausslobatto(7)
+        @test x1 ≈ x2 && w1 ≈ w2
+        @test x1[1] == -1
+        @test x1[end] == 1
+
+        I0, I1 = COP.ChebyshevInterval(), COP.UnitInterval()
+        P = Legendre()[COP.affine(I1, I0), :]
+        x1, w1 = COP.gausslobatto(P, 18)
+        x2, w2 = FGQ.gausslobatto(20)
+        @test 2x1 .- 1 ≈ x2 && 2w1 ≈ w2
+    end
+
+    @testset "Some numerical integration" begin
+        f = x -> 2x + 7x^2 + 10x^3 + exp(-x)
+        x, w = COP.gausslobatto(Chebyshev(), 10)
+        @test dot(f.(x), w) ≈ 14.97303754807069897
+        @test x[1] == -1
+        @test x[end] == 1 
+        @test length(x) == 12
+        
+        f = x -> -1.0 + 5x^6
+        x, w = COP.gausslobatto(Jacobi(-1/2, -1/2), 4)
+        @test dot(f.(x), w) ≈ 9π/16
+        @test x[1]==-1
+        @test x[end] == 1 
+        @test length(x) == 6
+    end
+end

test/test_monic.jl

@@ -0,0 +1,58 @@
+using ClassicalOrthogonalPolynomials
+using Test
+using ClassicalOrthogonalPolynomials: Monic, _p0, orthogonalityweight, recurrencecoefficients
+
+@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 "evaluation" 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.2, 5] # no longer inferred
+    end
+end