Skip to content

Commit

Permalink
Make Caputo and RL type stable
Browse files Browse the repository at this point in the history 
Signed-off-by: ErikQQY <[email protected]>
  • Loading branch information
@ErikQQY
ErikQQY committed Mar 23, 2024
1 parent 679f1e5 commit e4f1c5f
Show file tree
Hide file tree
Showing 2 changed files with 71 additions and 71 deletions.
70 changes: 35 additions & 35 deletions src/Derivative/Caputo.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -190,11 +190,11 @@ struct CaputoL2C <: Caputo end


function fracdiff(f::Union{Function, Number},
α::Float64,
α::T,
start_point::Number,
end_point::Real,
h::Float64,
::CaputoDirect) :: Float64
h::T,
::CaputoDirect) where {T <: Real}
#checks(f, α, start_point, end_point)
isa(f, Real) ? (end_point == 0 ? 0 : f/sqrt(pi*end_point)) : nothing
#Using complex step differentiation to calculate the first order differentiation
Expand All @@ -207,33 +207,33 @@ end
fracdiff(f, α, end_point, h, ::CaputoDirect) = fracdiff(f, α, 0, end_point, h, CaputoDirect())

function fracdiff(f::FunctionAndNumber,
α::Float64,
α::T,
start_point::Real,
end_point::AbstractArray,
h::Float64,
::CaputoDirect)::Vector
h::T,
::CaputoDirect)::Vector where {T <: Real}
result = map(x->fracdiff(f, α, start_point, x, h, CaputoDirect())[1], end_point)
return result
end

function fracdiff(f::Union{Function, Number},
α::Float64,
α::T,
end_point::Real,
h::Float64,
::CaputoDirect)
h::T,
::CaputoDirect) where {T <: Real}
fracdiff(f, α, 0, end_point, h, CaputoDirect())
end

function fracdiff(f::FunctionAndNumber,
α::Float64,
α::T,
end_point::Real,
h::Float64,
::CaputoTrap)
h::T,
::CaputoTrap) where {T <: Real}
isa(f, Number) ? (end_point == 0 ? (return 0) : (return f/sqrt(pi*end_point))) : nothing
end_point == 0 ? (return 0) : nothing
m = floor(Int, α)+1

summation = zero(Float64)
summation = zero(T)
n = floor(Int, end_point/h)

@fastmath @inbounds @simd for i in 0:n
Expand Down Expand Up @@ -272,13 +272,13 @@ Caputo Diethelm algorithm
=#
function fracdiff(f::FunctionAndNumber,
α::T,
end_point::T,
end_point::P,
h::T,
::CaputoDiethelm) where {T <: Real}
::CaputoDiethelm) where {T <: Real, P <: Number}
#checks(f, α, 0, end_point)
isa(f, Real) ? (end_point == 0 ? 0 : f/sqrt(pi*end_point)) : nothing
N = round(Int, end_point/h)
result = zero(Float64)
result = zero(T)

@fastmath @inbounds @simd for i in 0:N
result += quadweights(i, N, α)*(f(end_point-i*h) - f(0))
Expand All @@ -299,9 +299,9 @@ end

function fracdiff(f::FunctionAndNumber,
α::T,
end_point::AbstractArray{T},
end_point::AbstractArray{P},
h::T,
::CaputoDiethelm) where {T <: Real}
::CaputoDiethelm) where {T <: Real, P <: Number}
result = map(x->fracdiff(f, α, x, h, CaputoDiethelm()), end_point)
return result
end
Expand Down Expand Up @@ -342,7 +342,7 @@ function fracdiff(f::FunctionAndNumber,
h::Float64,
p::Integer,
::CaputoHighOrder)
n::Int64 = round(Int, T/h)
n = round(Int, T/h)
A = wj(n, p, α)
B = fj(n, T, f)
C = ones(1, n)
Expand All @@ -351,8 +351,8 @@ function fracdiff(f::FunctionAndNumber,
return D[1]
end

function wj(n::Int64, r::Int64, α::Float64)
A = zeros(Float64, n, n+1)
function wj(n::P, r::P, α::T) where {P <: Integer, T <: Real}
A = zeros(T, n, n+1)
for iw=1:r-2
for jw in 1:iw+1
A[iw, jw]=w(iw+1-jw, iw+1, iw, n, α)
Expand All @@ -366,18 +366,18 @@ function wj(n::Int64, r::Int64, α::Float64)
return A
end

function w(i::Int64, r, j, n, a)
function w(i::P, r::P, j, n::P, a) where {P <: Integer}
ar = ones(r-1)
br = ones(r-1)
for lj in 1:r-2
jj = r-lj-1
kj = r-2
tj = i-1
aj = collect(Int64, 0:tj)
bj = collect(Int64, tj+2:kj+1)
aj = collect(Int, 0:tj)
bj = collect(Int, tj+2:kj+1)
cj = [aj; bj]
dj = collect(Int64, -1:tj-1)
ej = collect(Int64, tj+1:kj)
dj = collect(Int, -1:tj-1)
ej = collect(Int, tj+1:kj)
fj = [dj; ej]
yj = binomial.(cj, jj)
pj = binomial.(fj, jj)
Expand All @@ -404,7 +404,7 @@ function w(i::Int64, r, j, n, a)
return s
end

function fj(n::Int64, T, fy)
function fj(n::Integer, T, fy)
B = zeros(n+1)
for i2 = 1:n+1
B[i2] = fy(T*(i2-1)/n)
Expand All @@ -413,10 +413,10 @@ function fj(n::Int64, T, fy)
end

function fracdiff(f::FunctionAndNumber,
alpha::Float64,
alpha::T,
end_point::Real,
h::Float64,
::CaputoL1)
h::T,
::CaputoL1) where {T <: Real}
result = zero(Float64)
n = floor(Int, end_point/h)

Expand All @@ -432,11 +432,11 @@ end


function fracdiff(f::FunctionAndNumber,
alpha::Float64,
alpha::T,
end_point::Real,
h::Float64,
::CaputoL2)
result = zero(Float64)
h::T,
::CaputoL2) where {T <: Real}
result = zero(T)
n = round(Int, end_point/h)

for k=-1:n
Expand All @@ -445,7 +445,7 @@ function fracdiff(f::FunctionAndNumber,
return h^(-alpha)/gamma(3-alpha)*result
end

function wkcoeff(k::Int64, n::Int64, alpha::Float64)
function wkcoeff(k::P, n::P, alpha::T) where {P <: Integer, T <: Real}
if k == -1
return 1
elseif k == 0
Expand Down
72 changes: 36 additions & 36 deletions src/Derivative/RL.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,15 +111,15 @@ struct RLDiffL2C <: RLDiff end

# Riemann Liouville sense L1 method
function fracdiff(f::FunctionAndNumber,
α::Float64,
α::T,
end_point::Real,
h::Float64,
::RLDiffL1)
h::T,
::RLDiffL1) where {T <: Real}
#checks(f, α, 0, end_point)
typeof(f) <: Real ? (end_point == 0 ? (return 0) : (return f/sqrt(pi*end_point))) : nothing
end_point == 0 ? (return 0) : nothing

summation = 0
summation = zero(T)
n = floor(Int, end_point/h)

@fastmath @inbounds @simd for i ∈ 0:n-1
Expand All @@ -130,7 +130,7 @@ function fracdiff(f::FunctionAndNumber,
return result
end

function fracdiff(f::FunctionAndNumber, α::Float64, end_point::AbstractArray, h::Float64, ::RLDiffL1)::Vector
function fracdiff(f::FunctionAndNumber, α::T, end_point::AbstractArray, h::T, ::RLDiffL1)::Vector where {T <: Real}
result = map(x->fracdiff(f, α, x, h, RLDiffL1()), end_point)
return result
end
Expand All @@ -139,10 +139,10 @@ end


function fracdiff(f::FunctionAndNumber,
α::Float64,
α::T,
end_point::Real,
h::Float64,
::RLDiffMatrix)
h::T,
::RLDiffMatrix) where {T <: Real}
N = round(Int, end_point/h+1)
@views tspan = collect(0:h:end_point)
return B(N, α, h)*f.(tspan)
Expand All @@ -165,7 +165,7 @@ function B(N, p)
return result
end
# Multiple dispatch for assigning step size *h*.
function B(N, p, h::Float64)
function B(N, p, h::T) where {T <: Real}
result = B(N, p)

return h^(-p)*result
Expand All @@ -181,21 +181,21 @@ Page 57
Linear Spline Interpolation
=#
function fracdiff(f::FunctionAndNumber,
α::Float64,
α::T,
x::Real,
h::Float64,
::RLLinearSplineInterp)
h::T,
::RLLinearSplineInterp) where {T <: Real}
typeof(f) <: Real ? (x == 0 ? (return 0) : (return f/sqrt(pi*x))) : nothing
x == 0 ? (return 0) : nothing
N = round(Int, x/h)
result = 0
result = zero(T)
@fastmath @inbounds @simd for k = 0:(N+1)
result += z̄ₘₖ(N, k, α)*f(k*h)
end
return 1/(gamma(4-α)h^α)*result
end

function z̄ₘₖ(m, k, α)
function z̄ₘₖ(m::P, k::P, α::T) where {P <: Integer, T <: Real}
if k ≤ m-1
return c̄ⱼₖ(m-1, k, α)-2*c̄ⱼₖ(m, k, α)+c̄ⱼₖ(m+1, k, α)
elseif k == m
Expand All @@ -217,43 +217,43 @@ function c̄ⱼₖ(j, k, α)
end
end

function fracdiff(f::FunctionAndNumber, α::Float64, end_point::AbstractArray, h::Float64, ::RLLinearSplineInterp)::Vector
function fracdiff(f::FunctionAndNumber, α::T, end_point::AbstractArray, h::T, ::RLLinearSplineInterp)::Vector where {T <: Real}
result = map(x->fracdiff(f, α, x, h, RLLinearSplineInterp()), end_point)
return result
end

function fracdiff(f::FunctionAndNumber,
α::Float64,
α::T,
start_point::Real,
end_point::Real,
h::Float64,
::RLG1)
h::T,
::RLG1) where {T <: Real}
typeof(f) <: Number ? (end_point == 0 ? (return 0) : (return f/sqrt(pi*end_point))) : nothing
end_point == 0 ? (return 0) : nothing

N = round(Int, (end_point-start_point)/h)

result = zero(Float64)
result = zero(T)
@fastmath @inbounds @simd for j = 0:N-1
result += gamma(j-α)/gamma(j+1)*f(end_point-j*h)
end

return h^(-α)/gamma(-α)*result
end

fracdiff(f::FunctionAndNumber, α, end_point, h::Float64, ::RLG1) = fracdiff(f::FunctionAndNumber, α, 0, end_point, h::Float64, RLG1())
fracdiff(f::FunctionAndNumber, α, end_point, h::T, ::RLG1) where {T <: Real} = fracdiff(f::FunctionAndNumber, α, 0, end_point, h::Float64, RLG1())

function fracdiff(f::FunctionAndNumber,
α::Float64,
α::T,
point::Real,
h::Float64,
::RLD)
h::T,
::RLD) where {T <: Real}
typeof(f) <: Real ? (point == 0 ? (return 0) : (return f/sqrt(pi*point))) : nothing
point == 0 ? (return 0) : nothing

N = round(Int, point/h)

result = zero(Float64)
result = zero(T)
@fastmath @inbounds @simd for k = 0:N
result += ωₖₙ(N, k, α)*f(point-k*h)
end
Expand All @@ -262,7 +262,7 @@ function fracdiff(f::FunctionAndNumber,

end

function ωₖₙ(n, k, α)
function ωₖₙ(n::P, k::P, α) where {P <: Integer}
temp = 0
if k == 0
temp = -1
Expand All @@ -275,18 +275,18 @@ function ωₖₙ(n, k, α)
end

function fracdiff(f::FunctionAndNumber,
α::Float64,
α::T,
point::Real,
h::Float64,
::RLDiffL2C)
h::T,
::RLDiffL2C) where {T <: Real}
N = round(Int, point/h)
temp = zero(Float64)
temp = zero(T)
for k=-1:N+1
temp += Ŵ(k, α, N)*f((N-k)*h)
end
end

function Ŵ(k::Int64, α::Float64, N::Int64)
function Ŵ(k::P, α::T, N::P) where {P <: Integer, T <: Real}
expo = 2-α
if k == -1
return 1
Expand All @@ -312,21 +312,21 @@ end


function fracdiff(f::FunctionAndNumber,
alpha::Float64,
alpha::T,
point::Real,
h::Float64,
::RLDiffL2)
h::T,
::RLDiffL2) where {T <: Real}

n = round(Int64, point/h)
result = zero(Float64)
n = round(Int, point/h)
result = zero(T)
for k=-1:n
result += WK(k, alpha, n)*f((n-k)*h)
end

return f(0)*point^(-alpha)/gamma(1-alpha) + derivative(f, 0)*point^(1-alpha)/gamma(2-alpha) + h^(-alpha)/gamma(3-alpha)*result
end

function WK(k::Int64, alpha::Float64, n::Int64)
function WK(k::P, alpha::T, n::P) where {P <: Integer, T <: Real}
if k == -1
return 1
elseif k == 0
Expand Down

0 comments on commit e4f1c5f

Please sign in to comment.