Jp/twoband (#64)

* Fix jacobimatrix (plus a test)

* Add support for derivative of bubble functions

* Fix tests on definition of TwoBand + add derivative and inner product of HalfWeighted{:ab}(TwoBand(\rho, 1,1,0))

* Type-stability and simplifications

* Fix first failing test

* Update test_twoband.jl

* v0.3.2

Co-authored-by: ioannisPApapadopoulos <[email protected]>
Co-authored-by: Sheehan Olver <[email protected]>

Author: 3 people

Project.toml

@@ -1,7 +1,7 @@
 name = "SemiclassicalOrthogonalPolynomials"
 uuid = "291c01f3-23f6-4eb6-aeb0-063a639b53f2"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/SemiclassicalOrthogonalPolynomials.jl

@@ -47,7 +47,8 @@ function getindex(P::SemiclassicalJacobiWeight, x::Real)
 end
 
 function sum(P::SemiclassicalJacobiWeight{T}) where T
-    (t,a,b,c) = map(big, map(float, (P.t,P.a,P.b,P.c)))
+    # (t,a,b,c) = map(big, map(float, (P.t,P.a,P.b,P.c)))
+    (t,a,b,c) = P.t, P.a, P.b, P.c
     return convert(T, t^c * beta(1+a,1+b) * _₂F₁general2(1+a,-c,2+a+b,1/t))
 end
 
@@ -166,7 +167,7 @@ SemiclassicalJacobi(t, a, b, c, P::SemiclassicalJacobi) = SemiclassicalJacobi(t,
 SemiclassicalJacobi(t, a, b, c) = SemiclassicalJacobi(t, a, b, c, semiclassical_jacobimatrix(t, a, b, c))
 SemiclassicalJacobi{T}(t, a, b, c) where T = SemiclassicalJacobi(convert(T,t), convert(T,a), convert(T,b), convert(T,c))
 
-WeightedSemiclassicalJacobi(t, a, b, c, P...) = SemiclassicalJacobiWeight(t, a, b, c) .* SemiclassicalJacobi(t, a, b, c, P...)
+WeightedSemiclassicalJacobi{T}(t, a, b, c, P...) where T = SemiclassicalJacobiWeight{T}(convert(T,t), convert(T,a), convert(T,b), convert(T,c)) .* SemiclassicalJacobi{T}(convert(T,t), convert(T,a), convert(T,b), convert(T,c), P...)
 
 
 
@@ -360,6 +361,7 @@ end
 
 \(A::SemiclassicalJacobi, w_B::WeightedSemiclassicalJacobi) = (SemiclassicalJacobiWeight(A.t,0,0,0) .* A) \ w_B
 
+\(w_A::Weighted{<:Any,<:SemiclassicalJacobi}, w_B::Weighted{<:Any,<:SemiclassicalJacobi}) = convert(WeightedBasis, w_A) \ convert(WeightedBasis, w_B)
 
 massmatrix(P::SemiclassicalJacobi) = Diagonal(Fill(sum(orthogonalityweight(P)),∞))
 
@@ -416,6 +418,7 @@ convert(::Type{WeightedBasis}, Q::HalfWeighted{:bc,T,<:SemiclassicalJacobi}) whe
 convert(::Type{WeightedBasis}, Q::HalfWeighted{:ac,T,<:SemiclassicalJacobi}) where T = SemiclassicalJacobiWeight(Q.P.t, Q.P.a,zero(T),Q.P.c) .* Q.P
 
 
+include("twoband.jl")
 include("derivatives.jl")
 
 
@@ -473,7 +476,6 @@ function LazyArrays.cache_filldata!(P::SemiclassicalJacobiFamily, inds::Abstract
     P
 end
 
-include("twoband.jl")
 include("lowering.jl")
 
 end

src/derivatives.jl

@@ -215,5 +215,4 @@ end
     t = P.t
     HQ = HalfWeighted{:bc}(SemiclassicalJacobi(t, P.a+1,P.b-1,P.c-1))
     HQ * divmul(HQ, D, HP)
-end
-
+end

src/twoband.jl

@@ -30,16 +30,15 @@ function sum(w::TwoBandWeight{T}) where T
     if 2w.a == 2w.b == -1 && 2w.c == 1
         convert(T,π)
     else
-        error("not implemented")
-#         1/2 \[Pi] (-((\[Rho]^(1 + 2 b + 2 c)
-#       Gamma[1 + b] Hypergeometric2F1Regularized[-a, 1/2 + c,
-#       3/2 + b + c, \[Rho]^2])/Gamma[1/2 - c]) + (
-#    Gamma[1 + a] Hypergeometric2F1Regularized[-b, -(1/2) - a - b - c,
-#      1/2 - b - c, \[Rho]^2])/
-#    Gamma[3/2 + a + b + c]) Sec[(b + c) \[Pi]]
+        a,b,c,ρ = w.a,w.b,w.c,w.ρ
+        convert(T,π)/2 * sec((b + c)*convert(T,π)) * (
+            - ρ^(one(T) + 2b + 2c) * gamma(one(T) + b) * _₂F₁(-a, one(T)/2 + c, convert(T,3)/2 + b + c, ρ^2) / gamma(one(T)/2 - c)
+            + gamma(one(T) + a) * _₂F₁(-b, -one(T)/2 - a - b - c, one(T)/2 - b - c, ρ^2) / gamma(convert(T,3)/2 + a + b + c)
+        )
     end
 end
 
+
 """
     TwoBandJacobi(ρ, a, b, c)
 
@@ -78,6 +77,8 @@ orthogonalityweight(Z::TwoBandJacobi) = TwoBandWeight(Z.ρ, Z.a, Z.b, Z.c)
 #     end
 # end
 
+weight(W::HalfWeighted{:ab,T,<:TwoBandJacobi}) where T = TwoBandWeight(W.P.ρ, W.P.a,W.P.b,zero(T))
+convert(::Type{WeightedBasis}, Q::HalfWeighted{:ab,T,<:TwoBandJacobi}) where T = TwoBandWeight(Q.P.ρ, Q.P.a,Q.P.b,zero(T)) .* Q.P
 
 struct Interlace{T,AA,BB} <: LazyVector{T}
     a::AA
@@ -94,7 +95,10 @@ getindex(A::Interlace{T}, k::Int) where T = convert(T, isodd(k) ? A.a[(k+1)÷2]
 function jacobimatrix(R::TwoBandJacobi{T}) where T
     ρ = R.ρ
     L = R.P \ (SemiclassicalJacobiWeight(R.P.t,0,0,1) .* R.Q)
+    # M = (L / L[1,1])' # equal to R.Q \ R.P
+
     Tridiagonal(Interlace(L.dv/L[1,1], (ρ^2-1) * L.ev), Zeros{T}(∞), Interlace((1-ρ^2) * L.dv,L.ev/(-L[1,1])))
+    # Tridiagonal(Interlace(L.dv/L[1,1], (1-ρ^2) * L.ev), Zeros{T}(∞), Interlace((1-ρ^2) * L.dv, L.ev/L[1,1]))
 end
 
 
@@ -140,3 +144,69 @@ end
     @assert axes(H,1) == axes(Q,1)
     Q * BandedMatrix(1 =>Fill(-A[1]*sum(w),∞))
 end
+
+
+##
+# Derivative of double weighted TwoBandJacobi
+##
+function divmul(R::TwoBandJacobi, D::Derivative, HP::HalfWeighted{:ab,<:Any,<:TwoBandJacobi})
+    T = promote_type(eltype(R), eltype(HP))
+    ρ=convert(T,R.ρ); t=inv(one(T)-ρ^2)
+    a,b,c = R.a,R.b,R.c
+
+    Dₑ = -2*(one(T)-ρ^2) .* (R.Q \ (Derivative(axes(R.Q,1))*HalfWeighted{:ab}(HP.P.P)))
+    D₀ = -2*(one(T)-ρ^2)^2 .* (Weighted(R.P) \ (Derivative(axes(R.P,1))*Weighted(HP.P.Q)))
+
+    (dₑ, dlₑ, d₀, dl₀) = Dₑ.data[1,:], Dₑ.data[2,:], D₀.data[1,:], D₀.data[2,:]
+    BandedMatrix(-1=>Interlace(dₑ, -d₀), -3=>Interlace(-dlₑ, dl₀))
+end
+
+@simplify function *(D::Derivative, HP::HalfWeighted{:ab,<:Any,<:TwoBandJacobi})
+    P = HP.P
+    ρ = P.ρ
+    @assert P.a == 1 && P.b == 1 && P.c == 0
+
+    R = TwoBandJacobi(ρ, P.a-1,P.b-1,P.c)
+    R * divmul(R, D, HP)
+end
+
+###
+# L^2 inner product of double weighted TwoBandJacobi
+###
+@simplify function *(A::QuasiAdjoint{<:Any,<:HalfWeighted{:ab,<:Any,<:TwoBandJacobi}}, B::HalfWeighted{:ab,<:Any,<:TwoBandJacobi})
+    T = promote_type(eltype(A), eltype(B))
+    P = B.P
+    a,b,c = P.a,P.b,P.c
+    
+    @assert A.parent.P.a == a && A.parent.P.b == b && A.parent.P.c == c
+    @assert a == 1 && b == 1 && c == 0
+
+    ρ = convert(T,P.ρ)
+    t = inv(one(T)-ρ^2)
+    Pₗ = SemiclassicalJacobi{T}(t,a-one(T),b-one(T),c-one(T)/2)
+    Qₗ = SemiclassicalJacobi{T}(t,a-one(T),b-one(T),c+one(T)/2)
+
+    Lₑ = Pₗ \ HalfWeighted{:ab}(P.P)
+    Lₒ = Weighted(Qₗ) \ Weighted(P.Q)
+    
+    mₑ = (one(T)-ρ^2)^4*(1-ρ^2)^(one(T)/2)*sum(orthogonalityweight(Pₗ))
+    mₒ = (one(T)-ρ^2)^4*(1-ρ^2)^(convert(T,3)/2)*sum(orthogonalityweight(Qₗ))
+
+    mₑ = Fill{T}(mₑ, ∞);  mₒ = Fill{T}(mₒ, ∞)
+
+    # Sum of entries in each column squared.
+    dₑ = mₑ.*((Lₑ .* Lₑ)' * Ones{T}(∞))
+    d₀ = mₒ.*((Lₒ .* Lₒ)' * Ones{T}(∞))
+
+    # Sum of entries x entries in next column.
+    dlₑ = mₑ.* (Lₑ[2:∞,:] .* Lₑ[2:∞,2:∞])' * Ones{T}(∞)
+    dl₀ = mₒ.* (Lₒ[2:∞,:] .* Lₒ[2:∞,2:∞])' * Ones{T}(∞)
+
+    # Sum of entries and entries in second column over.
+    dllₑ = mₑ.* (Lₑ[3:∞,:] .* Lₑ[3:∞,3:∞])' * Ones{T}(∞)
+    dll₀ = mₒ.* (Lₒ[3:∞,:] .* Lₒ[3:∞,3:∞])' * Ones{T}(∞)
+    
+    BandedMatrix(0=>Interlace(dₑ, d₀), 
+        -2=>Interlace(-dlₑ, -dl₀), 2=>Interlace(-dlₑ, -dl₀), 
+        -4=>Interlace(dllₑ, dll₀), 4=>Interlace(dllₑ, dll₀))
+end

test/runtests.jl

@@ -365,7 +365,7 @@ end
         @test 2φ(2t-1)*Base.unsafe_getindex(P,t,n)/(Base.unsafe_getindex(P,t,n+1)) ≈ 1 atol=1E-3
 
 
-        L1 = P \ WeightedSemiclassicalJacobi(t,0,0,1,P)
+        L1 = P \ Weighted(SemiclassicalJacobi(t,0,0,1,P))
         @test L1[n+1,n]/L1[n,n] ≈ -1/(2φ(2t-1)) atol=1E-3
     end
 

test/test_twoband.jl

@@ -1,5 +1,6 @@
-using SemiclassicalOrthogonalPolynomials, ClassicalOrthogonalPolynomials, BandedMatrices, Test
+using SemiclassicalOrthogonalPolynomials, ClassicalOrthogonalPolynomials, BandedMatrices, LinearAlgebra, Test
 import ClassicalOrthogonalPolynomials: orthogonalityweight, Weighted, associated, plotgrid
+import SemiclassicalOrthogonalPolynomials: Interlace, HalfWeighted, WeightedBasis
 
 @testset "Two Band" begin
     @testset "TwoBandWeight" begin
@@ -23,6 +24,10 @@ import ClassicalOrthogonalPolynomials: orthogonalityweight, Weighted, associated
         @test U ≠ T
         @test orthogonalityweight(T) == TwoBandWeight(ρ, -1/2, -1/2, 1/2)
 
+        R = TwoBandJacobi(ρ, 1, 1, 0)
+        x=0.6; τ=(1-x^2)*R.P.t; n=10
+        @test R[x, 1:n] ≈ Interlace((-1).^(2:n÷2+1).*R.P[τ, 1:n÷2], (-1).^(2:n÷2+1).*x.*R.Q[τ, 1:n÷2])[1:n]
+
         # bug
         @test !issymmetric(jacobimatrix(T)[1:10,1:10])
     end
@@ -71,4 +76,39 @@ import ClassicalOrthogonalPolynomials: orthogonalityweight, Weighted, associated
         x = axes(Q,1)
         @test norm((Q[:,Base.OneTo(30)] \ Vp.(b*x))[1]) ≤ 1E-13
     end
+
+    @testset "derivative half weighted twoband" begin    
+        ρ = 0.2
+        R = TwoBandJacobi(ρ, 0, 0, 0)
+        D = R \ Derivative(axes(R,1))*HalfWeighted{:ab}(TwoBandJacobi(ρ, 1, 1, 0))
+
+        @test (D.l, D.u) == (3, -1)
+        @test D[1:5,1] ≈ [0.0, -0.33858064516128994, 0.0, -1.030932383768154, 0.0]
+        @test D[1:5,2] ≈ [0.0, 0.0, -0.4609776443497252, 0.0, -0.34580926348402935]
+    end
+
+    @testset "L2 inner product" begin
+        ρ = 0.2
+        HP = HalfWeighted{:ab}(TwoBandJacobi(ρ, 1, 1, 0))
+        M = HP' * HP
+        
+        @test (M.l, M.u) == (4, 4)
+
+        Mm = M[1:20, 1:20]
+        @test Mm == Mm'
+
+        # The following is checked via QuadGK, e.g. 2*quadgk(x->HP[x,n]*HP[x,n], ρ, 1, rtol=1e-3)[1]
+        @test diag(Mm[1:5, 1:5]) ≈ [0.0399723828825396, 0.01926644161200575, 0.028656015290810938, 0.014202710309357196, 0.02719735154820007]
+        @test diag(Mm[2:6, 1:5]) ≈ [0, 0, 0, 0, 0]
+        @test diag(Mm[3:7, 1:5]) ≈ [0.003417521501385307, -0.0013645130261003458, 0.0008788935542028268, -0.0005306108105399927, 0.0003587270849050795]
+        @test diag(Mm[4:8, 1:5]) ≈ [0, 0, 0, 0, 0]
+        @test diag(Mm[5:9, 1:5]) ≈ [-0.011436406405959944, -0.004811623256062165, -0.012192463544540394, -0.005436023256389608, -0.012469708623429861]
+    end
+
+    @testset "halfweight" begin
+        ρ = 0.5
+        T = TwoBandJacobi(ρ, -1/2, -1/2, 1/2)
+        @test HalfWeighted{:ab}(T)[0.6,1:10] ≈ convert(WeightedBasis, HalfWeighted{:ab}(T))[0.6,1:10] ≈ (1-0.6^2)^(-1/2) * (0.6^2-ρ^2)^(-1/2) * T[0.6,1:10]
+        
+    end
 end