Change formula for weighted conversion

src/SemiclassicalOrthogonalPolynomials.jl

@@ -17,7 +17,8 @@ import ClassicalOrthogonalPolynomials: OrthogonalPolynomial, recurrencecoefficie
 import InfiniteArrays: OneToInf, InfUnitRange
 import ContinuumArrays: basis, Weight, @simplify, AbstractBasisLayout, BasisLayout, MappedBasisLayout, grid, plotgrid, equals_layout, ExpansionLayout
 import FillArrays: SquareEye
-import HypergeometricFunctions: _₂F₁general2
+import HypergeometricFunctions: _₂F₁general2, _₂F₁
+import SpecialFunctions: beta
 
 export LanczosPolynomial, Legendre, Normalized, normalize, SemiclassicalJacobi, SemiclassicalJacobiWeight, WeightedSemiclassicalJacobi, OrthogonalPolynomialRatio
 
@@ -129,6 +130,17 @@ SemiclassicalJacobi{T}(t, a, b, c) where T = SemiclassicalJacobi(convert(T,t), c
 
 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...)
 
+# Returns α, β such that P1(x) = β(x-α)
+function _linear_coefficients(t, a, b, c)
+    # beta(a + 1, b + 1) * t^c * _₂F₁(a + 1, -c, a + b + 2, 1/t) is the integral ∫₀¹ wᵗ⁽ᵃᵇᶜ⁾(x) dx 
+    Γᵃ = beta(a + 1, b + 1) * _₂F₁(a + 1, -c, a + b + 2, 1/t)
+    Γᵃ⁺¹ = beta(a + 2, b + 1) * _₂F₁(a + 2, -c, a + b + 3, 1/t)
+    Γᵃ⁺² = beta(a + 3, b + 1) * _₂F₁(a + 3, -c, a + b + 4, 1/t)
+    α = Γᵃ⁺¹/Γᵃ
+    β = sqrt(Γᵃ / (Γᵃ⁺² - 2α*Γᵃ⁺¹ + α^2*Γᵃ))
+    return α, β
+end
+
 function semiclassical_jacobimatrix(t, a, b, c)
     T = float(promote_type(typeof(t), typeof(a), typeof(b), typeof(c)))
     if iszero(a) && iszero(b) && c == -one(T)
@@ -471,13 +483,13 @@ function \(w_A::WeightedSemiclassicalJacobi{T}, w_B::WeightedSemiclassicalJacobi
         wP = Weighted(P) 
         wQ = Weighted(Q)
         R22 = wP \ wQ
-        r11 = A.t - 1 
-        qw0 = SemiclassicalJacobiWeight(Q.t, Q.a, zero(Q.b), Q.c)
-        pw0 = SemiclassicalJacobiWeight(P.t, P.a, zero(P.b), P.c)
-        r21 = wP[:, 1:2] \ (qw0 .- r11 .* pw0)
-        d0 = Vcat(r11, R22[band(0)])
-        d1 = Vcat(r21[begin], R22[band(-1)])
-        d2 = Vcat(r21[begin+1], R22[band(-2)])
+        α, β = _linear_coefficients(P.t, P.a, P.b, P.c)
+        ℓ₁ = A.t - 1 
+        ℓ₂ = 1 + α - A.t 
+        ℓ₃ = inv(β)
+        d0 = Vcat(ℓ₁, R22[band(0)])
+        d1 = Vcat(ℓ₂, R22[band(-1)])
+        d2 = Vcat(ℓ₃, R22[band(-2)])
         data = Hcat(d0, d1, d2)
         return _BandedMatrix(data', 1:∞, 2, 0)
     elseif (wA.a == A.a) && (wA.b == A.b) && (wA.c == A.c) && (wB.a == B.a) && (wB.b == B.b) && (wB.c == B.c) && isinteger(A.a) && isinteger(A.b) && isinteger(A.c) && isinteger(B.a) && isinteger(B.b) && isinteger(B.c)

test/runtests.jl

@@ -673,6 +673,13 @@ end
     end
 end
 
+@testset "_linear_coefficients" begin
+    t, a, b, c = 1.2, 2.3, 0.5, 0.2
+    P = SemiclassicalJacobi(t, a, b, c)
+    α, β = SemiclassicalOrthogonalPolynomials._linear_coefficients(t, a, b, c)
+    @test P[0.2, 2] ≈ β * (0.2 - α)
+end
+
 include("test_derivative.jl")
 include("test_neg1c.jl")
 include("test_neg1b.jl")

test/test_neg1b.jl

@@ -246,6 +246,7 @@ end
     @test Q2.X[1:100, 1:100] ≈ R2.X[1:100, 1:100]
     @test Q3.X[1:100, 1:100] ≈ R3.X[1:100, 1:100]
 end
+
 @testset "Weighted conversion between b=-1" begin
     for (t, a, b, c) in ((2.0, 1 / 2, -1.0, 1 / 2), (2.5, 3 / 2, -1.0, 1 / 2), (2.5, 1.0, -1.0, 2.0))
         Q = SemiclassicalJacobi(t, a, b, c)