Use a better method for the derivative of the $b = -1$ polynomials (#122)

* Just use the half-weighted formula

* More efficient

* use latest

* Update ci.yml

* Update Project.toml

Author: DanielVandH

.github/workflows/ci.yml

@@ -20,7 +20,9 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.10'
+          - '1'
+          - 'lts'
+          - 'pre'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

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

src/SemiclassicalOrthogonalPolynomials.jl

@@ -18,7 +18,6 @@ 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 InfiniteLinearAlgebra: BidiagonalConjugation
 
 export LanczosPolynomial, Legendre, Normalized, normalize, SemiclassicalJacobi, SemiclassicalJacobiWeight, WeightedSemiclassicalJacobi, OrthogonalPolynomialRatio
 

src/derivatives.jl

@@ -49,16 +49,13 @@ function diff(P::SemiclassicalJacobi{T}; dims=1) where {T}
         Q = SemiclassicalJacobi(P.t, P.a+1,P.b+1,P.c+1,P)
         Q * divdiff(Q, P)
     elseif P.b == -1
-        Pᵗᵃ⁰ᶜ = SemiclassicalJacobi(P.t, P.a, zero(P.b), P.c)
-        Pᵗᵃ¹ᶜ = SemiclassicalJacobi(P.t, P.a, one(P.b), P.c, Pᵗᵃ⁰ᶜ)
-        Rᵦₐ₁ᵪᵗᵃ⁰ᶜ = Weighted(Pᵗᵃ⁰ᶜ) \ Weighted(Pᵗᵃ¹ᶜ)
-        Dₐ₀ᵪᵃ⁺¹¹ᶜ⁺¹ = diff(Pᵗᵃ⁰ᶜ)
-        Pᵗᵃ⁺¹¹ᶜ⁺¹ = Dₐ₀ᵪᵃ⁺¹¹ᶜ⁺¹.args[1]
-        Pᵗᵃ⁺¹⁰ᶜ⁺¹ = SemiclassicalJacobi(P.t, P.a + 1, zero(P.b), P.c + 1, Pᵗᵃ⁰ᶜ)
-        Rₐ₊₁₀ᵪ₊₁ᵗᵃ⁺¹¹ᶜ⁺¹ = Pᵗᵃ⁺¹¹ᶜ⁺¹ \ Pᵗᵃ⁺¹⁰ᶜ⁺¹
-        Dₐ₋₁ᵪᵃ⁺¹⁰ᶜ⁺¹ = BidiagonalConjugation(Rₐ₊₁₀ᵪ₊₁ᵗᵃ⁺¹¹ᶜ⁺¹, coefficients(Dₐ₀ᵪᵃ⁺¹¹ᶜ⁺¹), Rᵦₐ₁ᵪᵗᵃ⁰ᶜ, 'U')
-        b2 = Vcat(zero(T), zero(T), supdiagonaldata(Dₐ₋₁ᵪᵃ⁺¹⁰ᶜ⁺¹))
-        b1 = Vcat(zero(T), diagonaldata(Dₐ₋₁ᵪᵃ⁺¹⁰ᶜ⁺¹))
+        P1 = SemiclassicalJacobi(P.t, P.a, one(P.b), P.c, P)
+        WP1 = HalfWeighted{:b}(P1)
+        D = diff(WP1)
+        Pᵗᵃ⁺¹⁰ᶜ⁺¹ = D.args[1].P
+        Dmat = D.args[2]
+        b2 = Vcat(zero(T), zero(T), Dmat[band(1)])
+        b1 = Vcat(zero(T), Dmat[band(0)])
         data = Hcat(b2, b1)'
         D = _BandedMatrix(data, ∞, -1, 2)
         return Pᵗᵃ⁺¹⁰ᶜ⁺¹ * D