radau/lobatto

Author: DanielVandH

Project.toml

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

src/ClassicalOrthogonalPolynomials.jl

@@ -318,6 +318,76 @@ 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) 
+    endpt = T(endpt)
+    p0 = zero(T) 
+    p1 = one(T)
+    for i in 1:n 
+        # Evaluate the monic polynomials πₙ₋₁(endpt), πₙ(endpt)
+        _p1 = p0 
+        p0 = p1 
+        p1 = (endpt - α[i]) * p0 
+        i > 1 && (p1 -= β[i - 1]^2 * _p1)
+    end 
+    a = endpt - β[end]^2 * p0 / p1 
+    α′ = vcat(@view(α[begin:end-1]), a)
+    J′ = SymTridiagonal(α′, β)
+    x, w = golubwelsch(J′)
+    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
+
+function gausslobatto(P, n::Integer)
+    ## 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 = zero(T) 
+    p1a = one(T) 
+    p0b = zero(T) 
+    p1b = one(T) 
+    for i in 1:(n+1)
+        # Evaluate the monic polynomials πₙ₋₁(a), πₙ₋₁(b), πₙ(a), πₙ(b)
+        _p1a = p0a 
+        _p1b = p0b 
+        p0a = p1a 
+        p0b = p1b
+        p1a = (a - α[i]) * p0a 
+        p1b = (b - α[i]) * p0b 
+        if i > 1 
+            p1a -= β[i - 1]^2 * _p1a 
+            p1b -= β[i - 1]^2 * _p1b
+        end
+    end
+    # 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
+
 # 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

test/runtests.jl

@@ -43,6 +43,7 @@ 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 +84,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-x62)
+        @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