Move fractional calculus from ApproxFun.jl (#61)

* move functions from ApproxFun

* fix imported functions

* allow SpecialFunctions v1

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunSingularities"
 uuid = "f8fcb915-6b99-5be2-b79a-d6dbef8e6e7e"
-version = "0.3.9"
+version = "0.3.10"
 
 [deps]
 ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"
@@ -9,6 +9,7 @@ DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
@@ -20,7 +21,7 @@ Aqua = "0.5"
 DomainSets = "0.4, 0.5, 0.6"
 IntervalSets = "0.5, 0.6, 0.7"
 Reexport = "0.2, 1"
-SpecialFunctions = "2"
+SpecialFunctions = "1, 2"
 StaticArrays = "1"
 julia = "1.6"
 

src/ApproxFunSingularities.jl

@@ -1,5 +1,5 @@
 module ApproxFunSingularities
-using Base, LinearAlgebra, Reexport, IntervalSets, DomainSets, Statistics
+using Base, LinearAlgebra, Reexport, IntervalSets, DomainSets, Statistics, SpecialFunctions
 
 @reexport using ApproxFunBase
 @reexport using ApproxFunOrthogonalPolynomials
@@ -25,7 +25,9 @@ import ApproxFunBase: Fun, SumSpace, SubSpace, WeightSpace, NoSpace,
             coefficients, isconvertible, SpaceOperator, cfstype, mobius, roots,
             splitatroots, domaintype, rangetype, weight, isapproxinteger,
             dotu, components, promoterangespace, ∞, gamma,
-            assert_integer, SpecialEvalPtType, isleftendpoint, isrightendpoint, evaluation_point
+            assert_integer, SpecialEvalPtType, isleftendpoint, isrightendpoint, evaluation_point,
+            @calculus_operator, ConcreteConversion, InterlaceOperator_Diagonal, UnsetSpace,
+            choosedomainspace
 
 import ApproxFunOrthogonalPolynomials: order
 
@@ -52,6 +54,7 @@ include("JacobiWeightOperators.jl")
 include("JacobiWeightChebyshev.jl")
 include("LogWeight.jl")
 include("ExpWeight.jl")
+include("fractionalcalculus.jl")
 
 /(c::Number,f::Fun{<:Ultraspherical}) = c/Fun(f,Chebyshev(domain(f)))
 /(c::Number,f::Fun{<:PolynomialSpace{<:IntervalOrSegment}}) = c/Fun(f,Chebyshev(domain(f)))

src/fractionalcalculus.jl

@@ -0,0 +1,219 @@
+####
+# Support for Fractional derivatives
+#
+# defined as
+####
+
+
+
+# export LeftIntegral,LeftDerivative, RightDerivative, RightIntegral
+
+@calculus_operator LeftIntegral
+@calculus_operator LeftDerivative
+
+@calculus_operator RightIntegral
+@calculus_operator RightDerivative
+
+# DLMF18.17.9 with μ=0.5 and α=β=0
+
+function LeftIntegral(S::Jacobi,k)
+    if S.b==0
+        ConcreteLeftIntegral(S,k)
+    else
+        J=Jacobi(zero(S.a),S.a,domain(S))
+        CLI = ConcreteLeftIntegral(J,k)
+        LeftIntegralWrapper(CLI*Conversion(S,J),k,S,rangespace(CLI))
+    end
+end
+
+function LeftIntegral(S::Ultraspherical,k)
+    J = Jacobi(S)
+    LI = LeftIntegral(J,k)
+    LeftIntegralWrapper(LI*ConcreteConversion(S,J),0.5,S,rangespace(LI))
+end
+
+function LeftIntegral(S::Chebyshev,k)
+    SU = Ultraspherical(Legendre(domain(S)))
+    LI = LeftIntegral(SU,k)
+    LeftIntegralWrapper(LI*ConcreteConversion(S,SU),0.5,S,rangespace(LI))
+end
+
+
+function rangespace(Q::ConcreteLeftIntegral{<:Jacobi{<:IntervalOrSegment}})
+    μ=Q.order
+    S=domainspace(Q)
+
+    JacobiWeight(S.b+μ,zero(S.b),Jacobi(S.b+μ,S.a-μ,domain(S)))
+end
+
+function RightIntegral(S::Jacobi,k)
+    if S.a==0
+        ConcreteRightIntegral(S,k)
+    else
+        J=Jacobi(S.b,zero(S.b),domain(S))
+        CRI = ConcreteRightIntegral(J,k)
+        RightIntegralWrapper(CRI*Conversion(S,J),k,S,rangespace(CRI))
+    end
+end
+
+function rangespace(Q::ConcreteRightIntegral{<:Jacobi{<:IntervalOrSegment}})
+    μ=Q.order
+    S=domainspace(Q)
+    @assert S.a==0
+
+    JacobiWeight(zero(S.a),S.a+μ,Jacobi(S.b-μ,S.a+μ,domain(S)))
+end
+
+for TYP in (:ConcreteLeftIntegral,:ConcreteRightIntegral)
+    @eval begin
+        bandwidths(Q::$TYP{<:Jacobi{<:IntervalOrSegment}}) = (0,0)
+        getindex(Q::$TYP{<:Jacobi{<:IntervalOrSegment},<:Any,T},k::Integer,j::Integer) where {T} =
+            convert(T, jacobi_frac_getindex(domain(Q),zero(T),Q.order,k,j))::T
+    end
+end
+
+
+jacobi_frac_getindex(d::IntervalOrSegment,α,μ,k::Integer,j::Integer) =
+    jacobi_frac_getindex((arclength(d)/2)^μ,α,μ,k,j)
+
+jacobi_frac_getindex(c::Number,α,μ,k::Integer,j::Integer) =
+    k==j ? c*exp(loggamma(α+k)-loggamma(α+μ+k)) : zero(promote_type(typeof(c),typeof(α),typeof(μ)))
+
+
+# jacobi_frac_addentries!(d::IntervalOrSegment,α,μ,A,kr::UnitRange)=
+#     jacobi_frac_addentries!((arclength(d)/2)^μ,α,μ,A,kr)
+# function jacobi_frac_addentries!(c::Number,α,μ,A,kr::UnitRange)
+#     γ=c*gamma(α+1)/gamma(α+1+μ)
+#     for k=1:first(kr)-1
+#         γ*=(α+k)/(α+μ+k)
+#     end
+#     for k=kr
+#        A[k,k]+=γ
+#        γ*=(α+k)/(α+μ+k)
+#         #should be gamma(S.α+k)/gamma(S.α+μ+k)=
+#         #(S.α+k-1)/(S.α+μ+k-1)*gamma(S.α+k-1)/gamma(S.α+μ+k-1)
+#     end
+#     A
+# end
+
+
+function LeftIntegral(S::JacobiWeight{<:Jacobi}, k)
+    J=S.space
+    @assert S.α ≈ 0
+    @assert S.β ≈ J.b
+    ConcreteLeftIntegral(S,k)
+end
+
+function RightIntegral(S::JacobiWeight{<:Jacobi}, k)
+    J=S.space
+    @assert S.α ≈ J.a
+    @assert S.β ≈ 0
+    ConcreteRightIntegral(S,k)
+end
+
+
+
+function LeftIntegral(S::JacobiWeight{<:Chebyshev},k)
+    # convert to Jacobi
+    Q=LeftIntegral(JacobiWeight(S.β,S.α,Jacobi(S.space)),k)
+    LeftIntegralWrapper(Q*Conversion(S,domainspace(Q)),k,S,rangespace(Q))
+end
+
+function RightIntegral(S::JacobiWeight{<:Chebyshev},k)
+    # convert to Jacobi
+    Q=RightIntegral(JacobiWeight(S.β,S.α,Jacobi(S.space)),k)
+    RightIntegralWrapper(Q*Conversion(S,domainspace(Q)),k,S,rangespace(Q))
+end
+
+for (TYP,WRAP) in ((:LeftIntegral,:LeftIntegralWrapper),
+                    (:RightIntegral,:RightIntegralWrapper))
+
+    @eval function $TYP(S::JacobiWeight{<:PolynomialSpace}, k)
+        JS = JacobiWeight(S.β,S.α,Jacobi(S.space))
+        IOP = $TYP(JS,k)
+        $WRAP(IOP*Conversion(S,JS),k,S,rangespace(IOP))
+    end
+end
+
+
+#DLMF18.17.9
+function rangespace(Q::ConcreteLeftIntegral{<:JacobiWeight{<:Jacobi{DD,RR},DD,RR}}) where {DD<:IntervalOrSegment,RR}
+    μ=Q.order
+    S=domainspace(Q)
+    J=S.space
+    if isapprox(S.β,-μ)
+        Jacobi(J.b+μ,J.a-μ,domain(J))
+    else
+        JacobiWeight(S.β+μ,zero(S.β),Jacobi(J.b+μ,J.a-μ,domain(J)))
+    end
+end
+
+function rangespace(Q::ConcreteRightIntegral{<:JacobiWeight{<:Jacobi{DD,RR},DD,RR}}) where {DD<:IntervalOrSegment,RR}
+    μ=Q.order
+    S=domainspace(Q)
+    J=S.space
+    @assert S.α==J.a
+    @assert S.β==0
+    if isapprox(S.β,-μ)
+        Jacobi(J.b-μ,J.a+μ,domain(J))
+    else
+        JacobiWeight(zero(S.α),S.α+μ,Jacobi(J.b-μ,J.a+μ,domain(J)))
+    end
+end
+
+for TYP in (:ConcreteLeftIntegral,:ConcreteRightIntegral)
+    @eval bandwidths(Q::$TYP{<:JacobiWeight{<:Jacobi{DD,RR},DD,RR}}) where {DD<:IntervalOrSegment,RR} = (0,0)
+end
+
+getindex(Q::ConcreteLeftIntegral{<:JacobiWeight{<:Jacobi{DD,RR},DD,RR},<:Any,T},k::Integer,j::Integer) where {DD<:IntervalOrSegment,RR,T} =
+    convert(T,jacobi_frac_getindex(domain(Q),domainspace(Q).β,Q.order,k,j))::T
+
+getindex(Q::ConcreteRightIntegral{<:JacobiWeight{<:Jacobi{DD,RR},DD,RR},<:Any,T},k::Integer,j::Integer) where {DD<:IntervalOrSegment,RR,T} =
+    convert(T,jacobi_frac_getindex(domain(Q),domainspace(Q).α,Q.order,k,j))::T
+
+function choosedomainspace(Q::LeftIntegral{UnsetSpace}, sp::JacobiWeight)
+    @assert sp.β>0 && isapproxinteger(sp.α)
+
+    if isapprox(Q.order,sp.β)
+        Legendre(domain(sp))
+    else
+        JacobiWeight(sp.β-Q.order,0.,Jacobi(sp.β-Q.order,0.,domain(sp)))
+    end
+end
+
+choosedomainspace(Q::LeftIntegral{UnsetSpace,T}, sp::PolynomialSpace) where {T} =
+                JacobiWeight(-Q.order,zero(T),Jacobi(-Q.order,-Q.order,domain(sp)))
+
+
+
+function choosedomainspace(Q::RightIntegral{UnsetSpace,T},sp::JacobiWeight) where T
+    @assert sp.α>0  && isapproxinteger(sp.β)
+
+    if isapprox(Q.order,sp.α)
+        Legendre(domain(sp))
+    else
+        JacobiWeight(zero(T),sp.α-Q.order,Jacobi(0.,sp.α-Q.order))
+    end
+end
+choosedomainspace(Q::RightIntegral{UnsetSpace,T},sp::PolynomialSpace) where {T}=
+    JacobiWeight(zero(T),-Q.order,Jacobi(-Q.order,-Q.order,domain(sp)))
+
+
+
+
+# Define Left/RightDerivative
+for (DTYP,QTYP,DWRAP,QWRAP) in ((:LeftDerivative,:LeftIntegral,:LeftDerivativeWrapper,:LeftIntegralWrapper),
+                            (:RightDerivative,:RightIntegral,:RightDerivativeWrapper,:RightIntegralWrapper))
+    @eval begin
+        function $DTYP(S::Space,k::Real)
+            i=ceil(Int,k)
+            r=i-k
+            $DWRAP(i<0 ? $QTYP(S,-k) : Derivative(i)*$QTYP(S,r),k,S)
+        end
+        function $QTYP(S::SumSpace,k)
+            t = map(s->$QTYP(s,k),components(S))
+            IOP = InterlaceOperator_Diagonal(t,S)
+            $QWRAP(IOP,k)
+        end
+    end
+end