From 5472eb484a478659bad3c8933e654e31cdf36711 Mon Sep 17 00:00:00 2001 From: ioannisPApapadopoulos Date: Thu, 24 Nov 2022 15:19:09 +0000 Subject: [PATCH 1/5] Implement Hilbert transform of LaguerreWeight (and thus all of Laguerre) --- src/classical/laguerre.jl | 3 ++- src/stieltjes.jl | 20 ++++++++++++++++++-- test/test_stieltjes.jl | 9 +++++++++ 3 files changed, 29 insertions(+), 3 deletions(-) diff --git a/src/classical/laguerre.jl b/src/classical/laguerre.jl index b61ed7be..d7b4eb07 100644 --- a/src/classical/laguerre.jl +++ b/src/classical/laguerre.jl @@ -9,7 +9,8 @@ end LaguerreWeight{T}() where T = LaguerreWeight{T}(zero(T)) LaguerreWeight() = LaguerreWeight{Float64}() -axes(::LaguerreWeight{T}) where T = (Inclusion(ℝ),) +# axes(::LaguerreWeight{T}) where T = (Inclusion(ℝ),) +axes(::LaguerreWeight{T}) where T = (Inclusion(HalfLine{T}()),) function getindex(w::LaguerreWeight, x::Number) x ∈ axes(w,1) || throw(BoundsError()) x^w.α * exp(-x) diff --git a/src/stieltjes.jl b/src/stieltjes.jl index b51be050..4bd516d5 100644 --- a/src/stieltjes.jl +++ b/src/stieltjes.jl @@ -94,9 +94,25 @@ end ChebyshevT{T}() * _BandedMatrix(Fill(convert(T,π),1,∞), ℵ₀, 1, -1) end +@simplify function *(H::Hilbert{<:Any,<:HalfLine,<:HalfLine}, w::LaguerreWeight) + T = promote_type(eltype(H), eltype(w)) + α = convert(T, w.α) + z = axes(w,1) + + α == 0 && return real(exp.(log.(expinti.(z) .+ zero(T)*im) .-z)) + + β = -1 < α < 0 ? α : α - ceil(α) + # S = (2* π * im) * (-z)^β*exp(-z)*gamma(-β,-z)/(gamma(-β)*(exp(im*π*β)-exp(-im*π*β))) + S = (2 * convert(T,π) * im) .* exp.(β*log.(-z .+ zero(T)*im) .- z + loggamma.(-β, -z .+ zero(T)*im)) ./ (gamma(-β)*(exp(im*convert(T,π)*β)-exp(-im*convert(T,π)*β))) + if α > 0 + for j in 1:Int(ceil(α)) + S = -gamma(β+j) .+ z .* S + end + end + real(S) +end - -@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wP::Weighted{<:Any,<:OrthogonalPolynomial}) +@simplify function *(H::Hilbert{<:Any,<:Any,<:Any}, wP::Weighted{<:Any,<:OrthogonalPolynomial}) P = wP.P w = orthogonalityweight(P) A = recurrencecoefficients(P)[1] diff --git a/test/test_stieltjes.jl b/test/test_stieltjes.jl index e4afa3fe..c2653f16 100644 --- a/test/test_stieltjes.jl +++ b/test/test_stieltjes.jl @@ -141,6 +141,15 @@ end @test H * ChebyshevTWeight() ≡ QuasiZeros{Float64}((x,)) @test H * ChebyshevUWeight() == π*x @test (H * LegendreWeight())[0.1] ≈ log((0.1+1)/(1-0.1)) + + x = axes(Laguerre(), 1) + @test (inv.(x .- x') * LaguerreWeight())[0.1] ≈ exp(-0.1)*expinti(0.1) + @test isapprox((inv.(x .- x') * LaguerreWeight(-0.5))[0.1], 3.317769, rtol=1e-4) # Mathematica + @test isapprox((inv.(x .- x') * LaguerreWeight(-0.5))[2.1], 1.0829, rtol=1e-4) # Mathematica + @test isapprox((inv.(x .- x') * LaguerreWeight(0.5))[0.1], -1.44067, rtol=1e-4) # Mathematica + @test isapprox((inv.(x .- x') * LaguerreWeight(0.5))[2.1], 0.5017, rtol=1e-3) # Mathematica + @test isapprox((inv.(x .- x') * LaguerreWeight(3.5))[0.1], -3.46658, rtol=1e-4) # Mathematica + @test isapprox((inv.(x .- x') * LaguerreWeight(3.5))[2.1], -5.37663, rtol=1e-4) # Mathematica end @test (Ultraspherical(1) \ (H*wT))[1:10,1:10] == diagm(1 => fill(-π,9)) From 35195f626ef3ebd9fd1ec8e0ddb9495debcd49cf Mon Sep 17 00:00:00 2001 From: ioannisPApapadopoulos Date: Thu, 24 Nov 2022 16:50:29 +0000 Subject: [PATCH 2/5] add conversion operators between (weighted) Laguerre polys --- src/classical/laguerre.jl | 16 ++++++++++++++++ test/test_laguerre.jl | 13 +++++++++++++ 2 files changed, 29 insertions(+) diff --git a/src/classical/laguerre.jl b/src/classical/laguerre.jl index d7b4eb07..846bea84 100644 --- a/src/classical/laguerre.jl +++ b/src/classical/laguerre.jl @@ -65,4 +65,20 @@ recurrencecoefficients(L::Laguerre{T}) where T = ((-one(T)) ./ (1:∞), ((L.α+1 T = promote_type(eltype(D),eltype(L)) D = _BandedMatrix(Fill(-one(T),1,∞), ∞, -1,1) Laguerre(L.α+1)*D +end + +########## +# Conversion +########## + +function \(L::Laguerre, K::Laguerre) + T = promote_type(eltype(L), eltype(K)) + @assert L.α ≈ K.α + one(T) + BandedMatrix(0=>Fill{T}(one(T), ∞), 1=>Fill{T}(-one(T), ∞)) +end + +function \(w_A::Weighted{<:Any,<:Laguerre}, w_B::Weighted{<:Any,<:Laguerre}) + T = promote_type(eltype(w_A), eltype(w_B)) + @assert w_A.P.α + one(T) ≈ w_B.P.α + BandedMatrix(0=>w_B.P.α:∞, -1=>-1:-1:-∞) end \ No newline at end of file diff --git a/test/test_laguerre.jl b/test/test_laguerre.jl index 6017685d..da96134e 100644 --- a/test/test_laguerre.jl +++ b/test/test_laguerre.jl @@ -26,4 +26,17 @@ import ClassicalOrthogonalPolynomials: orthogonalityweight X = L \ (x .* L) @test 0.1*L[0.1,1:10]' ≈ L[0.1,1:11]'*X[1:11,1:10] end + + @testset "Conversion" begin + α = 0.3 + L = Laguerre(α) + P = Laguerre(α-1) + + @test P[0.1, 1:10] ≈ (L[0.1,1:20]'*(L \ P)[1:20,1:20])[1:10] + + wL = Weighted(L) + wP = Weighted(P) + + @test wL[0.1, 1:10] ≈ (wP[0.1,1:20]'*(wP \ wL)[1:20,1:20])[1:10] + end end \ No newline at end of file From e608c5f88994b3fd9861b0ccf933540708cd3149 Mon Sep 17 00:00:00 2001 From: ioannisPApapadopoulos Date: Thu, 24 Nov 2022 17:02:54 +0000 Subject: [PATCH 3/5] improve sig figs. in tests --- test/test_stieltjes.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/test_stieltjes.jl b/test/test_stieltjes.jl index c2653f16..7c248941 100644 --- a/test/test_stieltjes.jl +++ b/test/test_stieltjes.jl @@ -144,12 +144,12 @@ end x = axes(Laguerre(), 1) @test (inv.(x .- x') * LaguerreWeight())[0.1] ≈ exp(-0.1)*expinti(0.1) - @test isapprox((inv.(x .- x') * LaguerreWeight(-0.5))[0.1], 3.317769, rtol=1e-4) # Mathematica - @test isapprox((inv.(x .- x') * LaguerreWeight(-0.5))[2.1], 1.0829, rtol=1e-4) # Mathematica - @test isapprox((inv.(x .- x') * LaguerreWeight(0.5))[0.1], -1.44067, rtol=1e-4) # Mathematica - @test isapprox((inv.(x .- x') * LaguerreWeight(0.5))[2.1], 0.5017, rtol=1e-3) # Mathematica - @test isapprox((inv.(x .- x') * LaguerreWeight(3.5))[0.1], -3.46658, rtol=1e-4) # Mathematica - @test isapprox((inv.(x .- x') * LaguerreWeight(3.5))[2.1], -5.37663, rtol=1e-4) # Mathematica + @test (inv.(x .- x') * LaguerreWeight(-0.5))[0.1] ≈ 3.3177694149902 # Mathematica + @test (inv.(x .- x') * LaguerreWeight(-0.5))[2.1] ≈ 1.08294830124810 # Mathematica + @test (inv.(x .- x') * LaguerreWeight(0.5))[0.1] ≈ -1.44067690942561 # Mathematica + @test (inv.(x .- x') * LaguerreWeight(0.5))[2.1] ≈ 0.501737581714500 # Mathematica + @test (inv.(x .- x') * LaguerreWeight(3.5))[0.1] ≈ -3.46658795669242 # Mathematica + @test (inv.(x .- x') * LaguerreWeight(3.5))[2.1] ≈ -5.37663478253357 # Mathematica end @test (Ultraspherical(1) \ (H*wT))[1:10,1:10] == diagm(1 => fill(-π,9)) From b36ded7d4e73fd3f4650a59b1e1d4d4c0400bce8 Mon Sep 17 00:00:00 2001 From: ioannisPApapadopoulos Date: Fri, 25 Nov 2022 11:34:44 +0000 Subject: [PATCH 4/5] fix tests --- src/classical/laguerre.jl | 18 ++++++++++++++---- test/test_laguerre.jl | 10 ++++++---- 2 files changed, 20 insertions(+), 8 deletions(-) diff --git a/src/classical/laguerre.jl b/src/classical/laguerre.jl index 846bea84..99e084b5 100644 --- a/src/classical/laguerre.jl +++ b/src/classical/laguerre.jl @@ -73,12 +73,22 @@ end function \(L::Laguerre, K::Laguerre) T = promote_type(eltype(L), eltype(K)) - @assert L.α ≈ K.α + one(T) - BandedMatrix(0=>Fill{T}(one(T), ∞), 1=>Fill{T}(-one(T), ∞)) + if L.α ≈ K.α + Eye{T}(∞) + elseif L.α ≈ K.α + 1 + BandedMatrix(0=>Fill{T}(one(T), ∞), 1=>Fill{T}(-one(T), ∞)) + else + error("Not implemented for this choice of L.α and K.α.") + end end function \(w_A::Weighted{<:Any,<:Laguerre}, w_B::Weighted{<:Any,<:Laguerre}) T = promote_type(eltype(w_A), eltype(w_B)) - @assert w_A.P.α + one(T) ≈ w_B.P.α - BandedMatrix(0=>w_B.P.α:∞, -1=>-1:-1:-∞) + if w_A.P.α ≈ w_B.P.α + Eye{T}(∞) + elseif w_A.P.α + 1 ≈ w_B.P.α + BandedMatrix(0=>w_B.P.α:∞, -1=>-one(T):-one(T):-∞) + else + error("Not implemented for this choice of w_A.P.α and w_B.P.α.") + end end \ No newline at end of file diff --git a/test/test_laguerre.jl b/test/test_laguerre.jl index da96134e..f0b47d83 100644 --- a/test/test_laguerre.jl +++ b/test/test_laguerre.jl @@ -31,12 +31,14 @@ import ClassicalOrthogonalPolynomials: orthogonalityweight α = 0.3 L = Laguerre(α) P = Laguerre(α-1) - - @test P[0.1, 1:10] ≈ (L[0.1,1:20]'*(L \ P)[1:20,1:20])[1:10] + + @test L \ L == Eye(∞) + @test P[0.1, 1:10] ≈ (L[0.1,1:10]'*(L \ P)[1:10,1:10])' wL = Weighted(L) wP = Weighted(P) - - @test wL[0.1, 1:10] ≈ (wP[0.1,1:20]'*(wP \ wL)[1:20,1:20])[1:10] + + @test wL \ wL == Eye(∞) + @test wL[0.1, 1:10] ≈ (wP[0.1,1:11]'*(wP \ wL)[1:11,1:10])' end end \ No newline at end of file From 42631c67e7e8b99a9d45f2ab38d443878deac36f Mon Sep 17 00:00:00 2001 From: ioannisPApapadopoulos Date: Fri, 25 Nov 2022 12:01:48 +0000 Subject: [PATCH 5/5] added derivatives for weighted Laguerre --- src/classical/laguerre.jl | 6 ++++++ test/test_laguerre.jl | 9 ++++++++- 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/src/classical/laguerre.jl b/src/classical/laguerre.jl index 99e084b5..65de8af3 100644 --- a/src/classical/laguerre.jl +++ b/src/classical/laguerre.jl @@ -67,6 +67,12 @@ recurrencecoefficients(L::Laguerre{T}) where T = ((-one(T)) ./ (1:∞), ((L.α+1 Laguerre(L.α+1)*D end +@simplify function *(D::Derivative, w_A::Weighted{<:Any,<:Laguerre}) + T = promote_type(eltype(D),eltype(w_A)) + D = BandedMatrix(-1=>one(T):∞) + Weighted(Laguerre{T}(w_A.P.α-1))*D +end + ########## # Conversion ########## diff --git a/test/test_laguerre.jl b/test/test_laguerre.jl index f0b47d83..6352254b 100644 --- a/test/test_laguerre.jl +++ b/test/test_laguerre.jl @@ -1,5 +1,6 @@ -using ClassicalOrthogonalPolynomials, Test +using ClassicalOrthogonalPolynomials, BandedMatrices, Test import ClassicalOrthogonalPolynomials: orthogonalityweight +using ForwardDiff @testset "Laguerre" begin @testset "Laguerre weight" begin @@ -22,6 +23,12 @@ import ClassicalOrthogonalPolynomials: orthogonalityweight L = Laguerre(1/2) @test (D*L)[x,1:4] ≈ [0,-1,x-20/8,-x^2/2 + 7/2*x - 35/8] + α = 0.3 + wL = Weighted(Laguerre(α)) + wP = Weighted(Laguerre(α-1)) + Dₓ = wP \ (D * wL) + ForwardDiff.derivative(x -> (Weighted(Laguerre{eltype(x)}(α)))[x, 1:10],0.1) ≈ (wP[x, 1:11]'*Dₓ[1:11,1:10])' + x = axes(L,1) X = L \ (x .* L) @test 0.1*L[0.1,1:10]' ≈ L[0.1,1:11]'*X[1:11,1:10]