Refactor FDDE Solvers

README.md

@@ -62,9 +62,9 @@ So we can use FractionalDiffEq.jl to solve the problem:
 ```julia
 using FractionalDiffEq, Plots
 fun(u, p, t) = 1-u
-u0 = [0, 0]; tspan = (0, 20); h = 0.001;
+u0 = [0, 0]; tspan = (0, 20)
 prob = SingleTermFODEProblem(fun, 1.8, u0, tspan)
-sol = solve(prob, h, PECE())
+sol = solve(prob, PECE(), dt = 0.001)
 plot(sol)
 ```
 
@@ -115,10 +115,10 @@ function chua!(du, x, p, t)
 end
 α = [0.93, 0.99, 0.92];
 x0 = [0.2; -0.1; 0.1];
-h = 0.01; tspan = (0, 100);
+tspan = (0, 100);
 p = [10.725, 10.593, 0.268, -1.1726, -0.7872]
 prob = FODESystem(chua!, α, x0, tspan, p)
-sol = solve(prob, h, NonLinearAlg())
+sol = solve(prob, NonLinearAlg(), dt = 0.01)
 plot(sol, vars=(1, 2), title="Chua System", legend=:bottomright)
 ```
 

src/delay/DelayABM.jl

@@ -1,19 +1,3 @@
-"""
-    solve(FDDE::FDDEProblem, dt, DelayABM())
-
-Use the Adams-Bashforth-Moulton method to solve fractional delayed differential equations.
-
-### References
-
-```tex
-@inproceedings{Bhalekar2011APS,
-  title={A PREDICTOR-CORRECTOR SCHEME FOR SOLVING NONLINEAR DELAY DIFFERENTIAL EQUATIONS OF FRACTIONAL ORDER},
-  author={Sachin Bhalekar and Varsha Daftardar-Gejji},
-  year={2011}
-}
-```
-"""
-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, h, y) or f(t, y, h)?

src/delay/DelayPECE.jl

@@ -34,7 +34,7 @@ Capable of solving both single term FDDE and multiple FDDE, support time varying
 """
 struct DelayPECE <: FDDEAlgorithm end
 #FIXME: What if we have an FDDE with both variable order and time varying lag??
-function solve(FDDE::FDDEProblem, step_size, ::DelayPECE)
+function solve(FDDE::FDDEProblem, step_size, alg::DelayPECE)
     # If the delays are time varying, we need to specify single delay and multiple delay
     if  FDDE.constant_lags[1] isa Function
         # Here is the PECE solver for single time varying lag
@@ -63,57 +63,59 @@ function solve(FDDE::FDDEProblem, step_size, ::DelayPECE)
     end
 end
 
-function solve_fdde_with_single_lag(FDDE::FDDEProblem, step_size)
-    @unpack f, h, order, constant_lags, p, tspan = FDDE
+function solve_fdde_with_single_lag(prob::FDDEProblem, step_size)
+    @unpack f, order, u0, h, constant_lags, p, tspan = prob
     τ = constant_lags[1]
-    iip = SciMLBase.isinplace(FDDE)
-    T = tspan[2]
-    t = collect(0:step_size:T)
-    maxn = size(t, 1)
-    yp = collect(Float64, 0:step_size:T+step_size)
-    y = copy(t)
-    y[1] = h(p, 0)
-
-    for n in 1:maxn-1
-        yp[n+1] = 0
+    iip = SciMLBase.isinplace(prob)
+    t = collect(tspan[1]:step_size:tspan[2])
+    N = round(Int, (tspan[2]-tspan[1])/step_size)
+    M = length(u0)
+    y = [zeros(eltype(u0), M) for i in 1:(N+1)]
+    yp = [zeros(eltype(u0), M) for i in 1:(N+1)]
+    fill!(y[1], h(p, 0))
+
+    for n in 1:N-1
         for j = 1:n
             if iip
-                tmp = zeros(length(yp[1]))
+                tmp = zeros(eltype(u0), M)
                 f(tmp, y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j])
-                yp[n+1] = yp[n+1]+generalized_binomials(j-1, n-1, order, step_size)*tmp
+                yp[n+1] += @. generalized_binomials(j-1, n-1, order, step_size)*tmp
             else
-                yp[n+1] = yp[n+1]+generalized_binomials(j-1, n-1, order, step_size)*f(y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j])
+                yp[n+1] += generalized_binomials(j-1, n-1, order, step_size)*f(y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j])
             end
         end
-        yp[n+1] = yp[n+1]/gamma(order)+h(p, 0)
+        yp[n+1] = yp[n+1]/gamma(order) .+ h(p, 0)
 
-        y[n+1] = 0
+        fill!(y[n+1], zero(eltype(u0)))
 
         @fastmath @inbounds @simd for j=1:n
             if iip
-                tmp = zeros(length(y[1]))
+                tmp = zeros(eltype(u0), M)
                 f(tmp, y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j])
-                y[n+1] = y[n+1]+a(j-1, n-1, order, step_size)*tmp
+                y[n+1] += a(j-1, n-1, order, step_size)*tmp
             else
-                y[n+1] = y[n+1]+a(j-1, n-1, order, step_size)*f(y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j])
+                y[n+1] += y[n+1] + a(j-1, n-1, order, step_size)*f(y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j])
             end
         end
 
         if iip
-            tmp = zeros(y[1])
+            tmp = zeros(eltype(u0), M)
             f(tmp, yp[n+1], v(h, n+1, τ, step_size, y, yp, t, p), p, t[n+1])
-            y[n+1] = y[n+1]/gamma(order)+step_size^order*tmp/gamma(order+2)+h(p, 0)
+            y[n+1] = @. y[n+1]/gamma(order)+step_size^order*tmp/gamma(order+2)+h(p, 0)
         else
-            y[n+1] = y[n+1]/gamma(order)+step_size^order*f(yp[n+1], v(h, n+1, τ, step_size, y, yp, t, p), p, t[n+1])/gamma(order+2)+h(p, 0)
+            y[n+1] = y[n+1] ./gamma(order) .+ step_size^order*f(yp[n+1], v(h, n+1, τ, step_size, y, yp, t, p), p, t[n+1])/gamma(order+2) .+ h(p, 0)
         end
     end
 
-    V = copy(t)
-    @fastmath @inbounds @simd for n = 1:maxn-1
+    V = copy(y)
+    @fastmath @inbounds @simd for n = 1:N-1
         # The delay term
         V[n] = v(h, n, τ, step_size, y, yp, t, p)
     end
+
     return V, y
+
+    #return DiffEqBase.build_solution(prob, alg, t, y)
 end
 
 function a(j::Int, n::Int, order, step_size)
@@ -160,16 +162,17 @@ function v(h, n, τ, step_size, y, yp, t, p)
 end
 
 
-function solve_fdde_with_multiple_lags(FDDE::FDDEProblem, step_size)
-    @unpack f, h, order, constant_lags, p, tspan = FDDE
+function solve_fdde_with_multiple_lags(prob::FDDEProblem, step_size)
+    @unpack f, h, order, u0, constant_lags, p, tspan = prob
     τ = constant_lags[1]
-    t = collect(0:step_size:tspan[2])
-    maxn = length(t)
-    yp = zeros(Float64, maxn)
+    t = collect(tspan[1]:step_size:tspan[2])
+    N = length(t)
+    yp = [similar(u0) for i in 1:N]
+    yp = similar(yp)
     y = copy(t)
-    y[1] = h(p, 0)
+    y[1] = [h(p, 0) for i in 1:length(u0)]
 
-    for n in 1:maxn-1
+    for n in 1:N-1
         yp[n+1] = 0
         for j = 1:n
             yp[n+1] = yp[n+1]+generalized_binomials(j-1, n-1, order, step_size)*f(t[j], y[j], multiv(h, j, τ, step_size, y, yp, p)...)
@@ -186,14 +189,15 @@ function solve_fdde_with_multiple_lags(FDDE::FDDEProblem, step_size)
     end
 
     V = []
-    for n = 1:maxn
+    for n = 1:N
         push!(V, multiv(h, n, τ, step_size, y, yp, p))
     end
-
+#=
     delayed = zeros(length(τ), length(V))
     for i=1:length(V)
         delayed[:, i] = V[i]
     end    
+    =#
     return delayed, y
 end
 
@@ -242,12 +246,12 @@ function solve_fdde_with_single_lag_and_variable_order(FDDE::FDDEProblem, step_s
     τ = constant_lags[1]
     T = tspan[2]
     t = collect(0:step_size:T)
-    maxn = size(t, 1)
+    N = size(t, 1)
     yp = collect(0:step_size:T+step_size)
     y = copy(t)
     y[1] = h(p, 0)
 
-    for n in 1:maxn-1
+    for n in 1:N-1
         yp[n+1] = 0
         for j = 1:n
             if iip
@@ -282,7 +286,7 @@ function solve_fdde_with_single_lag_and_variable_order(FDDE::FDDEProblem, step_s
     end
 
     V = copy(t)
-    @fastmath @inbounds @simd for n = 1:maxn-1
+    @fastmath @inbounds @simd for n = 1:N-1
         V[n] = v(h, n, τ, step_size, y, yp, t, p)
     end
     return V, y
@@ -293,12 +297,12 @@ function solve_fdde_with_multiple_lags_and_variable_order(FDDE::FDDEProblem, ste
     @unpack f, h, order, constant_lags, p, tspan = FDDE
     τ = constant_lags[1]
     t = collect(0:step_size:tspan[2])
-    maxn = length(t)
-    yp = zeros(maxn)
+    N = length(t)
+    yp = zeros(N)
     y = copy(t)
     y[1] = h(p, 0)
 
-    for n in 1:maxn-1
+    for n in 1:N-1
         yp[n+1] = 0
         for j = 1:n
             yp[n+1] = yp[n+1]+generalized_binomials(j-1, n-1, order(t[n+1]), step_size)*f(t[j], y[j], multiv(h, j, τ, step_size, y, yp, p)...)
@@ -313,15 +317,16 @@ function solve_fdde_with_multiple_lags_and_variable_order(FDDE::FDDEProblem, ste
 
         y[n+1] = y[n+1]/gamma(order(t[n+1]))+step_size^order(t[n+1])*f(t[n+1], yp[n+1], multiv(h, n+1, τ, step_size, y, yp, p)...)/gamma(order(t[n+1])+2) + h(p, 0)
     end
-
+#=
     V = []
-    for n = 1:maxn
+    for n = 1:N
         push!(V, multiv(h, n, τ, step_size, y, yp, p))
     end
 
     delayed = zeros(Float, length(τ), length(V))
     for i in eachindex(V)
         delayed[:, i] = V[i]
     end    
+    =#
     return delayed, y
 end

src/delay/product_integral.jl

@@ -1,25 +1,3 @@
-#=
-# Usage
-
-    solve(FDDE::FDDEProblem, dt, DelayPI())
-
-Use explicit rectangular product integration algorithm to solve an FDDE problem.
-
-### References
-
-```tex
-@article{2020,
-   title={On initial conditions for fractional delay differential equations},
-   ISSN={1007-5704},
-   url={http://dx.doi.org/10.1016/j.cnsns.2020.105359},
-   DOI={10.1016/j.cnsns.2020.105359},
-   journal={Communications in Nonlinear Science and Numerical Simulation},
-   author={Garrappa, Roberto and Kaslik, Eva},
-   year={2020},
-}
-```
-=#
-
 function solve(FDDE::FDDEProblem, dt, ::PIEX)
     @unpack f, order, h, tspan, p, constant_lags = FDDE
     τ = constant_lags[1]

src/types/algorithms.jl

@@ -185,3 +185,23 @@ struct Euler <: FODESystemAlgorithm end
 Solve Atangana-Baleanu fractional order differential equations using Newton Polynomials.
 """
 struct AtanganaSedaAB <: FODESystemAlgorithm end
+
+
+###################### FDDE ######################
+
+"""
+    solve(FDDE::FDDEProblem, dt, DelayABM())
+
+Use the Adams-Bashforth-Moulton method to solve fractional delayed differential equations.
+
+### References
+
+```tex
+@inproceedings{Bhalekar2011APS,
+  title={A PREDICTOR-CORRECTOR SCHEME FOR SOLVING NONLINEAR DELAY DIFFERENTIAL EQUATIONS OF FRACTIONAL ORDER},
+  author={Sachin Bhalekar and Varsha Daftardar-Gejji},
+  year={2011}
+}
+```
+"""
+struct DelayABM <: FDDEAlgorithm end

src/types/problems.jl

@@ -194,7 +194,7 @@ TruncatedStacktraces.@truncate_stacktrace FDDEProblem 5 1 2
 FDDEProblem(f, args...; kwargs...) = FDDEProblem(DDEFunction(f), args...; kwargs...)
 
 function FDDEProblem(f::SciMLBase.AbstractDDEFunction, args...; kwargs...)
-    FDDEProblem{SciMLBase.isinplace(f)}(f, args...; kwargs...)
+    FDDEProblem{SciMLBase.isinplace(f, 5)}(f, args...; kwargs...)
 end
 
 """

test/fdde.jl

@@ -7,7 +7,7 @@
         end
     end
 
-    f(y, ϕ, p, t) = 3.5*y*(1-ϕ/19)
+    f(y, ϕ, p, t) =@. 3.5*y*(1-ϕ/19)
 
     dt = 0.5
     α = 0.97
@@ -112,32 +112,10 @@ end
         dy[4] = y[3]-0.5*y[4]
     end
     q = [60, 10, 10, 20]; α = [0.95, 0.95, 0.95, 0.95]
-    prob = FDDESystem(EnzymeKinetics!, q, α, 4, 6)
+    prob = FDDEProblem(EnzymeKinetics!, q, α, 4, 6)
     #sold, sol = solve(prob, 0.1, DelayABM())
     sol = solve(prob, 0.1, DelayABM())
-#=
-    @test isapprox(sold, [ 56.3     11.3272   11.4284  22.4025
-    56.3229  11.2293   11.4106  22.4194
-    56.3468  11.141    11.3849  22.4332
-    56.3731  11.0592   11.3533  22.4434
-    56.4026  10.9824   11.3171  22.4495
-    56.4358  10.9096   11.2773  22.4516
-    56.4729  10.8398   11.2346  22.4493
-    56.5144  10.7725   11.1896  22.4429
-    56.5604  10.707    11.1429  22.4322
-    56.6111  10.6431   11.0946  22.4175
-    56.6668  10.5803   11.0452  22.3987
-    56.7275  10.5186   10.9947  22.376
-    56.7932  10.4577   10.9435  22.3495
-    56.864   10.3975   10.8916  22.3194
-    56.9398  10.3382   10.8392  22.2857
-    57.0205  10.2795   10.7863  22.2488
-    57.1059  10.2217   10.7332  22.2086
-    57.1961  10.1647   10.6799  22.1653
-    57.2908  10.1087   10.6265  22.1191
-    57.3897  10.0537   10.5731  22.0702
-    57.4928   9.99987  10.5198  22.0187]; atol=1e-3)
-=#
+
     @test isapprox(sol.u, [60.5329  10.8225  10.6355  20.6245
     60.3612  10.9555  10.6638  20.6616
     60.1978  11.0679  10.7015  20.6988