Support higher order derivatives (#211)

* Support higher order derivatives

* test not-implemented conversion

Author: dlfivefifty

src/ClassicalOrthogonalPolynomials.jl

@@ -39,7 +39,7 @@ import ContinuumArrays: Basis, Weight, basis_axes, @simplify, Identity, Abstract
     AffineQuasiVector, AffineMap, AbstractWeightLayout, AbstractWeightedBasisLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
     checkpoints, weight, unweighted, MappedBasisLayouts, sum_layout, invmap, plan_ldiv, layout_broadcasted, MappedBasisLayout, SubBasisLayout, broadcastbasis_layout,
     plan_grid_transform, plan_transform, MAX_PLOT_POINTS, MulPlan, grammatrix, AdjointBasisLayout, grammatrix_layout, plan_transform_layout, _cumsum
-import FastTransforms: Λ, ChebyshevGrid, chebyshevpoints, Plan, ScaledPlan, th_cheb2leg
+import FastTransforms: Λ, ChebyshevGrid, chebyshevpoints, Plan, ScaledPlan, th_cheb2leg, pochhammer
 import RecurrenceRelationships: forwardrecurrence, forwardrecurrence!, clenshaw, clenshaw!,
                         check_clenshaw_recurrences
 import RecurrenceRelationshipArrays: initiateforwardrecurrence, Clenshaw

src/classical/chebyshev.jl

@@ -194,17 +194,16 @@ end
 ##########
 
 # Ultraspherical(1)\(D*Chebyshev())
-function diff(S::ChebyshevT{T}; dims=1) where T
+function diff(::ChebyshevT{T}; dims=1) where T
     D = _BandedMatrix((zero(T):∞)', ℵ₀, -1,1)
     ApplyQuasiMatrix(*, ChebyshevU{T}(), D)
 end
 
-function diff(W::Weighted{T,<:ChebyshevU}; dims=1) where T
+function diff(::Weighted{T,<:ChebyshevU}; dims=1) where T
     D =  _BandedMatrix((-one(T):-one(T):(-∞))', ℵ₀, 1,-1)
     ApplyQuasiMatrix(*, Weighted(ChebyshevT{T}()), D)
 end
 
-
 #####
 # Conversion
 #####

src/classical/jacobi.jl

@@ -316,12 +316,18 @@ function _jacobi_convert_b(a, b, k, T) # Jacobi(a, b+k) \ Jacobi(a, b)
     end
 end
 
+isapproxinteger(x) = isinteger(x) || isapprox(x,round(Int,x)) || isapprox(x+1,round(Int,x+1))
+
+
 function \(A::Jacobi, B::Jacobi)
     T = promote_type(eltype(A), eltype(B))
     aa, ab = A.a, A.b
     ba, bb = B.a, B.b
-    ka = Integer(aa-ba)
-    kb = Integer(ab-bb)
+    if !isapproxinteger(aa-ba) || !isapproxinteger(ab-bb)
+        throw(ArgumentError("non-integer conversions not supported"))
+    end
+    ka = round(Integer, aa-ba)
+    kb = round(Integer, ab-bb)
     if ka >= 0
         C1 = _jacobi_convert_a(ba, ab, ka, T)
         if kb >= 0
@@ -419,6 +425,11 @@ end
 # Jacobi(a+1,b+1)\(D*Jacobi(a,b))
 diff(S::Jacobi; dims=1) = ApplyQuasiMatrix(*, Jacobi(S.a+1,S.b+1), _BandedMatrix((((1:∞) .+ (S.a + S.b))/2)', ℵ₀, -1,1))
 
+function diff(S::Jacobi{T}, m::Integer; dims=1) where T
+    D = _BandedMatrix((pochhammer.((S.a + S.b+1):∞, m)/convert(T, 2)^m)', ℵ₀, -m, m)
+    ApplyQuasiMatrix(*, Jacobi(S.a+m,S.b+m), D)
+end
+
 
 #L_6^t
 function diff(WS::HalfWeighted{:a,<:Any,<:Jacobi}; dims=1)

src/classical/ultraspherical.jl

@@ -119,9 +119,12 @@ end
 ##########
 
 # Ultraspherical(1)\(D*Chebyshev())
-diff(S::ChebyshevU; dims=1) = diff(Ultraspherical(S))
+diff(S::ChebyshevU, m...; dims=1) = diff(Ultraspherical(S), m...; dims)
+diff(S::Legendre, m...; dims=1) = diff(Ultraspherical(S), m...; dims)
+
 
 # Ultraspherical(1/2)\(D*Legendre())
+# Special cased as its a Ones
 function diff(S::Legendre{T}; dims=1) where T
     A = _BandedMatrix(Ones{T}(1,∞), ℵ₀, -1,1)
     ApplyQuasiMatrix(*, Ultraspherical{T}(convert(T,3)/2), A)
@@ -134,6 +137,20 @@ function diff(S::Ultraspherical{T}; dims=1) where T
     ApplyQuasiMatrix(*, Ultraspherical{T}(S.λ+1), A)
 end
 
+# higher order 
+
+function diff(::ChebyshevT{T}, m::Integer; dims=1) where T
+    μ = pochhammer(one(T),m-1)*convert(T,2)^(m-1)
+    D = _BandedMatrix((μ * (0:∞))', ℵ₀, -m, m)
+    ApplyQuasiMatrix(*, Ultraspherical{T}(m), D)
+end
+
+function diff(C::Ultraspherical{T}, m::Integer; dims=1) where T
+    μ = pochhammer(convert(T,C.λ),m)*convert(T,2)^m
+    D = _BandedMatrix(Fill(μ,1,∞), ℵ₀, -m, m)
+    ApplyQuasiMatrix(*, Ultraspherical{T}(C.λ+m), D)
+end
+
 # Ultraspherical(λ-1)\ (D*wUltraspherical(λ))
 function diff(WS::Weighted{T,<:Ultraspherical}; dims=1) where T
     S = WS.P

test/test_jacobi.jl

@@ -520,4 +520,21 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
             @test Weighted(Q) \ (w .* Q) isa Diagonal
         end
     end
+
+    @testset "higher order diff" begin
+        P = Jacobi(0.1,0.2)
+        P¹ = Jacobi(1.1,1.2)
+        P² = Jacobi(2.1,2.2)
+        P³ = Jacobi(3.1,3.2)
+        @test (P¹ \ diff(P,1))[1:10,1:10] == (P¹ \ diff(P))[1:10,1:10]
+        @test (P² \ diff(P,2))[1:10,1:10] ≈ (P² \ diff(diff(P)))[1:10,1:10]
+        @test (P³ \ diff(P,3))[1:10,1:10] ≈ (P³ \ diff(diff(diff(P))))[1:10,1:10]
+
+        @test (P² \ diff(P¹,1))[1:10,1:10] ≈ (P² \ diff(P¹))[1:10,1:10]
+        @test (P³ \ diff(P¹,2))[1:10,1:10] ≈ (P³ \ diff(diff(P¹)))[1:10,1:10]
+    end
+
+    @testset "conversion not implemented" begin 
+        @test_throws ArgumentError Jacobi(0,0) \ Jacobi(1.1,2.1)
+    end
 end

test/test_legendre.jl

@@ -208,4 +208,14 @@ import QuasiArrays: MulQuasiArray
         @test @inferred(T \ P[:,1]) == Vcat([1], Zeros(∞))
         @test @inferred(P \ P[:,1]) == Vcat([1], Zeros(∞))
     end
+
+    @testset "higher order diff" begin
+        P = Legendre()
+        P¹ = Ultraspherical(3/2)
+        P² = Ultraspherical(5/2)
+        P³ = Ultraspherical(7/2)
+        @test (P¹ \ diff(P,1))[1:10,1:10] == (P¹ \ diff(P))[1:10,1:10]
+        @test (P² \ diff(P,2))[1:10,1:10] ≈ (P² \ diff(diff(P)))[1:10,1:10]
+        @test (P³ \ diff(P,3))[1:10,1:10] ≈ (P³ \ diff(diff(diff(P))))[1:10,1:10]
+    end
 end

test/test_ultraspherical.jl

@@ -191,4 +191,17 @@ using ClassicalOrthogonalPolynomials: grammatrix
         x = axes(P,1)
         @test sum(@.(ultrasphericalc(1, -1/2, x)/(z-x))) ≈ sum(C[:,2] ./ (z .- x))
     end
+
+    @testset "higher order diff" begin
+        T = ChebyshevT()
+        U = ChebyshevU()
+        C = Ultraspherical(2)
+        C³ = Ultraspherical(3)
+        @test (U \ diff(T,1))[1:10,1:10] == (U \ diff(T))[1:10,1:10]
+        @test (C \ diff(T,2))[1:10,1:10] == (C \ diff(diff(T)))[1:10,1:10]
+        @test (C³ \ diff(T,3))[1:10,1:10] == (C³ \ diff(diff(diff(T))))[1:10,1:10]
+
+        @test (C \ diff(U,1))[1:10,1:10] == (C \ diff(U))[1:10,1:10]
+        @test (C³ \ diff(U,2))[1:10,1:10] == (C³ \ diff(diff(U)))[1:10,1:10]
+    end
 end