From e12006cf6523d64069c2bd843000ac4cc750f7ad Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Mon, 22 Jul 2024 22:06:07 +0100 Subject: [PATCH 1/7] radau/lobatto --- Project.toml | 2 +- src/ClassicalOrthogonalPolynomials.jl | 70 +++++++++++++++++++++ test/runtests.jl | 3 +- test/test_jacobi.jl | 2 +- test/test_lobattoradau.jl | 91 +++++++++++++++++++++++++++ 5 files changed, 165 insertions(+), 3 deletions(-) create mode 100644 test/test_lobattoradau.jl 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..ac4ffae 100644 --- a/src/ClassicalOrthogonalPolynomials.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index cd75d99..bc6a554 100644 --- a/test/runtests.jl +++ b/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 \ 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..8859169 --- /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-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 From d6f93759d7802a6d7a81f62bebbe4021778d56ef Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Mon, 22 Jul 2024 22:14:49 +0100 Subject: [PATCH 2/7] Take the first iteration out of the loop to avoid if --- src/ClassicalOrthogonalPolynomials.jl | 27 +++++++++++---------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/src/ClassicalOrthogonalPolynomials.jl b/src/ClassicalOrthogonalPolynomials.jl index ac4ffae..6f507dc 100644 --- a/src/ClassicalOrthogonalPolynomials.jl +++ b/src/ClassicalOrthogonalPolynomials.jl @@ -326,14 +326,13 @@ function gaussradau(P, n::Integer, endpt) β = supdiagonaldata(J) T = eltype(P) endpt = T(endpt) - p0 = zero(T) - p1 = one(T) - for i in 1:n + 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 && (p1 -= β[i - 1]^2 * _p1) + p1 = (endpt - α[i]) * p0 - β[i - 1]^2 * _p1 end a = endpt - β[end]^2 * p0 / p1 α′ = vcat(@view(α[begin:end-1]), a) @@ -357,22 +356,18 @@ function gausslobatto(P, n::Integer) β = 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) + 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 - p1b = (b - α[i]) * p0b - if i > 1 - p1a -= β[i - 1]^2 * _p1a - p1b -= β[i - 1]^2 * _p1b - end + p1a = (a - α[i]) * p0a - β[i - 1]^2 * _p1a + p1b = (b - α[i]) * p0b - β[i - 1]^2 * _p1b end # Solve Eq. 3.1.2.8 Δ = p1a * p0b - p1b * p0a # This could underflow/overflow for large n From 27fc2e2ae1861ce1660c5f476ddbac96e4bd6121 Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Tue, 23 Jul 2024 08:33:54 +0100 Subject: [PATCH 3/7] Implement monic --- src/ClassicalOrthogonalPolynomials.jl | 67 ++++++++++----------------- src/monic.jl | 42 +++++++++++++++++ test/runtests.jl | 1 + test/test_monic.jl | 58 +++++++++++++++++++++++ 4 files changed, 125 insertions(+), 43 deletions(-) create mode 100644 src/monic.jl create mode 100644 test/test_monic.jl diff --git a/src/ClassicalOrthogonalPolynomials.jl b/src/ClassicalOrthogonalPolynomials.jl index 6f507dc..7edc4ef 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,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. diff --git a/src/monic.jl b/src/monic.jl new file mode 100644 index 0000000..475051a --- /dev/null +++ b/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? + + diff --git a/test/runtests.jl b/test/runtests.jl index bc6a554..8148153 100644 --- a/test/runtests.jl +++ b/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") diff --git a/test/test_monic.jl b/test/test_monic.jl new file mode 100644 index 0000000..8d782b1 --- /dev/null +++ b/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 \ No newline at end of file From b739598f7b88b9f23c78fb7244693bdb32068b93 Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Tue, 23 Jul 2024 08:35:48 +0100 Subject: [PATCH 4/7] typo --- test/test_lobattoradau.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_lobattoradau.jl b/test/test_lobattoradau.jl index 8859169..30d66b3 100644 --- a/test/test_lobattoradau.jl +++ b/test/test_lobattoradau.jl @@ -47,7 +47,7 @@ using LinearAlgebra @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 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 From 1698661d694ddf12a67d0ead36faccdae8d2cd6d Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Tue, 23 Jul 2024 11:43:56 +0100 Subject: [PATCH 5/7] Recurrence coefficients --- src/ClassicalOrthogonalPolynomials.jl | 7 ++--- src/monic.jl | 42 ++++++++++++--------------- test/test_monic.jl | 6 ++-- 3 files changed, 25 insertions(+), 30 deletions(-) diff --git a/src/ClassicalOrthogonalPolynomials.jl b/src/ClassicalOrthogonalPolynomials.jl index 7edc4ef..279c247 100644 --- a/src/ClassicalOrthogonalPolynomials.jl +++ b/src/ClassicalOrthogonalPolynomials.jl @@ -325,8 +325,7 @@ function gaussradau(P::Monic{T}, n::Integer, endpt) where {T} J = jacobimatrix(P.P, n + 1) α, β = diagonaldata(J), supdiagonaldata(J) endpt = T(endpt) - 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() + p0, p1 = P[endpt,n:n+1] a = (endpt - β[end]^2 * p0 / p1)::T α′ = vcat(@view(α[begin:end-1]), a) J′ = SymTridiagonal(α′, β) @@ -347,8 +346,8 @@ function gausslobatto(P::Monic{T}, n::Integer) where {T} a, b = axes(P, 1)[begin], axes(P, 1)[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] + p0a, p0b = P[a, (n+1):(n+2)] + p1a, 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) / Δ diff --git a/src/monic.jl b/src/monic.jl index 475051a..70c10a7 100644 --- a/src/monic.jl +++ b/src/monic.jl @@ -1,13 +1,17 @@ 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. + α::AbstractVector{T} # diagonal of Jacobi matrix of P + β::AbstractVector{T} # squared supdiagonal of Jacobi matrix of P 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 +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 @@ -20,23 +24,15 @@ _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 +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 -# Should a method be written that makes this more efficient when requesting multiple n? - diff --git a/test/test_monic.jl b/test/test_monic.jl index 8d782b1..deae45d 100644 --- a/test/test_monic.jl +++ b/test/test_monic.jl @@ -1,6 +1,6 @@ using ClassicalOrthogonalPolynomials using Test -using ClassicalOrthogonalPolynomials: Monic, _p0, orthogonalityweight +using ClassicalOrthogonalPolynomials: Monic, _p0, orthogonalityweight, recurrencecoefficients @testset "Basic definition" begin P1 = Legendre() @@ -15,7 +15,7 @@ using ClassicalOrthogonalPolynomials: Monic, _p0, orthogonalityweight @test sprint(show, MIME"text/plain"(), P3) == "Monic(Legendre())" end -@testset "getindex" begin +@testset "evaluation" begin function _pochhammer(x, n) y = one(x) for i in 0:(n-1) @@ -53,6 +53,6 @@ end @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] + # @inferred Q[0.2, 5] # no longer inferred end end \ No newline at end of file From 251e6912b9a399f452f43d8253413a1f0e13bfbd Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Tue, 23 Jul 2024 11:45:09 +0100 Subject: [PATCH 6/7] Doesn't really do much --- src/ClassicalOrthogonalPolynomials.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ClassicalOrthogonalPolynomials.jl b/src/ClassicalOrthogonalPolynomials.jl index 279c247..83e0c97 100644 --- a/src/ClassicalOrthogonalPolynomials.jl +++ b/src/ClassicalOrthogonalPolynomials.jl @@ -326,7 +326,7 @@ function gaussradau(P::Monic{T}, n::Integer, endpt) where {T} α, β = diagonaldata(J), supdiagonaldata(J) endpt = T(endpt) p0, p1 = P[endpt,n:n+1] - a = (endpt - β[end]^2 * p0 / p1)::T + a = endpt - β[end]^2 * p0 / p1 α′ = vcat(@view(α[begin:end-1]), a) J′ = SymTridiagonal(α′, β) x, w = golubwelsch(J′) # not inferred From dcb37b762f7530ab4dad7fc13eb8faa8b657dc39 Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Tue, 23 Jul 2024 11:47:24 +0100 Subject: [PATCH 7/7] Typo --- src/ClassicalOrthogonalPolynomials.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ClassicalOrthogonalPolynomials.jl b/src/ClassicalOrthogonalPolynomials.jl index 83e0c97..57180bf 100644 --- a/src/ClassicalOrthogonalPolynomials.jl +++ b/src/ClassicalOrthogonalPolynomials.jl @@ -346,8 +346,8 @@ function gausslobatto(P::Monic{T}, n::Integer) where {T} a, b = axes(P, 1)[begin], axes(P, 1)[end] J = jacobimatrix(P.P, n + 2) α, β = diagonaldata(J), supdiagonaldata(J) - p0a, p0b = P[a, (n+1):(n+2)] - p1a, p1b = P[b, (n+1):(n+2)] + 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) / Δ