diff --git a/Project.toml b/Project.toml index 162d780..bb2f0dd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ClassicalOrthogonalPolynomials" uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd" authors = ["Sheehan Olver "] -version = "0.13.4" +version = "0.13.5" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" diff --git a/src/ClassicalOrthogonalPolynomials.jl b/src/ClassicalOrthogonalPolynomials.jl index 4db1af5..57180bf 100644 --- a/src/ClassicalOrthogonalPolynomials.jl +++ b/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 diff --git a/src/monic.jl b/src/monic.jl new file mode 100644 index 0000000..70c10a7 --- /dev/null +++ b/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 + diff --git a/test/runtests.jl b/test/runtests.jl index cd75d99..8148153 100644 --- a/test/runtests.jl +++ b/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 \ No newline at end of file diff --git a/test/test_jacobi.jl b/test/test_jacobi.jl index e24676a..7ad1866 100644 --- a/test/test_jacobi.jl +++ b/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) diff --git a/test/test_lobattoradau.jl b/test/test_lobattoradau.jl new file mode 100644 index 0000000..30d66b3 --- /dev/null +++ b/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 diff --git a/test/test_monic.jl b/test/test_monic.jl new file mode 100644 index 0000000..deae45d --- /dev/null +++ b/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 \ No newline at end of file