Move to ClassicalOrthogonalPolynomials.jl (#12)

* Move to ClassicalOrthogonalPolynomials.jl

* updates

* tests pass

* R asymptotics for P^(1,1,1)

* Update runtests.jl

* v0.0.3

Author: dlfivefifty

Project.toml

@@ -1,31 +1,31 @@
 name = "SemiclassicalOrthogonalPolynomials"
 uuid = "291c01f3-23f6-4eb6-aeb0-063a639b53f2"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.0.2"
+version = "0.0.3"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+ClassicalOrthogonalPolynomials = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 ContinuumArrays = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-OrthogonalPolynomialsQuasi = "aa41a628-2c43-45df-899b-83ab96621781"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 ArrayLayouts = "0.5.1"
 BandedMatrices = "0.16"
-ContinuumArrays = "0.4"
+ClassicalOrthogonalPolynomials = "0.1"
+ContinuumArrays = "0.5"
 FillArrays = "0.11"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.8, 0.9"
 LazyArrays = "0.20"
-OrthogonalPolynomialsQuasi = "0.4"
-QuasiArrays = "0.3.8"
+QuasiArrays = "0.4"
 SpecialFunctions = "0.10, 1.0"
 julia = "1.5"
 

README.md

@@ -7,7 +7,7 @@ A Julia repository for semiclassical orthogonal polynomials
 
 This package implements `SemiclassicalJacobi`, a family of orthogonal 
 polynomials orthogonal with respect to the weight `x^a * (1-x)^b * (t-x)^c`. 
-This builds on top of [OrthogonalPolynomialsQuasi.jl](https://github.com/JuliaApproximation/OrthogonalPolynomialsQuasi.jl) and usage is similar.
+This builds on top of [ClassicalOrthogonalPolynomials.jl](https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl) and usage is similar.
 
 For example, the following gives a half-range Chebyshev polynomial:
 ```julia

src/SemiclassicalOrthogonalPolynomials.jl

@@ -1,13 +1,13 @@
 module SemiclassicalOrthogonalPolynomials
-using OrthogonalPolynomialsQuasi, FillArrays, LazyArrays, ArrayLayouts, QuasiArrays, InfiniteArrays, ContinuumArrays, LinearAlgebra, BandedMatrices, 
+using ClassicalOrthogonalPolynomials, FillArrays, LazyArrays, ArrayLayouts, QuasiArrays, InfiniteArrays, ContinuumArrays, LinearAlgebra, BandedMatrices, 
         SpecialFunctions, HypergeometricFunctions
 
 import Base: getindex, axes, size, \, /, *, +, -, summary, ==, copy, sum, unsafe_getindex
 
 import ArrayLayouts: MemoryLayout, ldiv
 import BandedMatrices: bandwidths, _BandedMatrix, AbstractBandedMatrix, BandedLayout
 import LazyArrays: resizedata!, paddeddata, CachedVector, CachedMatrix, LazyMatrix, LazyVector, arguments, ApplyLayout, colsupport, AbstractCachedVector
-import OrthogonalPolynomialsQuasi: OrthogonalPolynomial, recurrencecoefficients, jacobimatrix, normalize, _p0, UnitInterval, orthogonalityweight
+import ClassicalOrthogonalPolynomials: OrthogonalPolynomial, recurrencecoefficients, jacobimatrix, normalize, _p0, UnitInterval, orthogonalityweight
 import InfiniteArrays: OneToInf, InfUnitRange
 import ContinuumArrays: basis, Weight, @simplify
 import FillArrays: SquareEye
@@ -58,20 +58,19 @@ end
 
 is a quasi-matrix for the  orthogonal polynomials w.r.t. x^a * (1-x)^b * (t-x)^c on [0,1]
 """
-struct SemiclassicalJacobi{T,PP<:LanczosPolynomial} <: OrthogonalPolynomial{T}
+struct SemiclassicalJacobi{T,PP} <: OrthogonalPolynomial{T}
     t::T
     a::T
     b::T
     c::T
     P::PP # We need to store the basic case where ã,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)
+    SemiclassicalJacobi{T,PP}(t::T,a::T,b::T,c::T,P::PP) where {T,PP} = new{T,PP}(t,a,b,c,P)
 end
 
 const WeightedSemiclassicalJacobi{T} = WeightedBasis{T,<:SemiclassicalJacobiWeight,<:SemiclassicalJacobi}
 
-function SemiclassicalJacobi(t, a, b, c, P::LanczosPolynomial)
+function SemiclassicalJacobi(t, a, b, c, P::OrthogonalPolynomial)
     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
@@ -89,6 +88,14 @@ function SemiclassicalJacobi(t, a, b, c)
     SemiclassicalJacobi(t, a, b, c, LanczosPolynomial(@.(x^ã * (1-x)^b̃ * (t-x)^c̃), P))
 end
 
+
+function SemiclassicalJacobi(t, a::Int, b::Int, c::Int)
+    T = promote_type(typeof(t), typeof(a), typeof(b), typeof(c))
+    P = Normalized(legendre(UnitInterval{T}()))
+    x = axes(P,1)
+    SemiclassicalJacobi(t, a, b, c, P)
+end
+
 copy(P::SemiclassicalJacobi) = P
 axes(P::SemiclassicalJacobi) = axes(P.P)
 

test/runtests.jl

@@ -1,7 +1,7 @@
-using SemiclassicalOrthogonalPolynomials, OrthogonalPolynomialsQuasi, ContinuumArrays, BandedMatrices, QuasiArrays, Test, LazyArrays, LinearAlgebra
+using SemiclassicalOrthogonalPolynomials, ClassicalOrthogonalPolynomials, ContinuumArrays, BandedMatrices, QuasiArrays, Test, LazyArrays, LinearAlgebra
 import BandedMatrices: _BandedMatrix
 import SemiclassicalOrthogonalPolynomials: op_lowering
-import OrthogonalPolynomialsQuasi: recurrencecoefficients, orthogonalityweight
+import ClassicalOrthogonalPolynomials: recurrencecoefficients, orthogonalityweight
 
 @testset "OrthogonalPolynomialRatio" begin
     P = Legendre()
@@ -298,6 +298,17 @@ end
     end
 end
 
+@testset "Legendre" begin
+    t = 2
+    P =   SemiclassicalJacobi(t, 0, 0, 0)
+    P¹¹ = SemiclassicalJacobi(t, 1, 1, 0)
+    Q = SemiclassicalJacobi(t, 0, 0, 1)
+    @test P[0.1,1:10] ≈ Normalized(legendre(0..1))[0.1,1:10]
+    @test Normalized(P¹¹)[0.1,1:10] ≈ 2Normalized(jacobi(1,1,0..1))[0.1,1:10]
+    x = axes(P,1)
+    @test LanczosPolynomial(t .- x, legendre(0..1))[0.1,1:10] ≈ Normalized(Q)[0.1,1:10]
+end
+
 @testset "Semiclassical operator asymptotics" begin
     t = 2.2
     P = SemiclassicalJacobi(t, 0, 0, 0)
@@ -337,10 +348,18 @@ end
         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
+        L =T \ (SemiclassicalJacobiWeight(t, 1, 0, 1) .* U)
+        @test L[n+1,n]/L[n,n] ≈ 1+c atol=1E-3
+        @test L[n+2,n]/L[n,n] ≈ c atol=1E-3
+    end
 
-        # x = z -> 
+    @testset "P,Q" begin
+        Q = SemiclassicalJacobi(t, 1, 1, 1, P)
+        R = Q \ P
+        c = -1/(2*φ(2t-1))
+        # (1 + c*z)*(1-z^2) == 1 + c*z - z^2 - c*z^2
+        @test R[200,201]/R[200,200] ≈ c atol=1e-2
+        @test R[200,202]/R[200,200] ≈ -1 atol=1e-2
+        @test R[200,203]/R[200,200] ≈ -c atol=1e-2
     end
 end