HalfWeighted derivatives (#36)

* Support HalfWeighted Derivatives

* Half weighted a derivative

* Half derivative b

* c half derivative

* Update test_derivative.jl

Author: dlfivefifty

src/derivatives.jl

@@ -44,4 +44,55 @@ end
     Q = wQ.P
     P = SemiclassicalJacobi(Q.t, Q.a-1,Q.b-1,Q.c-1)
     Weighted(P) * ((-sum(orthogonalityweight(Q))/sum(orthogonalityweight(P))) * (Q \ (D * P))')
+end
+
+@simplify function *(D::Derivative, HP::HalfWeighted{:a,<:Any,<:SemiclassicalJacobi})
+    P = HP.P
+    Q = SemiclassicalJacobi(P.t, P.a-1, P.b+1, P.c+1)
+    a = Q.a
+    A,B,C = recurrencecoefficients(P)
+    α,β,γ = recurrencecoefficients(Q)
+    d = AccumulateAbstractVector(*, A ./ α)
+    v1 = MulAddAccumulate(Vcat(0,0,α[2:∞] ./ α), Vcat(0,β))
+    v2 = MulAddAccumulate(Vcat(0,0,A[2:∞] ./ α), Vcat(0,B[1], B[2:end] .* d))
+
+    HalfWeighted{:a}(Q) * _BandedMatrix(
+        Vcat(
+        ((a:∞) .* v2 .- ((a+1):∞) .* Vcat(1,v1[2:end] .* d))',
+        (((a+1):∞) .* Vcat(1,d))'), ℵ₀, 0,1)
+end
+
+@simplify function *(D::Derivative, HP::HalfWeighted{:b,<:Any,<:SemiclassicalJacobi})
+    P = HP.P
+    Q = SemiclassicalJacobi(P.t, P.a+1, P.b-1, P.c+1)
+    b = Q.b
+    A,B,C = recurrencecoefficients(P)
+    α,β,γ = recurrencecoefficients(Q)
+    d = AccumulateAbstractVector(*, A ./ α)
+    d2 = AccumulateAbstractVector(*, A ./ Vcat(1,α))
+    v1 = MulAddAccumulate(Vcat(0,0,α[2:∞] ./ α), Vcat(0,β))
+    v2 = MulAddAccumulate(Vcat(0,0,A[2:∞] ./ α), Vcat(0,B[1], B[2:end] .* d))
+
+    HalfWeighted{:b}(Q) * _BandedMatrix(
+        Vcat(
+        (-(b:∞) .* v2 .+ ((b+1):∞) .* Vcat(1,v1[2:end] .* d) .+ Vcat(0,(1:∞) .* d2))',
+        (-((b+1):∞) .* Vcat(1,d))'), ℵ₀, 0,1)
+end
+
+@simplify function *(D::Derivative, HP::HalfWeighted{:c,<:Any,<:SemiclassicalJacobi})
+    P = HP.P
+    t = P.t
+    Q = SemiclassicalJacobi(t, P.a+1, P.b+1, P.c-1)
+    c = Q.c
+    A,B,C = recurrencecoefficients(P)
+    α,β,γ = recurrencecoefficients(Q)
+    d = AccumulateAbstractVector(*, A ./ α)
+    d2 = AccumulateAbstractVector(*, A ./ Vcat(1,α))
+    v1 = MulAddAccumulate(Vcat(0,0,α[2:∞] ./ α), Vcat(0,β))
+    v2 = MulAddAccumulate(Vcat(0,0,A[2:∞] ./ α), Vcat(0,B[1], B[2:end] .* d))
+
+    HalfWeighted{:c}(Q) * _BandedMatrix(
+        Vcat(
+        (-(c:∞) .* v2 .+ ((c+1):∞) .* Vcat(1,v1[2:end] .* d) .+ Vcat(0,(t:t:∞) .* d2))',
+        (-((c+1):∞) .* Vcat(1,d))'), ℵ₀, 0,1)
 end

test/test_derivative.jl

@@ -146,51 +146,20 @@ import SemiclassicalOrthogonalPolynomials: MulAddAccumulate, HalfWeighted
 
         t,a,b,c = 2,0.1,0.2,0.3
         P = SemiclassicalJacobi(t, a+1, b, c)
-        Q = SemiclassicalJacobi(t,a,b+1,c+1)
+        # Q = SemiclassicalJacobi(t,a,b+1,c+1)
         HP = HalfWeighted{:a}(P)
-        HQ = HalfWeighted{:a}(Q)
+        h = 0.000001
+        @test (D * HP)[0.1,1:10] ≈ (HP[0.1+h,1:10]-HP[0.1,1:10])/h atol=200h
 
-        A,B,C = recurrencecoefficients(P);
-        α,β,γ = recurrencecoefficients(Q);
-
-        k = cumprod(A);
-        κ = cumprod(α);
-        j = Vector{Float64}(undef, 100)
-        j[1] = B[1]
-        for n = 1:length(j)-1
-            j[n+1] = A[n+1]*j[n] + B[n+1]*k[n]
-        end
-        ξ = Vector{Float64}(undef, 100)
-        ξ[1] = β[1]
-        for n = 1:length(ξ)-1
-            ξ[n+1] = α[n+1]*ξ[n] + β[n+1]*κ[n]
-        end
-
-        h = 0.00001
-        n = 1
-        x = 0.1
-        @test (HP[x+h,n]-HP[x,n])/h ≈ (a+1) * HQ[x,n] atol=10h
-        n = 2
-        @test P[x,n] ≈ k[1]*x + j[1]
-        @test Q[x,n] ≈ κ[1]*x + ξ[1]
-        @test HP[x,n] ≈ x^(a+1) * P[x,n]
-        @test (HP[x+h,n]-HP[x,n])/h ≈ x^a * ((a+1)*(k[1]*x + j[1])+ x * k[1]) atol=10h
-        @test (HP[x+h,n]-HP[x,n])/h ≈ x^a * ((a+2)*k[1]*x + (a+1)*j[1]) atol=10h
-        @test (HP[x+h,n]-HP[x,n])/h ≈ x^a * ((a+2)*k[1]*(Q[x,n]-ξ[1])/κ[1] + (a+1)*j[1]) atol=10h
-        @test (HP[x+h,n]-HP[x,n])/h ≈ x^a * ((a+2)*k[1]/κ[1]*Q[x,n]+ (-(a+2)*k[1]*ξ[1]/κ[1] + (a+1)*j[1])*Q[x,n-1]) atol=10h
-        
-        n = 3
-        @test (HP[x+h,n]-HP[x,n])/h ≈ x^a * ((a+n)*k[n-1]/κ[n-1]*Q[x,n]+ (-(a+n)*k[n-1]*ξ[n-1]/κ[n-1] + (a+n-1)*j[n-1])/κ[n-2]*Q[x,n-1]) atol=20h
-        n = 4
-        @test (HP[x+h,n]-HP[x,n])/h ≈ x^a * ((a+n)*k[n-1]/κ[n-1]*Q[x,n]+ (-(a+n)*k[n-1]*ξ[n-1]/κ[n-1] + (a+n-1)*j[n-1])/κ[n-2]*Q[x,n-1]) atol=100h
-
-
-        d = AccumulateAbstractVector(*, A ./ α)
-        # v1 = AccumulateAbstractVector(+, Vcat(1,α) ./ B)
 
-        @test (HP[x+h,n]-HP[x,n])/h ≈ x^a * ((a+n)*d[n-1]*Q[x,n]+ (-(a+n)*k[n-1]*ξ[n-1]/(κ[n-2]κ[n-1]) + (a+n-1)*j[n-1]/κ[n-2])*Q[x,n-1]) atol=100h
+        P = SemiclassicalJacobi(t, a, b+1, c)
+        # Q = SemiclassicalJacobi(t,a+1,b,c+1)
+        HP = HalfWeighted{:b}(P)
+        @test (D * HP)[0.1,1:10] ≈ (HP[0.1+h,1:10]-HP[0.1,1:10])/h atol=1000h
 
-        
-        # Bidiagonal(((P.a+1):∞) .* d, , :U)
+        P = SemiclassicalJacobi(t, a, b, c+1)
+        # Q = SemiclassicalJacobi(t,a+1,b+1,c)
+        HP = HalfWeighted{:c}(P)
+        @test (D * HP)[0.1,1:10] ≈ (HP[0.1+h,1:10]-HP[0.1,1:10])/h atol=2000h
     end
 end