Add Ratio Asymptotics (#4)

* CompatHelper: bump compat for "FillArrays" to "0.11" * Update Project.toml * look at operator asymptotics * Update runtests.jl * ratio asyms * v0.0.2 * Update Project.toml * Update Project.toml * Update runtests.jl * Update runtests.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: Sheehan Olver <[email protected]>
JuliaApproximation · Jan 18, 2021 · 555dc0e · 555dc0e · dlfivefifty · Jan 18, 2021
1 parent 0c7913e
commit 555dc0e
Show file tree

Hide file tree

Showing 4 changed files with 131 additions and 32 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "SemiclassicalOrthogonalPolynomials"
 uuid = "291c01f3-23f6-4eb6-aeb0-063a639b53f2"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.0.1"
+version = "0.0.2"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -17,15 +17,15 @@ QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-ArrayLayouts = "0.4"
-BandedMatrices = "0.15.15"
-ContinuumArrays = "0.3"
-FillArrays = "0.10"
+ArrayLayouts = "0.5.1"
+BandedMatrices = "0.16"
+ContinuumArrays = "0.4"
+FillArrays = "0.11"
 HypergeometricFunctions = "0.3.4"
-InfiniteArrays = "0.8"
-LazyArrays = "0.19"
+InfiniteArrays = "0.8, 0.9"
+LazyArrays = "0.20"
 OrthogonalPolynomialsQuasi = "0.4"
-QuasiArrays = "0.3"
+QuasiArrays = "0.3.8"
 SpecialFunctions = "0.10, 1.0"
 julia = "1.5"
 

diff --git a/src/SemiclassicalOrthogonalPolynomials.jl b/src/SemiclassicalOrthogonalPolynomials.jl
@@ -6,15 +6,17 @@ import Base: getindex, axes, size, \, /, *, +, -, summary, ==, copy, sum, unsafe
 
 import ArrayLayouts: MemoryLayout, ldiv
 import BandedMatrices: bandwidths, _BandedMatrix, AbstractBandedMatrix, BandedLayout
-import LazyArrays: resizedata!, paddeddata, CachedVector, CachedMatrix, LazyMatrix, LazyVector, arguments, ApplyLayout, colsupport
-import OrthogonalPolynomialsQuasi: OrthogonalPolynomial, recurrencecoefficients, jacobimatrix, normalize, recurrencecoefficients, _p0, UnitInterval, orthogonalityweight
+import LazyArrays: resizedata!, paddeddata, CachedVector, CachedMatrix, LazyMatrix, LazyVector, arguments, ApplyLayout, colsupport, AbstractCachedVector
+import OrthogonalPolynomialsQuasi: OrthogonalPolynomial, recurrencecoefficients, jacobimatrix, normalize, _p0, UnitInterval, orthogonalityweight
 import InfiniteArrays: OneToInf, InfUnitRange
 import ContinuumArrays: basis, Weight, @simplify
 import FillArrays: SquareEye
 
-export LanczosPolynomial, Legendre, Normalized, normalize, SemiclassicalJacobi, SemiclassicalJacobiWeight, ConjugateTridiagonal
+export LanczosPolynomial, Legendre, Normalized, normalize, SemiclassicalJacobi, SemiclassicalJacobiWeight, WeightedSemiclassicalJacobi, ConjugateTridiagonal, OrthogonalPolynomialRatio
 
 
+include("ratios.jl")
+
 """"
     SemiclassicalJacobiWeight(t, a, b, c)
 
@@ -63,21 +65,26 @@ struct SemiclassicalJacobi{T,PP<:LanczosPolynomial} <: OrthogonalPolynomial{T}
     c::T
     P::PP # We need to store the basic case where ã,b̃,c̃ = mod(a,-1),mod(b,-1),mod(c,-1)
           # in order to compute lowering operators, etc.
+    SemiclassicalJacobi{T,PP}(t::T,a::T,b::T,c::T,P::PP) where {T,PP<:LanczosPolynomial} = 
+        new{T,PP}(t,a,b,c,P)
 end
 
 const WeightedSemiclassicalJacobi{T} = WeightedBasis{T,<:SemiclassicalJacobiWeight,<:SemiclassicalJacobi}
 
 function SemiclassicalJacobi(t, a, b, c, P::LanczosPolynomial)
-    T = promote_type(typeof(t), typeof(a), typeof(b), typeof(c), eltype(P))
-    SemiclassicalJacobi(T(t), T(a), T(b), T(c), P)
+    T = float(promote_type(typeof(t), typeof(a), typeof(b), typeof(c), eltype(P)))
+    SemiclassicalJacobi{T,typeof(P)}(T(t), T(a), T(b), T(c), P)
 end
 
 SemiclassicalJacobi(t, a, b, c, P::SemiclassicalJacobi) = SemiclassicalJacobi(t, a, b, c, P.P)
 
+WeightedSemiclassicalJacobi(t, a, b, c, P...) = SemiclassicalJacobiWeight(t, a, b, c) .* SemiclassicalJacobi(t, a, b, c, P...)
+
 
 function SemiclassicalJacobi(t, a, b, c)
     ã,b̃,c̃ = mod(a,-1),mod(b,-1),mod(c,-1)
-    P = jacobi(b̃, ã, UnitInterval())
+    T = promote_type(typeof(t), typeof(a), typeof(b), typeof(c))
+    P = jacobi(b̃, ã, UnitInterval{T}())
     x = axes(P,1)
     SemiclassicalJacobi(t, a, b, c, LanczosPolynomial(@.(x^ã * (1-x)^b̃ * (t-x)^c̃), P))
 end
@@ -208,20 +215,14 @@ Gives the Lowering operator from OPs w.r.t. (x-y)*w(x) to Q
 as constructed from Chistoffel–Darboux
 """
 function op_lowering(Q, y)
-    X = jacobimatrix(Q)
-    M = massmatrix(Q)
-    b = X[band(1)] ./ M[band(0)][2:end] .* M[1,1]
-    _BandedMatrix(Vcat((-b .* Q[y,2:∞])', (b .* Q[y,1:∞])'), ∞, 1, 0)
+    R = OrthogonalPolynomialRatio(Q,y)
+    A,_,_ = recurrencecoefficients(Q)
+    # we first use Christoff-Darboux with d = 1
+    # But we want the first OP to be 1 so we rescale
+    d = inv(A[1]*_p0(Q)*R[1])
+    _BandedMatrix(Vcat(Fill(-d,1,∞), (d .* R)'), ∞, 1, 0)
 end
 
-function unsafe_op_lowering(Q, y)
-    X = jacobimatrix(Q)
-    M = massmatrix(Q)
-    b = X[band(1)] ./ M[band(0)][2:end] .* M[1,1]
-    _BandedMatrix(Vcat((-b .* unsafe_getindex(Q,y,2:∞))', (b .* unsafe_getindex(Q,y,1:∞))'), ∞, 1, 0)
-end
-
-
 function semijacobi_ldiv(Q, P::SemiclassicalJacobi)
     if P.a ≤ 0 && P.b ≤ 0 && P.c ≤ 0
         (Q \ P.P)/_p0(P.P)
@@ -268,7 +269,7 @@ function \(w_A::WeightedSemiclassicalJacobi, w_B::WeightedSemiclassicalJacobi)
         @assert A.a == B.a && A.b == B.b && A.c+1 == B.c
         # priority goes to lowering b
         if A.a ≤ 0 && A.b ≤ 0
-            -unsafe_op_lowering(A,A.t)
+            -op_lowering(A,A.t)
         elseif A.b ≤ 0 #lower then raise by inverting
             T = SemiclassicalJacobi(B.t, B.a-1, B.b, B.c-1)
             L = T \ (SemiclassicalJacobiWeight(B.t,1,0,1) .* B)

diff --git a/src/ratios.jl b/src/ratios.jl
@@ -0,0 +1,32 @@
+"""
+    OrthogonalPolynomialRatio(P,x)
+
+is a is equivalent to the vector `P[x,:] ./ P[x,2:end]`` but built from the
+recurrence coefficients of `P`.
+"""
+mutable struct OrthogonalPolynomialRatio{T, PP<:AbstractQuasiMatrix{T}} <: AbstractCachedVector{T}
+    P::PP # OPs
+    x::T
+    data::Vector{T}
+    datasize::Tuple{Int}
+
+    function OrthogonalPolynomialRatio{T, PP}(P::PP, x::T) where {T,PP<:AbstractQuasiMatrix{T}}
+        μ = inv(sqrt(sum(orthogonalityweight(P))))
+        new{T, PP}(P, x, [Base.unsafe_getindex(P,x,1)/Base.unsafe_getindex(P,x,2)], (1,))
+    end
+end
+
+OrthogonalPolynomialRatio(P::AbstractQuasiMatrix{T}, x) where T = OrthogonalPolynomialRatio{T,typeof(P)}(P, convert(T, x))
+
+size(K::OrthogonalPolynomialRatio) = (∞,)
+
+
+function LazyArrays.cache_filldata!(R::OrthogonalPolynomialRatio, inds)
+    A,B,C = recurrencecoefficients(R.P)
+    x = R.x
+    data = R.data
+    @inbounds for n in inds
+        data[n] = inv(A[n]*x + B[n] - C[n] * data[n-1])
+    end
+
+end
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -3,12 +3,22 @@ import BandedMatrices: _BandedMatrix
 import SemiclassicalOrthogonalPolynomials: op_lowering
 import OrthogonalPolynomialsQuasi: recurrencecoefficients, orthogonalityweight
 
+@testset "OrthogonalPolynomialRatio" begin
+    P = Legendre()
+    R = OrthogonalPolynomialRatio(P,0.1)
+    @test P[0.1,1:10] ./ P[0.1,2:11] ≈ R[1:10]
+end
+
 @testset "Jacobi" begin
     P = Normalized(Legendre())
     L = op_lowering(P,1)
     L̃ = P \ WeightedJacobi(1,0)
-    # off by diagonal
+    # off by scaling
     @test (L ./ L̃)[5,5] ≈ (L ./ L̃)[6,5]
+
+    t = 2
+    P̃ = Normalized(SemiclassicalJacobi(t, 0, 0, 0))
+    @test P̃[0.1,1:10] ≈ P[2*0.1-1,1:10]/P[0.1,1]
 end
 
 @testset "SemiclassicalJacobiWeight" begin
@@ -42,13 +52,13 @@ end
             W̃ = LanczosPolynomial(@.(sqrt(x)/sqrt(2-x)), P₊)
             A_W̃, B_W̃, C_W̃ = recurrencecoefficients(W̃)
 
-            @test Normalized(W)[0.1,1:10] ≈ W̃[0.1,1:10] .* (-1) .^ (0:9)
+            @test Normalized(W)[0.1,1:10] ≈ W̃[0.1,1:10]
 
             # this expresses W in terms of W̃
             kᵀ = cumprod(A)
             k_W̃ = cumprod(A_W̃)
             W̄ = (x,n) -> n == 1 ? one(x) : kᵀ[n]*L[n+1,n]/(W̃[0.1,1]k_W̃[n-1]) *W̃[x,n]
-            @test W̄(0.1,5) ≈ W[0.1,5]
+            @test W̄.(0.1,1:5) ≈ W[0.1,1:5]
 
             @testset "x*W == T*L tells us k_n for W" begin
                 @test T[0.1,1] == 1
@@ -91,7 +101,8 @@ end
             @test (T'*(w_T .* T))[1:10,1:10] ≈ sum(w_T)I
             M = W'*(w_W .* W)
             # emperical from mathematica with recurrence
-            @test M[1:3,1:3] ≈ Diagonal([0.5707963267948967,0.5600313808965515,0.5574362259623227])
+            # broken since I changed the scaling
+            @test_broken M[1:3,1:3] ≈ Diagonal([0.5707963267948967,0.5600313808965515,0.5574362259623227])
 
             R = W \ T;
             @test T[0.1,1:10]' ≈ W[0.1,1:10]' * R[1:10,1:10]
@@ -122,7 +133,7 @@ end
             k_Ṽ = cumprod(A_Ṽ)
             V̄ = (x,n) -> n == 1 ? one(x) : -kᵀ[n]*L[n+1,n]/(Ṽ[0.1,1]k_Ṽ[n-1]) *Ṽ[x,n]
 
-            @test V̄(0.1,5) ≈ V[0.1,5]
+            @test V̄.(0.1,1:5) ≈ V[0.1,1:5]
 
             @testset "x*W == T*L tells us k_n for W" begin
                 @test (2-0.1)*V̄(0.1,2) ≈ dot(T[0.1,2:3],L[2:3,2])
@@ -281,4 +292,59 @@ end
         X = Q \ (x .* Q)
         @time X[1:1000,1:1000];
     end
+
+    @testset "BigFloat" begin
+        # T = SemiclassicalJacobi(2, -BigFloat(1)/2, 0, -BigFloat(1)/2)
+    end
 end
+
+@testset "Semiclassical operator asymptotics" begin
+    t = 2
+    P = SemiclassicalJacobi(t, 0, 0, 0)
+    # ratio asymptotics
+    φ = z -> (z + sqrt(z-1)sqrt(z+1))/2
+    U = ChebyshevU()
+
+    @testset "ratio asymptotics" begin
+        n = 200; 
+        @test 2φ(t)*Base.unsafe_getindex(U,t,n)/(Base.unsafe_getindex(U,t,n+1)) ≈ 1
+        @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)
+        @test L1[n+1,n]/L1[n,n] ≈ -1/(2φ(2t-1)) atol=1E-3
+    end
+
+    U = ChebyshevU()
+    n = 20; Base.unsafe_getindex(P,t,n)/Base.unsafe_getindex(U,2t-1,n)
+
+    R_0 = Normalized(SemiclassicalJacobi(t, 1, 0, 0, P)) \ Normalized(P);
+    R_1 = Normalized(SemiclassicalJacobi(t, 0, 1, 0, P)) \ Normalized(P);
+    R_t = Normalized(SemiclassicalJacobi(t, 0, 0, 1, P)) \ Normalized(P);
+
+    @test R_0[999,999:1000] ≈ [0.5,0.5] atol=1e-2
+    @test R_1[999,999:1000] ≈ [0.5,-0.5] atol=1e-2
+    σ = (1 + sqrt(t))/2
+    @test R_t[200,200:201] ≈ [σ,1-σ] atol=1e-2
+
+    n = 100; R_t[n,n:n+1]
+
+    @testset "T,V,W,U" begin
+        T = SemiclassicalJacobi(2, -1/2, 0, -1/2)
+        V = SemiclassicalJacobi(2,  -1/2, 0, 1/2, T)
+        U = SemiclassicalJacobi(2,  1/2, 0, 1/2, T)
+
+        R_t = V \ T;
+        n = 1000
+        c = -1/(2φ(2*2-1))
+        @test R_t[n,n+1]/R_t[n,n] ≈ c atol=1E-3
+        R = U \ T;
+        @test R[n,n+1]/R[n,n] ≈ 1+c atol=1E-3
+        @test R[n,n+2]/R[n,n]  ≈ c atol=1E-3
+        L =T \ (SemiclassicalJacobiWeight(2, 1, 0, 1) .* U)
+        @test 2L[n+1,n] ≈ 1+c atol=1E-3
+        @test 2L[n+2,n] ≈ c atol=1E-3
+
+        # x = z -> 
+    end
+end