Refactor FDDEProblem and solvers

Project.toml

@@ -16,6 +16,7 @@ LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"

docs/make.jl

@@ -12,6 +12,8 @@ makedocs(;
         canonical="https://SciFracX.github.io/FractionalDiffEq.jl",
         assets=["assets/favicon.ico"],
     ),
+    warnonly = [:missing_docs, :cross_references],
+    doctest = false,
     pages = [
         "FractionalDiffEq.jl" => "index.md",
         "Get Started" => "get_started.md",

docs/src/algorithms.md

@@ -8,56 +8,46 @@ Pages = ["algorithms.md"]
 
 ### Single Term FODE
 
-```@docs
-FractionalDiffEq.PECE
-FractionalDiffEq.Euler
-FractionalDiffEq.PIEX
-FractionalDiffEq.AtanganaSeda
-```
+PECE
+Euler
+PIEX
+AtanganaSeda
+
 
 ### Multi-Term FODE
 
-```@docs
-FractionalDiffEq.FODEMatrixDiscrete
-FractionalDiffEq.ClosedForm
-FractionalDiffEq.ClosedFormHankelM
-FractionalDiffEq.PIPECE
-FractionalDiffEq.PIRect
-FractionalDiffEq.PITrap
-```
+FODEMatrixDiscrete
+ClosedForm
+ClosedFormHankelM
+PIPECE
+PIRect
+PITrap
 
 ### System of FODE
 
-```@docs
-FractionalDiffEq.NonLinearAlg
-FractionalDiffEq.PECE
-FractionalDiffEq.GL
-FractionalDiffEq.FLMMBDF
-FractionalDiffEq.FLMMNewtonGregory
-FractionalDiffEq.FLMMTrap
-FractionalDiffEq.PIEX
-FractionalDiffEq.NewtonPolynomial
-FractionalDiffEq.AtanganaSedaAB
-```
+NonLinearAlg
+PECE
+GL
+FLMMBDF
+FLMMNewtonGregory
+FLMMTrap
+PIEX
+NewtonPolynomial
+AtanganaSedaAB
 
-## Fractional Delay Differential Equatinos
+## Fractional Delay Differential Equations
+
+DelayPECE
+DelayPI
+MatrixForm
+DelayABM
 
-```@docs
-FractionalDiffEq.DelayPECE
-FractionalDiffEq.DelayPI
-FractionalDiffEq.MatrixForm
-FractionalDiffEq.DelayABM
-```
 
 ## Distributed Order Differential Equations
 
-```@docs
-FractionalDiffEq.DOMatrixDiscrete
-```
+DOMatrixDiscrete
 
 ## Fractional Differences Equations
 
-```@docs
-FractionalDiffEq.PECEDifference
-FractionalDiffEq.GL
-```
+PECEDifference
+GL

docs/src/problems.md

@@ -4,15 +4,11 @@
 Pages = ["problems.md"]
 ```
 
-The general fractional differential equations problem type:
-
-```@docs
-FractionalDiffEq.FDEProblem
-```
+The general fractional differential equations problem type.
 
 ## FODE Problem
 
-Fractional ordinary problems type, we can classified them as Single-Term and Multi-Term problems:
+Fractional ordinary problems type, we can classify them as Single-Term and Multi-Term problems:
 
 ```math
 D^{\alpha}y(t)=f(t, y)
@@ -22,49 +18,31 @@ D^{\alpha}y(t)=f(t, y)
 \sum_{s=0}^pc_sD^{\beta_s}y=f
 ```
 
-```@docs
-FractionalDiffEq.SingleTermFODEProblem
-FractionalDiffEq.MultiTermsFODEProblem
-FractionalDiffEq.FODESystem
-```
+SingleTermFODEProblem
+MultiTermsFODEProblem
+FODESystem
 
 ### FFODE Problem
 
-```@docs
-FractionalDiffEq.FFPODEProblem
-FractionalDiffEq.FFEODEProblem
-FractionalDiffEq.FFMODEProblem
-```
-
-## FPDE Problem
-
-```@docs
-FractionalDiffEq.FPDEProblem
-```
+FFPODEProblem
+FFEODEProblem
+FFMODEProblem
 
 ## FDDE Problem
 
-```@docs
-FractionalDiffEq.FDDEProblem
-FractionalDiffEq.FDDEMatrixProblem
-FractionalDiffEq.FDDESystem
-```
+FDDEProblem
+FDDEMatrixProblem
+FDDESystem
 
 ## DODE Problem
 
-```@docs
 FractionalDiffEq.DODEProblem
-```
 
 ## Fractional Difference Problem
 
-```@docs
-FractionalDiffEq.FractionalDifferenceProblem
-FractionalDiffEq.FractionalDifferenceSystem
-```
+FractionalDifferenceProblem
+FractionalDifferenceSystem
 
 ## FIE Problem
 
-```@docs
-FractionalDiffEq.FIEProblem
-```
+FIEProblem

src/FractionalDiffEq.jl

@@ -20,6 +20,8 @@ include("types/problems.jl")
 include("types/algorithms.jl")
 include("types/solutions.jl")
 
+include("types/problem_utils.jl")
+
 # Multi-terms fractional ordinary differential equations
 include("multitermsfode/matrix.jl")
 include("multitermsfode/hankelmatrix.jl")

src/delay/DelayABM.jl

@@ -1,5 +1,5 @@
 """
-    solve(FDDE::FDDEProblem, h, DelayABM())
+    solve(FDDE::FDDEProblem, dt, DelayABM())
 
 Use the Adams-Bashforth-Moulton method to solve fractional delayed differential equations.
 
@@ -16,44 +16,45 @@ Use the Adams-Bashforth-Moulton method to solve fractional delayed differential
 struct DelayABM <: FDDEAlgorithm end
 #FIXME: There are still some improvments about initial condition
 #FIXME: Fix DelayABM method for FDDESystem : https://www.researchgate.net/publication/245538900_A_Predictor-Corrector_Scheme_For_Solving_Nonlinear_Delay_Differential_Equations_Of_Fractional_Order
-#FIXME: Also the problem definition f(t, ϕ, y) or f(t, y, ϕ)?
-function solve(FDDE::FDDEProblem, h, ::DelayABM)
-    @unpack f, ϕ, α, τ, tspan = FDDE
-    N::Int = round(Int, tspan/h)
-    Ndelay::Int = round(Int, τ/h)
+#FIXME: Also the problem definition f(t, h, y) or f(t, y, h)?
+function solve(FDDE::FDDEProblem, dt, ::DelayABM)
+    @unpack f, h, order, constant_lags, tspan, p = FDDE
+    τ = constant_lags[1]
+    N::Int = round(Int, tspan[2]/dt)
+    Ndelay::Int = round(Int, τ/dt)
     x1 = zeros(Float64, Ndelay+N+1)
     x = zeros(Float64, Ndelay+N+1)
     #x1[Ndelay+N+1] = 0
 
     #x[Ndelay+N+1] = 0
 
     # History function handling
-    #FIXME: When the value of history function ϕ is different with the initial value?
-    if typeof(ϕ) <: Number
-        x[1:Ndelay] = ϕ*ones(Ndelay)
-    elseif typeof(ϕ) <: Function
-        x[Ndelay] = ϕ(0)
-        x[1:Ndelay-1] .= ϕ(collect(Float64, -h*(Ndelay-1):h:(-h)))
+    #FIXME: When the value of history function h is different with the initial value?
+    if typeof(h) <: Number
+        x[1:Ndelay] = h*ones(Ndelay)
+    elseif typeof(h) <: Function
+        x[Ndelay] = h(p, 0)
+        x[1:Ndelay-1] .= h(p, collect(Float64, -dt*(Ndelay-1):dt:(-dt)))
     end
 
     x0 = copy(x[Ndelay])
 
 
-    x1[Ndelay+1] = x0 + h^α*f(0, x[1], x0)/(gamma(α)*α)
+    x1[Ndelay+1] = x0 + dt^order*f(x[1], x0, p, 0)/(gamma(order)*order)
 
-    x[Ndelay+1] = x0 + h^α*(f(0, x[1], x[Ndelay+1]) + α*f(0, x[1], x0))/gamma(α+2)
+    x[Ndelay+1] = x0 + dt^order*(f(x[1], x[Ndelay+1], p, 0) + order*f(x[1], x0, p, 0))/gamma(order+2)
 
     @fastmath @inbounds @simd for n=1:N 
-        M1=(n^(α+1)-(n-α)*(n+1)^α)*f(0, x[1], x0)
+        M1=(n^(order+1)-(n-order)*(n+1)^order)*f(x[1], x0, p, 0)
 
-        N1=((n+1)^α-n^α)*f(0, x[1], x0)
+        N1=((n+1)^order-n^order)*f(x[1], x0, p, 0)
 
         @fastmath @inbounds @simd for j=1:n   
-            M1 = M1+((n-j+2)^(α+1)+(n-j)^(α+1)-2*(n-j+1)^(α+1))*f(0, x[j], x[Ndelay+j])
-            N1 = N1+((n-j+1)^α-(n-j)^α)*f(0, x[j], x[Ndelay+j])
+            M1 = M1+((n-j+2)^(order+1)+(n-j)^(order+1)-2*(n-j+1)^(order+1))*f(x[j], x[Ndelay+j], p, 0)
+            N1 = N1+((n-j+1)^order-(n-j)^order)*f(x[j], x[Ndelay+j], p, 0)
         end
-        x1[Ndelay+n+1] = x0+h^α*N1/(gamma(α)*α)
-        x[Ndelay+n+1] = x0+h^α*(f(0, x[n+1], x[Ndelay+n+1])+M1)/gamma(α+2)
+        x1[Ndelay+n+1] = x0+dt^order*N1/(gamma(order)*order)
+        x[Ndelay+n+1] = x0+dt^order*(f(x[n+1], x[Ndelay+n+1], p, 0)+M1)/gamma(order+2)
     end
 
     xresult = zeros(N-Ndelay+1)
@@ -83,12 +84,12 @@ DelayABM method for system of fractional delay differential equations.
 ```
 =#
 
-function solve(FDDESys::FDDESystem, h, ::DelayABM)
+function solve(FDDESys::FDDESystem, dt, ::DelayABM)
     @unpack f, ϕ, α, τ, T = FDDESys
     len = length(ϕ)
-    N::Int = round(Int, T/h)
-    Ndelay = round(Int, τ/h)
-    t = collect(Float64, 0:h:(T-τ))
+    N::Int = round(Int, T/dt)
+    Ndelay = round(Int, τ/dt)
+    t = collect(Float64, 0:dt:(T-τ))
     x1 = zeros(Ndelay+N+1, len)
     x = zeros(Ndelay+N+1, len)
     du = zeros(len)
@@ -103,15 +104,15 @@ function solve(FDDESys::FDDESystem, h, ::DelayABM)
     for i=1:len
         αi = α[i]
         f(du, x0, x[1, :], 0)
-        x1[Ndelay+1, i] = x0[i] + h^αi*du[i]/(gamma(αi)*αi)
+        x1[Ndelay+1, i] = x0[i] + dt^αi*du[i]/(gamma(αi)*αi)
     end
 
     for i=1:len
         αi = α[i]
         f(du, x[Ndelay+1, :], x[1, :], 0)
         du1 = copy(du)
         f(du, x0, x[1, :], 0)
-        x[Ndelay+1, i] = x0[i] + h^αi*(du1[i] + αi*du[i])/gamma(αi+2)
+        x[Ndelay+1, i] = x0[i] + dt^αi*(du1[i] + αi*du[i])/gamma(αi+2)
     end
 
     M1=zeros(len); N1=zeros(len)
@@ -134,9 +135,9 @@ function solve(FDDESys::FDDESystem, h, ::DelayABM)
 
         for i=1:len
             αi = α[i]
-            x1[Ndelay+n+1, i] = x0[i]+h^αi*N1[i]/(gamma(αi)*αi)
+            x1[Ndelay+n+1, i] = x0[i]+dt^αi*N1[i]/(gamma(αi)*αi)
             f(du, x[Ndelay+n+1, :], x[n+1, :], 0)
-            x[Ndelay+n+1, i] = x0[i]+h^αi*(du[i]+M1[i])/gamma(αi+2)
+            x[Ndelay+n+1, i] = x0[i]+dt^αi*(du[i]+M1[i])/gamma(αi+2)
         end
     end