Fix plot recipe for FODESystemSolution

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

README.md

@@ -107,20 +107,18 @@ By using the ```NonLinearAlg``` algorithms to solve this problem:
 ```julia
 using FractionalDiffEq, Plots
 function chua!(du, x, p, t)
- a = 10.725; b = 10.593
- c = 0.268
- m0 = -1.1726
- m1 = -0.7872
+ a, b, c, m0, m1 = p
  du[1] = a*(x[2]-x[1]-(m1*x[1]+0.5*(m0-m1)*(abs(x[1]+1)-abs(x[1]-1))))
  du[2] = x[1]-x[2]+x[3]
  du[3] = -b*x[2]-c*x[3]
 end
 α = [0.93, 0.99, 0.92];
 x0 = [0.2; -0.1; 0.1];
 h = 0.01; tspan = (0, 100);
-prob = FODESystem(chua!, α, x0, tspan)
-result = solve(prob, h, NonLinearAlg())
-plot(result[:, 1], result[:, 2], title="Chua System", legend=:bottomright)
+p = [10.725, 10.593, 0.268, -1.1726, -0.7872]
+prob = FODESystem(chua!, α, x0, tspan, p)
+sol = solve(prob, h, NonLinearAlg())
+plot(sol, vars=(1, 2), title="Chua System", legend=:bottomright)
 ```
 
 And plot the result:

docs/src/fdesystem.md

@@ -28,17 +28,18 @@ Use the ```NonLinearAlg``` algorithm in FractionalDiffEq.jl to solve the Chua sy
 ```julia
 using FractionalDiffEq, Plots
 function chua!(du, x, p, t)
- a = 10.725; b = 10.593; c = 0.268; m0 = -1.1726; m1 = -0.7872;
+ a, b, c, m0, m1 = p
  du[1] = a*(x[2]-x[1]-(m1*x[1]+0.5*(m0-m1)*(abs(x[1]+1)-abs(x[1]-1))))
  du[2] = x[1]-x[2]+x[3]
  du[3] = -b*x[2]-c*x[3]
 end
 α = [0.93, 0.99, 0.92];
 x0 = [0.2; -0.1; 0.1];
 h = 0.01; tspan = (0, 50);
-prob = FODESystem(chua!, α, x0, tspan)
-result = solve(prob, h, NonLinearAlg())
-plot(result[:, 1], result[:, 2], title="Chua System", legend=:bottomright)
+p = [10.725, 10.593, 0.268, -1.1726, -0.7872]
+prob = FODESystem(chua!, α, x0, tspan, p)
+sol = solve(prob, h, NonLinearAlg())
+plot(sol, vars=(1, 2) title="Chua System", legend=:bottomright)
 ```
 
 ![Chua](./assets/chua.png)

examples/chua.jl

@@ -1,17 +1,15 @@
 using FractionalDiffEq, Plots
 function chua!(du, x, p, t)
- a = 10.725; b = 10.593
- c = 0.268
- m0 = -1.1726
- m1 = -0.7872
+ a, b, c, m0, m1 = p
  du[1] = a*(x[2]-x[1]-(m1*x[1]+0.5*(m0-m1)*(abs(x[1]+1)-abs(x[1]-1))))
  du[2] = x[1]-x[2]+x[3]
  du[3] = -b*x[2]-c*x[3]
 end
 α = [0.93, 0.99, 0.92]
 x0 = [0.2; -0.1; 0.1]
-h = 0.001; tspan = (0, 20)
-prob = FODESystem(chua!, α, x0, tspan)
-result = solve(prob, h, NonLinearAlg())
+h = 0.001; tspan = (0, 50)
+p = [10.725, 10.593, 0.268, -1.1726, -0.7872]
+prob = FODESystem(chua!, α, x0, tspan, p)
+sol = solve(prob, h, NonLinearAlg())
 
-plot(result[1, :], result[2, :], title="Chua System", legend=:bottomright)
+plot(sol, vars=(1, 2), title="Fractional Chua System", legend=:bottomright)

src/fodesystem/NonLinear.jl

@@ -14,6 +14,7 @@ struct NonLinearAlg <: AbstractFDEAlgorithm end
 function solve(prob::FODESystem, h, ::NonLinearAlg, L0=1e10)
  @unpack f, α, u0, tspan, p = prob
  t0 = tspan[1]; T = tspan[2]
+ time = collect(t0:h:T)
  n = length(u0)
  m::Int = round(Int, (T-t0)/h)+1
  g = genfun(1)
@@ -44,7 +45,7 @@ function solve(prob::FODESystem, h, ::NonLinearAlg, L0=1e10)
  z[:, k] = x1 - u0
  end
  result = z + repeat(u0, 1, m)
- return result
+ return FODESystemSolution(time, result)
 end
 
 """

src/singletermfode/PECE.jl

@@ -241,6 +241,10 @@ struct FODESolution <: AbstractFDESolution
  u::AbstractArray
 end
 
+struct FODESystemSolution <: AbstractFDESolution
+ t::AbstractArray
+ u::AbstractArray
+end
 
 struct FDifferenceSolution <: AbstractFDESolution
  t::AbstractArray

src/utils.jl

@@ -113,4 +113,57 @@ Fractional differential equation solutions visulization hooks.
 
 @recipe f(sol::DODESolution) = sol.t, sol.u
 
-@recipe f(sol::FFMODESolution) = sol.t, sol.u
+@recipe f(sol::FFMODESolution) = sol.t, sol.u
+
+@recipe function f(sol::FODESystemSolution; 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
+ 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
+ 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
+ sol.u[index[1], :], sol.u[index[2], :]
+ else
+ sol.t, sol.u[index[index0[1]], :]
+ end
+ else
+ error("Plot variable index is not correct.")
+ end
+=#
+ end
+end

test/FODETests.jl

@@ -134,7 +134,7 @@ end
  prob = FODESystem(chua!, α, x0, tspan)
  result = solve(prob, h, NonLinearAlg())
 
- @test isapprox(result, [ 0.2 0.11749 0.074388 0.0733938 0.117483 0.210073
+ @test isapprox(result.u, [ 0.2 0.11749 0.074388 0.0733938 0.117483 0.210073
  -0.1 -0.0590683 -0.018475 0.0192931 0.0534393 0.0844175
  0.1 0.224134 0.282208 0.286636 0.246248 0.168693]; atol=1e-3)
 end