Merge pull request #81 from ErikQQY/refactor

Add DiffEqBase and SciMLBase

Project.toml

@@ -4,6 +4,7 @@ authors = ["Qingyu Qu <[email protected]>"]
 version = "0.3.1"
 
 [deps]
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
@@ -12,10 +13,12 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 SpecialMatrices = "928aab9d-ef52-54ac-8ca1-acd7ca42c160"
 ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
+TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]

src/FractionalDiffEq.jl

@@ -1,6 +1,7 @@
 module FractionalDiffEq
 
 using LinearAlgebra
+using SciMLBase, DiffEqBase
 using SpecialFunctions
 using SparseArrays
 using InvertedIndices: Not
@@ -13,19 +14,12 @@ using ToeplitzMatrices
 using RecipesBase
 using ForwardDiff
 using Polynomials: Polynomial
+using TruncatedStacktraces
 
 include("types/problems.jl")
 include("types/algorithms.jl")
 include("types/solutions.jl")
 
-
-# Single-term fractional ordinary differential equations
-include("singletermfode/PECE.jl")
-include("singletermfode/product_integral.jl")
-include("singletermfode/GL.jl")
-include("singletermfode/atangana_seda.jl")
-include("singletermfode/Euler.jl")
-
 # Multi-terms fractional ordinary differential equations
 include("multitermsfode/matrix.jl")
 include("multitermsfode/hankelmatrix.jl")
@@ -80,7 +74,7 @@ export solve, FDEProblem
 export FDDEProblem, FDDESystem, FDDEMatrixProblem
 
 # FODE problems
-export SingleTermFODEProblem, MultiTermsFODEProblem, FODESystem, DODEProblem, FFPODEProblem, FFEODEProblem, FFMODEProblem
+export FODEProblem, MultiTermsFODEProblem, DODEProblem, FFPODEProblem, FFEODEProblem, FFMODEProblem
 
 # Fractional Discrete probelms
 export FractionalDiscreteProblem, FractionalDiscreteSystem
@@ -94,14 +88,15 @@ export FODESolution, FDifferenceSolution, DODESolution, FFMODESolution
 export FODESystemSolution, FDDESystemSolution, FFMODESystem
 
 # FODE solvers
-export PIEX, PIPECE, PIRect, PITrap
+export PIPECE, PIRect, PITrap
 export PECE, FODEMatrixDiscrete, ClosedForm, ClosedFormHankelM, GL
-export AtanganaSeda, AtanganaSedaAB
-export Euler
+export AtanganaSedaAB
+#export Euler
 
 # System of FODE solvers
 export NonLinearAlg, FLMMBDF, FLMMNewtonGregory, FLMMTrap, PIEX, NewtonPolynomial
 export AtanganaSedaCF
+export AtanganaSeda
 
 # FDDE solvers
 export DelayPECE, DelayABM, MatrixForm
@@ -129,6 +124,4 @@ export FOLyapunov, FOLE
 # Distributed order auxiliary SpecialFunctions
 export DOB, DOF, DORANORT, isFunction
 
-export ourfft, ourifft, rowfft, FastConv
-
 end

src/ffode/atangana_seda.jl

@@ -1,3 +1,5 @@
+struct AtanganaSeda <: FODESystemAlgorithm end
+
 function solve(prob::Union{FFMODEProblem, FFMODESystem}, h::Float64, ::AtanganaSeda)
     if typeof(prob.order[2]) <: Function
         solve_cf_variable_ffodesystem(prob, h)

src/fodesystem/Euler.jl

@@ -1,7 +1,9 @@
-function solve(prob::FODESystem, h, ::Euler)
-    @unpack f, α, u0, tspan, p = prob
+struct Euler <: FODESystemAlgorithm end
+
+function solve(prob::FODEProblem, h, ::Euler)
+    @unpack f, order, u0, tspan, p = prob
     t0=tspan[1]; tfinal=tspan[2]
-    α = α[1]
+    order = order[1]
     t = collect(Float64, t0:h:tfinal)
     N::Int = ceil(Int, (tfinal-t0)/h)
     l = length(u0)
@@ -12,9 +14,9 @@ function solve(prob::FODESystem, h, ::Euler)
         tmp = zeros(Float64, l)
         for j=1:n
             f(tmp, result[:, j], p, t[j])
-            temp = ((n-j+1)^α-(n-j)^α)*tmp
+            temp = ((n-j+1)^order-(n-j)^order)*tmp
         end
-        result[:, n+1] = result[:, 1] + h^α/gamma(α+1)*temp
+        result[:, n+1] = result[:, 1] + h^order/gamma(order+1)*temp
     end
     return result
 end

src/fodesystem/GLWithMemory.jl

@@ -1,7 +1,7 @@
 """
 # Usage
 
-    solve(prob::FODESystem, h, GL())
+    solve(prob::FODEProblem, h, GL())
 
 Use Grunwald Letnikov difference method to solve system of system of FODE.
 
@@ -18,29 +18,30 @@ doi={10.1109/MOCAST.2019.8742063}}
 
 Python version by https://github.com/DClementeL/Grunwald_Letnikov
 """
-# Grunwald Letnikov discretization method dispatch for FODESystem
+# Grunwald Letnikov discretization method dispatch for FODEProblem
 # struct GLWithMemory <: FractionalDiffEqAlgorithm end
+struct GL <: FODESystemAlgorithm end
 
-function solve(prob::FODESystem, h, ::GL)
+function solve(prob::FODEProblem, h, ::GL)
     # GL method is only for same order FODE
-    @unpack f, α, u0, tspan, p = prob
+    @unpack f, order, u0, tspan, p = prob
     t0 = tspan[1]; T = tspan[2]
     t = collect(Float64, t0:h:T)
-    α = α[1]
-    hα = h^α[1]
+    order = order[1]
+    horder = h^order[1]
     n::Int64 = floor(Int64, (T-t0)/h)+1
     l = length(u0)
 
     # Initialize solution
     result = zeros(Float64, length(u0), n)
     result[:, 1] = u0
 
-    # generating generalized binomial Cα
-    Cα = zeros(Float64, n)
-    Cα[1] = 1
+    # generating generalized binomial Corder
+    Corder = zeros(Float64, n)
+    Corder[1] = 1
 
     @fastmath @inbounds @simd for j in range(2, n, step=1)
-        Cα[j] = (1-(1+α)/(j-1))*Cα[j-1]
+        Corder[j] = (1-(1+order)/(j-1))*Corder[j-1]
     end
 
     du = zeros(Float64, l)
@@ -50,12 +51,12 @@ function solve(prob::FODESystem, h, ::GL)
 
         @fastmath @inbounds @simd for j in range(1, k-1, step=1)
             for i in eachindex(summation)
-                summation[i] += Cα[j+1]*result[i, k-j]
+                summation[i] += Corder[j+1]*result[i, k-j]
             end
         end
 
         f(du, result[:, k-1], p, t0+(k-1)*h)
-        result[:, k] = @. hα*du-summation
+        result[:, k] = @. horder*du-summation
     end
     return FODESystemSolution(t, result)
 end

src/fodesystem/NonLinear.jl

@@ -1,7 +1,7 @@
 """
 # Usage
 
-    solve(prob::FODESystem, h, NonLinearAlg())
+    solve(prob::FODEProblem, h, NonLinearAlg())
 
 Nonlinear algorithm for nonlinear fractional differential equations.
 
@@ -11,16 +11,16 @@ Dingyu Xue, Northeastern University, China ISBN:9787030543981
 """
 struct NonLinearAlg <: FODESystemAlgorithm end
 
-function solve(prob::FODESystem, h, ::NonLinearAlg, L0=1e10)
-    @unpack f, α, u0, tspan, p = prob
+function solve(prob::FODEProblem, h, ::NonLinearAlg, L0=1e10)
+    @unpack f, order, u0, tspan, p = prob
     t0 = tspan[1]; T = tspan[2]
     time = collect(t0:h:T)
     n = length(u0)
     m = round(Int, (T-t0)/h)+1
     g = genfun(1)
     g = g[:]
     u0 = u0[:]
-    ha = h.^α
+    ha = h.^order
     z = zeros(Float64, n, m)
     x1::Vector{Float64} = copy(u0)
 
@@ -29,7 +29,7 @@ function solve(prob::FODESystem, h, ::NonLinearAlg, L0=1e10)
     W = zeros(n, SetMemoryEffect) #Initializing W a n*m matrix
 
     @fastmath @inbounds @simd for i = 1:n
-        W[i, :] = getvec(α[i], SetMemoryEffect, g)
+        W[i, :] = getvec(order[i], SetMemoryEffect, g)
     end
 
     du = zeros(n)
@@ -60,12 +60,12 @@ function genfun(p)
     return (1 .-a')*inv(A')
 end
 
-function getvec(α, n, g)
+function getvec(order, n, g)
     p = length(g)-1
-    b = 1 + α
+    b = 1 + order
     g0 = g[1]
     w = Float64[]
-    push!(w, g[1]^α)
+    push!(w, g[1]^order)
 
     @fastmath @inbounds @simd for m = 2:p
         M = m-1

src/fodesystem/PIPECE.jl

@@ -1,6 +1,6 @@
 #=
 """
-    solve(prob::FODESystem, h, PECE())
+    solve(prob::FODEProblem, h, PECE())
 
 Use the Adams-Bashforth-Moulton method to solve the system of FODEs.
 
@@ -29,27 +29,30 @@ mutable struct M
     an_fft
     bn_fft
 end
+
+struct PECE <: FODESystemAlgorithm end
+
 #TODO: Rename as PIPECE
-function solve(prob::FODESystem, h, ::PECE)
-    @unpack f, α, u0, tspan, p = prob
+function solve(prob::FODEProblem, h, ::PECE)
+    @unpack f, order, u0, tspan, p = prob
     t0 = tspan[1]; T = tspan[2]
     METH = M(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)# Initialization
     mu_tol = 1.0e-6
     mu = 1
-    α = α[:]
+    order = order[:]
 
     # issue [#64](https://github.com/SciFracX/FractionalDiffEq.jl/issues/64)
-    max_order = findmax(α)[1]
+    max_order = findmax(order)[1]
     if max_order > 1
         @error "This method doesn't support high order FDEs"
     end
 
 
     # Check compatibility size of the problem with number of fractional orders
-    alpha_length = length(α)
+    alpha_length = length(order)
     problem_size = size(u0, 1)
 
-    m_alpha = ceil.(Int, α)
+    m_alpha = ceil.(Int, order)
     m_alpha_factorial = zeros(alpha_length, maximum(m_alpha))
     for i = 1 : alpha_length
         for j = 0 : m_alpha[i]-1
@@ -78,7 +81,7 @@ function solve(prob::FODESystem, h, ::PECE)
     bn = zeros(alpha_length, NNr+1); an = copy(bn); a0 = copy(bn)
     for i_alpha = 1:alpha_length
         find_alpha = Float64[]
-        if α[i_alpha] == α[1:i_alpha-1]
+        if order[i_alpha] == order[1:i_alpha-1]
             push!(find_alpha, i_alpha)
         end
 
@@ -87,16 +90,16 @@ function solve(prob::FODESystem, h, ::PECE)
             an[i_alpha, :] = an[find_alpha[1], :]
             a0[i_alpha, :] = a0[find_alpha[1], :]
         else
-            nalpha = nvett.^α[i_alpha]
+            nalpha = nvett.^order[i_alpha]
             nalpha1 = nalpha.*nvett
             bn[i_alpha, :] = nalpha[2:end] - nalpha[1:end-1]
             an[i_alpha, :] = [1; (nalpha1[1:end-2] - 2*nalpha1[2:end-1] + nalpha1[3:end]) ]
-            a0[i_alpha, :] = [0; nalpha1[1:end-2]-nalpha[2:end-1].*(nvett[2:end-1].-α[i_alpha].-1)]
+            a0[i_alpha, :] = [0; nalpha1[1:end-2]-nalpha[2:end-1].*(nvett[2:end-1].-order[i_alpha].-1)]
         end
     end
     METH.bn = bn; METH.an = an; METH.a0 = a0
-    METH.halpha1 = h.^α./gamma.(α.+1)
-    METH.halpha2 = h.^α./gamma.(α.+2)
+    METH.halpha1 = h.^order./gamma.(order.+1)
+    METH.halpha2 = h.^order./gamma.(order.+2)
     METH.mu = mu; METH.mu_tol = mu_tol
 
     # Evaluation of FFT of coefficients of the PECE method
@@ -116,7 +119,7 @@ function solve(prob::FODESystem, h, ::PECE)
             coef_end = 2^l*r
             for i_alpha = 1 : alpha_length
                 find_alpha = Float64[]
-                if α[i_alpha] == α[1:i_alpha-1]
+                if order[i_alpha] == order[1:i_alpha-1]
                     push!(find_alpha, i_alpha)
                 end
                 if isempty(find_alpha) == false

src/fodesystem/atangana_seda.jl

@@ -5,13 +5,13 @@ Solve Atangana-Baleanu fractional order differential equations using Newton Poly
 """
 struct AtanganaSedaAB <: FODESystemAlgorithm end
 
-function solve(prob::FODESystem, h, ::AtanganaSedaAB)
-    @unpack f, α, u0, tspan, p = prob
-    α = α[1]
+function solve(prob::FODEProblem, h, ::AtanganaSedaAB)
+    @unpack f, order, u0, tspan, p = prob
+    order = order[1]
     t0 = tspan[1]; tfinal = tspan[2]
     t = collect(t0:h:tfinal)
     N::Int = ceil(Int, (tfinal-t0)/h)
-    AB = 1-α+α/gamma(α)
+    AB = 1-order+order/gamma(order)
     t = collect(Float64, t0:h:tfinal)
 
     l = length(u0)
@@ -36,16 +36,16 @@ function solve(prob::FODESystem, h, ::AtanganaSedaAB)
         f(temp1, result[:, n-1], p, t[n-1])
         for j=3:n
             f(temp1, result[:, j-2], p, t[j-2])
-            temptemp1 += ((n+1-j)^α-(n-j)^α)*temp1
+            temptemp1 += ((n+1-j)^order-(n-j)^order)*temp1
             f(temp2, result[:, j-1], p, t[j-1])
-            temptemp2 += (temp2-temp1)*((n+1-j)^α*(n-j+3+2*α)-(n-j)^α*(n-j+3+3*α))
+            temptemp2 += (temp2-temp1)*((n+1-j)^order*(n-j+3+2*order)-(n-j)^order*(n-j+3+3*order))
             f(temp3, result[:, j], p, t[j])
-            temptemp3 += (temp3-2*temp2+temp1)*((n-j+1)^α*(2*(n-j)^2+(3*α+10)*(n-j)+2*(α)^2+9*α+12)-(n-j)^α*(2*(n-j)^2+(5*α+10)*(n-j)+6*α^2+18*α+12))
+            temptemp3 += (temp3-2*temp2+temp1)*((n-j+1)^order*(2*(n-j)^2+(3*order+10)*(n-j)+2*(order)^2+9*order+12)-(n-j)^order*(2*(n-j)^2+(5*order+10)*(n-j)+6*order^2+18*order+12))
         end
         temptemptemp = zeros(l)
         f(temptemptemp, result[:, n], p, t[n])
 
-        result[:, n+1] = u0+(1-α)/AB*temptemptemp+((h^α)*α/(AB*gamma(α+1)))*temptemp1 + ((h^α)*α/(AB*gamma(α+2)))*temptemp2+((h^α)*α/(2*AB*gamma(α+3)))*temptemp3
+        result[:, n+1] = u0+(1-order)/AB*temptemptemp+((h^order)*order/(AB*gamma(order+1)))*temptemp1 + ((h^order)*order/(AB*gamma(order+2)))*temptemp2+((h^order)*order/(2*AB*gamma(order+3)))*temptemp3
     end
     return FODESystemSolution(t, result)
 end
@@ -65,12 +65,12 @@ https://doi.org/10.1016/c2020-0-02711-8
 """
 struct AtanganaSedaCF <: FODESystemAlgorithm end
 #FIXME: Tests
-function solve(prob::FODESystem, h, ::AtanganaSedaCF)
-    @unpack f, α, u0, tspan, p = prob
+function solve(prob::FODEProblem, h, ::AtanganaSedaCF)
+    @unpack f, order, u0, tspan, p = prob
     t0 = tspan[1]; tfinal = tspan[2]
     t=collect(Float64, t0:h:tfinal)
-    α=α[1]
-    M=1-α+α/gamma(α)
+    order=order[1]
+    M=1-order+order/gamma(order)
     N=ceil(Int, (tfinal-t0)/h)
     l=length(u0)
     result = zeros(l, N+1)
@@ -85,7 +85,7 @@ function solve(prob::FODESystem, h, ::AtanganaSedaCF)
         f(tempn1, result[:, n-1], p, t[n-1])
 
         f(tempn2, result[:, n-2], p, t[n-2])
-        result[:, n+1] = result[:, n] + (1-α)/M*(tempn-tempn1)+α*M*h*(23/12*tempn-4/3*tempn1+5/12*tempn2)
+        result[:, n+1] = result[:, n] + (1-order)/M*(tempn-tempn1)+order*M*h*(23/12*tempn-4/3*tempn1+5/12*tempn2)
     end
     return FODESystemSolution(t, result)
 end