Fix plot recipe for FDDESystemSolution

Signed-off-by: ErikQQY <2283984853@qq.com>

docs/src/fdde.md

@@ -103,8 +103,8 @@ function delaychen!(dy, y, ϕ, t)
 	dy[3] = y[1]*y[2]-b*ϕ[3]
 end
 prob = FDDESystem(delaychen!, ϕ, α, τ, T)
-y, x=solve(prob, h, DelayABM())
-plot(x[:, 1], x[:, 2], x[:, 3], title="Fractional Order Chen Delayed System")
+sol=solve(prob, h, DelayABM())
+plot(sol, title="Fractional Order Chen Delayed System")
 ```
 
 ![FOChen](./assets/fodelaychen.png)
@@ -120,7 +120,7 @@ y(t)=[60, 10, 10, 20],\ t\leq0
 ```
 
 ```julia
-using FractionalDiffEq
+using FractionalDiffEq, Plots
 function EnzymeKinetics!(dy, y, ϕ, t)
     dy[1] = 10.5-y[1]/(1+0.0005*ϕ[4]^3)
     dy[2] = y[1]/(1+0.0005*ϕ[4]^3)-y[2]
@@ -129,13 +129,8 @@ function EnzymeKinetics!(dy, y, ϕ, t)
 end
 q = [60, 10, 10, 20]; α = [0.95, 0.95, 0.95, 0.95]
 prob = FDDESystem(EnzymeKinetics!, q, α, 4, 150)
-sold, sol = solve(prob, 0.01, DelayABM())
-tspan = collect(0:0.01:146)
-using Plots
-plot(tspan, sol[:, 1])
-plot!(tspan, sol[:, 2])
-plot!(tspan, sol[:, 3])
-plot!(tspan, sol[:, 4])
+sol = solve(prob, 0.01, DelayABM())
+plot(sol)
 ```
 
 ![Enzyme](./assets/enzyme_kinetics.png)

examples/delaychen.jl

@@ -7,5 +7,7 @@ function delaychen!(dy, y, ϕ, t)
 	dy[3] = y[1]*y[2]-b*ϕ[3]
 end
 prob = FDDESystem(delaychen!, ϕ, α, τ, T)
-y, x=solve(prob, h, DelayABM())
-plot(x[:, 1], x[:, 2], x[:, 3], title="Fractional Order Chen Delayed System")
+#y, x=solve(prob, h, DelayABM())
+#plot(x[:, 1], x[:, 2], x[:, 3], title="Fractional Order Chen Delayed System")
+sol=solve(prob, h, DelayABM())
+plot(sol)

examples/enzyme_kinetics.jl

@@ -1,16 +1,11 @@
-using FractionalDiffEq
+using FractionalDiffEq, Plots
 function EnzymeKinetics!(dy, y, ϕ, t)
     dy[1] = 10.5-y[1]/(1+0.0005*ϕ[4]^3)
     dy[2] = y[1]/(1+0.0005*ϕ[4]^3)-y[2]
     dy[3] = y[2]-y[3]
     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)
-sold, sol = solve(prob, 0.1, DelayABM())
-tspan = collect(0:0.1:2)
-using Plots
-plot(tspan, sol[:, 1])
-plot!(tspan, sol[:, 2])
-plot!(tspan, sol[:, 3])
-plot!(tspan, sol[:, 4])
+prob = FDDESystem(EnzymeKinetics!, q, α, 4, 50)
+sol = solve(prob, 0.01, DelayABM())
+plot(sol)

src/FDDE/DelayABMSystem.jl

@@ -16,6 +16,7 @@ function solve(FDDESys::FDDESystem, h, ::DelayABM)
     len = length(ϕ)
     N::Int = round(Int, T/h)
     Ndelay = round(Int, τ/h)
+    t = collect(Float64, 0:h:(T-τ))
     x1 = zeros(Ndelay+N+1, len)
     x = zeros(Ndelay+N+1, len)
     du = zeros(len)
@@ -78,5 +79,6 @@ function solve(FDDESys::FDDESystem, h, ::DelayABM)
        yresult[n-2*Ndelay, :] = x[n-Ndelay, :]
     end
 
-    return xresult, yresult
+    #return xresult, yresult
+    return FDDESystemSolution(t, yresult')
 end

src/FractionalDiffEq.jl

@@ -109,6 +109,7 @@ export FIEProblem
 
 export AbstractFDESolution
 export FODESolution, FIESolution, FDifferenceSolution, DODESolution, FFMODESolution
+export FODESystemSolution, FDDESystemSolution
 
 # FODE solvers
 export PIEX, PIPECE, PIIMRect, PIIMTrap

src/fodesystem/PIEX.jl

@@ -129,7 +129,6 @@ function solve(prob::FODESystem, h, ::PIEX)
 
     t = t[1:N+1]; y = y[:, 1:N+1]
     return FODESystemSolution(t, y)
-
 end
 
 

src/singletermfode/PECE.jl

@@ -246,6 +246,11 @@ struct FODESystemSolution <: AbstractFDESolution
     u::AbstractArray
 end
 
+struct FDDESystemSolution <: AbstractFDESolution
+    t::AbstractArray
+    u::AbstractArray
+end
+
 struct FDifferenceSolution <: AbstractFDESolution
     t::AbstractArray
     u::AbstractArray

src/utils.jl

@@ -151,19 +151,45 @@ Fractional differential equation solutions visulization hooks.
                 sol.t, sol.u[newindex[1], :]
             end
         end
-#=
-        if len == 3 # For 3 dimension
-            sol.u[index[1], :], sol.u[index[2], :], sol.u[index[3], :]
-        elseif len == 2 # For 2 dimension
-            
-            if length(index0) == 0
+    end
+end
+
+
+@recipe function f(sol::FDDESystemSolution; vars=nothing)
+    if typeof(vars) == Nothing # When vars is not specified, show time versus each variable.
+        l = size(sol.u, 1)
+        for i in 1:l
+            @series begin
+                sol.t, sol.u[i, :]
+            end
+        end
+    else
+        index = Int[]
+        for i in vars
+            append!(index, i)
+        end
+        len = length(index)
+        index0 = findall(x->x==0, index) # Find the 0 index, which is the time index.
+        if length(index0) == 0
+            if len == 3
+                sol.u[index[1], :], sol.u[index[2], :], sol.u[index[3], :]
+            elseif len == 2
                 sol.u[index[1], :], sol.u[index[2], :]
             else
-                sol.t, sol.u[index[index0[1]], :]
+                error("Plots variable index is not correct.")
+            end
+        elseif length(index0) == 1
+            newindex = deleteat!(index, index0)
+            newlen = length(newindex)
+            if newlen == 2
+                for i in newindex
+                    @series begin
+                        sol.t, sol.u[i, :]
+                    end
+                end
+            elseif newlen == 1
+                sol.t, sol.u[newindex[1], :]
             end
-        else
-            error("Plot variable index is not correct.")
         end
-=#
     end
 end

test/FDDETests.jl

@@ -105,8 +105,9 @@ end
     end
     q = [60, 10, 10, 20]; α = [0.95, 0.95, 0.95, 0.95]
     prob = FDDESystem(EnzymeKinetics!, q, α, 4, 6)
-    sold, sol = solve(prob, 0.1, DelayABM())
-
+    #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
@@ -128,8 +129,8 @@ end
     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, [60.5329  10.8225  10.6355  20.6245
+=#
+    @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
     60.0407  11.1632  10.7456  20.7379
@@ -149,7 +150,7 @@ end
     58.2489  11.6016  11.3331  21.4857
     58.1427  11.6013  11.357   21.5433
     58.0389  11.5992  11.3783  21.6001
-    57.9373  11.5953  11.3971  21.6557]; atol=1e-3)
+    57.9373  11.5953  11.3971  21.6557]'; atol=1e-3)
 end
 
 @testset "Test FDDE Matrix Form method" begin