Large error in two-band transform (#47)

* Large error in two-band transform

* Fix error in transform for associated OPs

* Update Project.toml

* Update equilibriummeasure.jl

* Update test_derivative.jl

* v0.2.2

Author: dlfivefifty

Project.toml

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

examples/equilibriummeasure.jl

@@ -40,7 +40,7 @@ function equilibriumcoefficients(T, b)
     W = Weighted(T)
     t = axes(W,1)
     H = @. inv(t - t')
-    H̃ = U \ (H*W)
+    H̃ = U \H*W
     [1/(b*sum(W[:,1])); 2H̃[:,2:end] \ ( U \ derivative.(V, b*t))]
 end
 function equilibriumendpointvalue(b::Number)
@@ -55,10 +55,24 @@ function equilibrium(b::Number)
     Weighted(U)[affine(-b..b,axes(T,1)),:] * ((Weighted(T) \ Weighted(U))[3:end,:] \ equilibriumcoefficients(T,b)[3:end])
 end
 
+μ = equilibrium(sqrt(2))
+
+T = Chebyshev()
+b = sqrt(2)
+μ = Weighted(T) * equilibriumcoefficients(T, b)
+x = axes(μ,1)
+
+plot(μ)
+
+xx = 0.7; 2b*(log.(abs.(x .- x'))*μ)[xx] - V(b*xx)
+
+
+
 b = 1.0 # initial guess
 for _ = 1:10
     b -= derivative(equilibriumendpointvalue,b) \ equilibriumendpointvalue(b)
 end
+b
 
 plot(equilibrium(b))
 
@@ -92,12 +106,15 @@ function equilibriumcoefficients(P,a,b)
     W = Weighted(P)
     Q = associated(P)
     t = axes(W,1)
-    H = @. inv(t - t')
-    H̃ = Q \ (H*W)
-    [1/(b*sum(W[:,1])); 2H̃[:,2:end] \( Q \ derivative.(V, b*t))]
+    x = axes(Q,1)
+    H = @. inv(x - t')
+    H̃ = Q \ H*W
+    [1/(b*sum(W[:,1])); 2H̃[:,2:end] \( Q \ derivative.(V, b*x))]
 end
 function equilibriumendpointvalues(ab::SVector{2})
     a,b = ab
+    # orthogonal polynomials w.r.t.
+    # abs(x) / (sqrt(1-x^2) * sqrt(x^2 - ρ^2))
     P = TwoBandJacobi(a/b, -one(a)/2, -one(a)/2, one(a)/2)
     Vector(P[[a/b,1],:] * equilibriumcoefficients(P,a,b))
 end
@@ -110,7 +127,11 @@ end
 
 ab = SVector(2.,3.)
 ab -= jacobian(equilibriumendpointvalues,ab) \ equilibriumendpointvalues(ab)
-
+a,b = ab
+xx = range(-4,4;length=1000)
+μ = equilibrium(ab)
+μx = x -> a < abs(x) < b ? μ[x/b] : 0.0
+plot!(xx, μx.(xx))
 
 plot(equilibrium(ab))
 

src/SemiclassicalOrthogonalPolynomials.jl

@@ -11,7 +11,7 @@ import LazyArrays: resizedata!, paddeddata, CachedVector, CachedMatrix, CachedAb
 import ClassicalOrthogonalPolynomials: OrthogonalPolynomial, recurrencecoefficients, jacobimatrix, normalize, _p0, UnitInterval, orthogonalityweight, NormalizedBasisLayout,
                                         Bidiagonal, Tridiagonal, SymTridiagonal, symtridiagonalize, normalizationconstant, LanczosPolynomial,
                                         OrthogonalPolynomialRatio, Weighted, Expansion, UnionDomain, oneto, Hilbert, WeightedOrthogonalPolynomial, HalfWeighted,
-                                        Associated
+                                        Associated, golubwelsch, associated
 import InfiniteArrays: OneToInf, InfUnitRange
 import ContinuumArrays: basis, Weight, @simplify, AbstractBasisLayout, BasisLayout, MappedBasisLayout, grid, plotgrid
 import FillArrays: SquareEye

src/twoband.jl

@@ -1,5 +1,5 @@
-twoband(ρ) = UnionDomain(-1..(-ρ), ρ..1)
-
+twoband(ρ::T) where T = UnionDomain(-one(T)..(-ρ), ρ..one(T))
+twoband_0(ρ::T) where T = UnionDomain(-one(T)..(-ρ), zero(T), ρ..one(T))
 
 """
     TwoBandWeight(ρ, a, b, c)
@@ -32,9 +32,9 @@ function sum(w::TwoBandWeight{T}) where T
     else
         error("not implemented")
 #         1/2 \[Pi] (-((\[Rho]^(1 + 2 b + 2 c)
-#       Gamma[1 + b] Hypergeometric2F1Regularized[-a, 1/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, 
+#    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]]
     end
@@ -69,14 +69,14 @@ copy(A::TwoBandJacobi) = A
 
 orthogonalityweight(Z::TwoBandJacobi) = TwoBandWeight(Z.ρ, Z.a, Z.b, Z.c)
 
-function  getindex(R::TwoBandJacobi, x::Real, j::Integer)
-    ρ = R.ρ
-    if isodd(j)
-        R.P[(1-x^2)/(1-ρ^2), (j+1)÷2]
-    else
-        x * R.Q[(1-x^2)/(1-ρ^2), j÷2]
-    end
-end
+# function  getindex(R::TwoBandJacobi, x::Real, j::Integer)
+#     ρ = R.ρ
+#     if isodd(j)
+#         R.P[(1-x^2)/(1-ρ^2), (j+1)÷2]
+#     else
+#         x * R.Q[(1-x^2)/(1-ρ^2), j÷2]
+#     end
+# end
 
 
 struct Interlace{T,AA,BB} <: LazyVector{T}
@@ -94,28 +94,49 @@ 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)
-    Tridiagonal(Interlace(L.dv/L[1,1], (1-ρ^2) * L.ev), Zeros{T}(∞), Interlace((1-ρ^2) * L.dv,L.ev/L[1,1]))
+    Tridiagonal(Interlace(L.dv/L[1,1], (ρ^2-1) * L.ev), Zeros{T}(∞), Interlace((1-ρ^2) * L.dv,L.ev/(-L[1,1])))
 end
 
 
-@simplify function *(H::Hilbert, w::TwoBandWeight)
-    if 2w.a == 2w.b == -1 && 2w.c == 1
-        zeros(promote_type(eltype(H),eltype(w)), axes(w,1))
+const ConvKernel2{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(-),Tuple{D1,QuasiAdjoint{T,D2}}}
+const Hilbert2{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{ConvKernel2{T,Inclusion{T,D1},Inclusion{T,D2}}}}
+
+@simplify function *(H::Hilbert2, w::TwoBandWeight)
+    if 2w.a == 2w.b == -1 && 2w.c == 1 && axes(H,2) == axes(w,1) && axes(H,1).domain ⊆ twoband_0(w.ρ)
+        zeros(promote_type(eltype(H),eltype(w)), axes(H,1))
     else
        error("Not Implemented")
     end
 end
 
-
-function grid(L::SubQuasiArray{T,2,<:Associated{<:Any,<:TwoBandJacobi},<:Tuple{Inclusion,OneTo}}) where T
-    g = grid(legendre(parent(L).P.ρ .. 1)[:,parentindices(L)[2]])
-    [-g; g]
-end
-    
 function plotgrid(L::SubQuasiArray{T,2,<:TwoBandJacobi,<:Tuple{Inclusion,AbstractUnitRange}}) where T
     g = plotgrid(legendre(parent(L).ρ .. 1)[:,parentindices(L)[2]])
     [-g; g]
 end
 
-LinearAlgebra.factorize(L::SubQuasiArray{<:Any,2,<:Associated{<:Any,<:TwoBandJacobi},<:Tuple{Inclusion,OneTo}}) =
-    ContinuumArrays._factorize(BasisLayout(), L)
+###
+# Associated
+###
+
+axes(Q::Associated{T,<:TwoBandJacobi}) where T = (Inclusion(twoband_0(Q.P.ρ)), axes(Q.P,2))
+function golubwelsch(V::SubQuasiArray{T,2,<:Normalized{<:Any,<:Associated{<:Any,<:TwoBandJacobi}},<:Tuple{Inclusion,AbstractUnitRange}}) where T
+    x,w = golubwelsch(jacobimatrix(V))
+    n = length(x)
+    if isodd(n)
+        x[(n+1)÷2] = zero(T) # make exactly zero
+    else
+        x[n÷2] = zero(T) # make exactly zero
+        x[(n÷2)+1] = zero(T) # make exactly zero
+    end
+    w .*= sum(orthogonalityweight(parent(V)))
+    x,w
+end
+
+@simplify function *(H::Hilbert2, wP::Weighted{<:Any,<:TwoBandJacobi}) 
+    P = wP.P
+    w = orthogonalityweight(P)
+    A = recurrencecoefficients(P)[1]
+    Q = associated(P)
+    @assert axes(H,1) == axes(Q,1)
+    Q * BandedMatrix(1 =>Fill(-A[1]*sum(w),∞))
+end

test/test_derivative.jl

@@ -35,7 +35,7 @@ import SemiclassicalOrthogonalPolynomials: MulAddAccumulate, HalfWeighted
         L = Q \ (D * P̃);
         # L is bidiagonal
         @test norm(triu(L[1:10,1:10],3)) ≤ 3E-12
-        @test L[:,5] isa Vcat
+        @test L[:,5][1:10] ≈ L[1:10,5]
 
         A,B,C = recurrencecoefficients(P);
         α,β,γ = recurrencecoefficients(Q);

test/test_twoband.jl

@@ -1,4 +1,4 @@
-using SemiclassicalOrthogonalPolynomials, ClassicalOrthogonalPolynomials, Test
+using SemiclassicalOrthogonalPolynomials, ClassicalOrthogonalPolynomials, BandedMatrices, Test
 import ClassicalOrthogonalPolynomials: orthogonalityweight, Weighted, associated, plotgrid
 
 @testset "Two Band" begin
@@ -27,21 +27,33 @@ import ClassicalOrthogonalPolynomials: orthogonalityweight, Weighted, associated
         @test !issymmetric(jacobimatrix(T)[1:10,1:10])
     end
 
+    @testset "Associated" begin
+        ρ = 0.5
+        T = TwoBandJacobi(ρ, -1/2, -1/2, 1/2)
+        Q = associated(T)
+        x = axes(Q,1)
+        @test 0 in x
+        @test (Q[0.6,1:101]' * (Q[:,Base.OneTo(101)] \ exp.(x))) ≈ exp(0.6)
+        @test (Q[-0.6,1:101]' * (Q[:,Base.OneTo(101)] \ exp.(x))) ≈ exp(-0.6)
+
+        @test (Q * (Q \ exp.(x)))[0.6] ≈ exp(0.6)
+    end
+
     @testset "Hilbert" begin
         ρ = 0.5
         w = TwoBandWeight(ρ, -1/2, -1/2, 1/2)
-        x = axes(w,1)
-        H = inv.(x .- x')
+        T = TwoBandJacobi(ρ, -1/2, -1/2, 1/2)
+        Q = associated(T)
+        t = axes(w,1)
+        x = axes(Q,1)
+        H = inv.(x .- t')
         @test iszero(H*w)
         @test sum(w) ≈ π
 
-        T = TwoBandJacobi(ρ, -1/2, -1/2, 1/2)
-        Q = associated(T)
-        @test (Q[0.6,1:100]' * (Q[:,Base.OneTo(100)] \ exp.(x))) ≈ exp(0.6)
-        @test (Q[-0.6,1:100]' * (Q[:,Base.OneTo(100)] \ exp.(x))) ≈ exp(-0.6)
+        B = Q \ H * Weighted(T)
+        @test B isa BandedMatrix
 
-        @test (Q * (Q \ exp.(x)))[0.6] ≈ exp(0.6)
-        @test_broken Q \ (H * Weighted(T)) # need to deal with Hcat
+        @test Q[0.6,:]'* (B * (T\exp.(t))) ≈ -4.701116657130821 # Mathematica
 
         @test_broken H*TwoBandWeight(ρ, 1/2, 1/2, -1/2)
     end
@@ -51,4 +63,12 @@ import ClassicalOrthogonalPolynomials: orthogonalityweight, Weighted, associated
         P = TwoBandJacobi(ρ, -1/2, -1/2, 1/2)
         @test all(x -> ρ ≤ abs(x) ≤ 1, plotgrid(P[:,1:5]))
     end
+
+    @testset "associated transform error" begin
+        a,b = 0.5, 0.9
+        Vp = x -> 4x^3 - 20x
+        Q = associated(TwoBandJacobi((a/b), -1/2, -1/2, 1/2))
+        x = axes(Q,1)
+        @test norm((Q[:,Base.OneTo(30)] \ Vp.(b*x))[1]) ≤ 1E-13
+    end
 end