Log kernel point evaluation for weighted ChebyshevU (#68)

* Fix bug in colsupport for Hilbert

* add comments

* add comments

* Fix within interval stieltjes evaluation

* Use ChebyshevT

* Update test_stieltjes.jl

* Log kernel real point evalulation

* Notation and Sheehan's fix for complex type

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

Author: ioannisPApapadopoulos

src/stieltjes.jl

@@ -58,6 +58,7 @@ associated(::ChebyshevU{T}) where T = ChebyshevU{T}()
 
 
 const StieltjesPoint{T,W<:Number,V,D} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{BroadcastQuasiMatrix{T,typeof(-),Tuple{W,QuasiAdjoint{V,Inclusion{V,D}}}}}}
+const LogKernelPoint{T<:Real,C,W<:Number,V,D} = BroadcastQuasiMatrix{T,typeof(log),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{BroadcastQuasiMatrix{C,typeof(-),Tuple{W,QuasiAdjoint{V,Inclusion{V,D}}}}}}}}
 const ConvKernel{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(-),Tuple{D1,QuasiAdjoint{T,D2}}}
 const Hilbert{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{ConvKernel{T,Inclusion{T,D1},Inclusion{T,D2}}}}
 const LogKernel{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(log),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{ConvKernel{T,Inclusion{T,D1},Inclusion{T,D2}}}}}}
@@ -79,18 +80,18 @@ end
     log.(x .+ one(T)) .- log.(one(T) .- x)
 end
 
-@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wT::Weighted{<:Any,<:ChebyshevT}) 
+@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wT::Weighted{<:Any,<:ChebyshevT})
     T = promote_type(eltype(H), eltype(wT))
     ChebyshevU{T}() * _BandedMatrix(Fill(-convert(T,π),1,∞), ℵ₀, -1, 1)
 end
 
-@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wU::Weighted{<:Any,<:ChebyshevU}) 
+@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wU::Weighted{<:Any,<:ChebyshevU})
     T = promote_type(eltype(H), eltype(wU))
     ChebyshevT{T}() * _BandedMatrix(Fill(convert(T,π),1,∞), ℵ₀, 1, -1)
 end
 
 
-@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wP::Weighted{<:Any,<:OrthogonalPolynomial}) 
+@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wP::Weighted{<:Any,<:OrthogonalPolynomial})
     P = wP.P
     w = orthogonalityweight(P)
     A = recurrencecoefficients(P)[1]
@@ -150,11 +151,11 @@ end
     PiecewiseInterlace(c,d)  * BlockBroadcastArray{promote_type(eltype(H),eltype(S))}(hvcat, 2, A, B, C, D)
 end
 
-### 
+###
 # LogKernel
 ###
 
-@simplify function *(L::LogKernel{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wT::Weighted{<:Any,<:ChebyshevT}) 
+@simplify function *(L::LogKernel{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wT::Weighted{<:Any,<:ChebyshevT})
     T = promote_type(eltype(L), eltype(wT))
     ChebyshevT{T}() * Diagonal(Vcat(-convert(T,π)*log(2*one(T)),-convert(T,π)./(1:∞)))
 end
@@ -172,11 +173,11 @@ end
 
 
 
-### 
+###
 # PowKernel
 ###
 
-@simplify function *(K::PowKernel, wT::Weighted{<:Any,<:Jacobi}) 
+@simplify function *(K::PowKernel, wT::Weighted{<:Any,<:Jacobi})
     T = promote_type(eltype(K), eltype(wT))
     cnv,α = K.args
     x,y = K.args[1].args[1].args
@@ -219,6 +220,29 @@ sqrtx2(x::Real) = sign(x)*sqrt(x^2-1)
     transpose(convert(T,π) * ξ.^oneto(∞))
 end
 
+####
+# LogKernelPoint
+####
+
+@simplify function *(L::LogKernelPoint, wP::Weighted{<:Any,<:ChebyshevU})
+    T = promote_type(eltype(L), eltype(wP))
+    z, xc = parent(L).args[1].args[1].args
+    if z in axes(wP,1)
+        Tn = Vcat(convert(T,π)*log(2*one(T)), convert(T,π)*ChebyshevT()[z,2:end]./oneto(∞))
+        return transpose((Tn[3:end]-Tn[1:end])/2)
+    else
+        # for U_k where k>=1
+        ξ = inv(z + sqrtx2(z))
+        ζ = (convert(T,π)*ξ.^oneto(∞))./oneto(∞)
+        ζ = (ζ[3:end]- ζ[1:end])/2
+
+        # for U_0
+        ζ = Vcat(convert(T,π)*(ξ^2/4 - (log.(abs.(ξ)) + log(2*one(T)))/2), ζ)
+        return transpose(ζ)
+    end
+
+end
+
 """
    HilbertVandermonde(M, data)
 
@@ -244,7 +268,7 @@ function cache_filldata!(H::HilbertVandermonde{T}, kr, jr) where T
     n,m = H.datasize
     isempty(jr) && return
     resize!(H.colsupport, max(length(H.colsupport), maximum(jr)))
-    
+
     isempty(kr) || (H.data[(n+1):maximum(kr),1:m] .= zero(T))
     for j in (m+1):maximum(jr)
         u = H.M * [H.data[:,j-1]; Zeros{T}(∞)]
@@ -258,7 +282,7 @@ end
     T̃ = chebyshevt(x)
     ψ_1 = T̃ \ inv.(x .+ sqrtx2.(x)) # same ψ_1 = x .- sqrt(x^2 - 1) but with relative accuracy as x -> ∞
     data = convert(eltype(H),π) * Matrix(reshape(paddeddata(ψ_1),:,1))
-    # Operator has columns π * ψ_1^k 
+    # Operator has columns π * ψ_1^k
     T̃ * HilbertVandermonde(Clenshaw(T̃ * ψ_1, T̃), data)
 end
 
@@ -273,14 +297,14 @@ end
     (inv.(z̃ .- x̃') * P)[:,parentindices(wT)[2]]
 end
 
-@simplify function *(H::Hilbert, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}) 
+@simplify function *(H::Hilbert, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
     P = parent(wT)
     x = axes(P,1)
     apply(*, inv.(x .- x'), P)[parentindices(wT)...]
 end
 
 
-@simplify function *(L::LogKernel, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}}) 
+@simplify function *(L::LogKernel, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}})
     V = promote_type(eltype(L), eltype(wT))
     wP = parent(wT)
     kr, jr = parentindices(wT)
@@ -295,7 +319,7 @@ end
 
 ### generic fallback
 for Op in (:Hilbert, :StieltjesPoint, :LogKernel, :PowKernel)
-    @eval @simplify function *(H::$Op, wP::WeightedBasis{<:Any,<:Weight,<:Any}) 
+    @eval @simplify function *(H::$Op, wP::WeightedBasis{<:Any,<:Weight,<:Any})
         w,P = wP.args
         Q = OrthogonalPolynomial(w)
         (H * Weighted(Q)) * (Q \ P)
@@ -497,4 +521,4 @@ function dot(v::AbstractVector{T}, W::PowerLawMatrix, q::AbstractVector{T}) wher
     vpad, qpad = paddeddata(v), paddeddata(q)
     vl, ql = length(vpad), length(qpad)
     return dot(vpad,W[1:vl,1:ql]*qpad)
-end
+end

test/test_interlace.jl

@@ -64,7 +64,7 @@ import ClassicalOrthogonalPolynomials: PiecewiseInterlace
             cμ = H̃[Block.(1:N), Block.(1:N)] \ Vp_cfs_N;
             c1,c2 = cμ[Block(1)]
             μ = W[:,Block.(1:N-1)] * cμ[Block.(2:N)]/2;
-            
+
             # H * μ == Vp(x) + c1 on first interval
             # H * μ == Vp(x) + c2 on second interval
         end

test/test_stieltjes.jl

@@ -1,5 +1,5 @@
 using ClassicalOrthogonalPolynomials, ContinuumArrays, QuasiArrays, BandedMatrices, Test
-import ClassicalOrthogonalPolynomials: Hilbert, StieltjesPoint, ChebyshevInterval, associated, Associated, orthogonalityweight, Weighted, gennormalizedpower, *, dot, PowerLawMatrix, PowKernelPoint
+import ClassicalOrthogonalPolynomials: Hilbert, StieltjesPoint, ChebyshevInterval, associated, Associated, orthogonalityweight, Weighted, gennormalizedpower, *, dot, PowerLawMatrix, PowKernelPoint, LogKernelPoint
 import InfiniteArrays: I
 
 @testset "Associated" begin
@@ -43,6 +43,28 @@ end
         @test (H*w_P)[0.1] ≈ (log(1+(-0.8)) - log(1-(-0.8)))
     end
 
+    @testset "LogKernelPoint" begin
+        wU = Weighted(ChebyshevU())
+        x = axes(wU,1)
+        z = 0.1+0.2im
+        L = log.(abs.(z.-x'))
+        @test L isa LogKernelPoint{Float64,ComplexF64,ComplexF64,Float64,ChebyshevInterval{Float64}}
+
+        @testset "Real point" begin
+            U = ChebyshevU()
+            x = axes(U,1)
+
+            t = 2.0
+            @test (log.(abs.(t .- x') )* Weighted(U))[1,1:3] ≈ [1.0362686329607178,-0.4108206734393296, -0.054364775221816465] #mathematica
+
+            t = 0.5
+            @test (log.(abs.(t .- x') )* Weighted(U))[1,1:3] ≈ [-1.4814921268505252, -1.308996938995747, 0.19634954084936207] #mathematica
+
+            t = 0.5+0im
+            @test (log.(abs.(t .- x') )* Weighted(U))[1,1:3] ≈ [-1.4814921268505252, -1.308996938995747, 0.19634954084936207] #mathematica
+        end
+    end
+
     @testset "Stieltjes" begin
         T = Chebyshev()
         wT = ChebyshevWeight() .* T
@@ -364,4 +386,4 @@ end
             @time ClassicalOrthogonalPolynomials.resizedata!(D,100);
         end
     end
-end
+end