diff --git a/src/classical/laguerre.jl b/src/classical/laguerre.jl index b61ed7be..65de8af3 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) @@ -64,4 +65,36 @@ 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 + +@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 +########## + +function \(L::Laguerre, K::Laguerre) + T = promote_type(eltype(L), eltype(K)) + 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)) + 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 6017685d..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,8 +23,29 @@ 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] end + + @testset "Conversion" begin + α = 0.3 + L = Laguerre(α) + P = Laguerre(α-1) + + @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 \ 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