Refactor FOLyapunov for intuitive usage

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

docs/make.jl

@@ -18,7 +18,7 @@ makedocs(;
         "Get Started" => "get_started.md",
         "Root Locus" => "rlocus.md",
         "Frequency Analysis" => "frequency.md",
-        "Lyapunov Exponents" => "flyapunovexp.md",
+        "Lyapunov Exponents Spectrum" => "flyapunovexp.md",
         "Conversation" => "conversation.md"
     ],
 )

docs/src/flyapunovexp.md

@@ -18,7 +18,6 @@ function RF(du, u, t)
     du[1] = u[2]*(u[3]-1+u[1]*u[1])+0.1*u[1];
     du[2] = u[1]*(3*u[3]+1-u[1]*u[1])+0.1*u[2];
     du[3] = -2*u[3]*(0.98+u[1]*u[2]);
-    du
 end
 LE, tspan = FOLyapunov(RF, 0.98, 0, 0.02, 300, [0.1; 0.1; 0.1], 0.005, 1000)
 ```
@@ -46,7 +45,7 @@ The output would be:
 The computed **LE** is the Lyapunov exponent of this system.
 
 ```julia-repl
-julia> LE
+julia> LE[end-2:end]
 3×1 Matrix{Float64}:
   0.06111650166568285
   0.0038981396237095034
@@ -84,7 +83,6 @@ function Danca(du, u, t)
     du[2] = -u[3]*sign(u[1])+u[4]
     du[3] = abs(u[1])-1
     du[4] = -0.5*u[2]
-    du
 end
 
 LE, tspan=FOLyapunov(Danca, 0.98, 0, 0.02, 300, [0.1; 0.1; 0.1; 0.1], 0.005, 1000)
@@ -95,6 +93,6 @@ plot!(tspan, LE[3, :])
 plot!(tspan, LE[4, :])
 ```
 
-By plot the Lyapunov exponent:
+By plot the Lyapunov exponent spectrum:
 
 ![PWC](./assets/PWCLE.png)

examples/PWCLE.jl

@@ -1,11 +1,10 @@
 using FractionalSystems, Plots
 
-function Danca(du, u, t)
+function Danca(du, u, p, t)
     du[1] = -u[1]+u[2]
     du[2] = -u[3]*sign(u[1])+u[4]
     du[3] = abs(u[1])-1
     du[4] = -0.5*u[2]
-    du
 end
 
 LE, tspan=FOLyapunov(Danca, 0.98, 0, 0.02, 300, [0.1; 0.1; 0.1; 0.1], 0.005, 1000)
@@ -14,4 +13,4 @@ LE, tspan=FOLyapunov(Danca, 0.98, 0, 0.02, 300, [0.1; 0.1; 0.1; 0.1], 0.005, 100
 plot(tspan, LE[1, :])
 plot!(tspan, LE[2, :])
 plot!(tspan, LE[3, :])
-plot!(tspan, LE[4, :])
+plot!(tspan, LE[4, :])

examples/RFLE.jl

@@ -1,13 +1,12 @@
-using FractionalSystems
+using FractionalSystems, Plots
 
-function RF(du, u, t)
-    du[1] = u[2]*(u[3]-1+u[1]*u[1])+0.1*u[1];
-    du[2] = u[1]*(3*u[3]+1-u[1]*u[1])+0.1*u[2];
-    du[3] = -2*u[3]*(0.98+u[1]*u[2]);
-    du
+function RF(du, u, p, t)
+    du[1] = u[2]*(u[3]-1+u[1]*u[1])+0.1*u[1]
+    du[2] = u[1]*(3*u[3]+1-u[1]*u[1])+0.1*u[2]
+    du[3] = -2*u[3]*(0.98+u[1]*u[2])
 end
 
-LE, tspan=FOLyapunov(RF, 0.98, 0, 0.02, 300, [0.1; 0.1; 0.1], 0.005, 1000)
+LE, tspan=FOLyapunov(RF, 0.999, 0, 0.02, 300, [0.1; 0.1; 0.1], 0.005, 1000)
 
 plot(tspan, LE[1, :])
 plot!(tspan, LE[2, :])

src/FOLE.jl

@@ -25,7 +25,8 @@ function FOLyapunov(fun, order, t_start, h_norm, t_end, u0, h, out)# TODO: Gener
     # Generate extend system with jacobian
     function extend_fun(t, temp)
         temp=reshape(temp, ne, ne+1)
-        result = fun(zeros(ne), temp[:, 1], t)
+        result = zeros(ne)
+        fun(result, temp[:, 1], nothing, t)
         for i=2:ne+1
             result = [result; Jfdefun(t, temp[:, 1])*temp[:, i]]
         end
@@ -103,7 +104,7 @@ function FOLyapunov(fun, order, t_start, h_norm, t_end, u0, h, out)# TODO: Gener
     return LE, tspan
 end
 
-
+FOLyapunov(sys::FODESystem, h_norm, h, out) = FOLyapunov(sys.f, sys.α, tspan[1], h_norm, tspan[2], sys.x0, h, out)
 
 mutable struct M
     an
@@ -220,7 +221,7 @@ function pc(alpha, f_fun, t0, T, y0, h)
     # Initializing solution and proces of computation
     t = t0 .+ collect(0:N)*h
     y[:, 1] = y0[:, 1]
-    fy[:, 1] = f_temp ;
+    fy[:, 1] = f_temp
     (y, fy) = Triangolo(1, r-1, t, y, fy, zn_pred, zn_corr, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, y0, t0, f_fun) ;
 
     # Main process of computation by means of the FFT algorithm
@@ -435,7 +436,7 @@ end
 function jacobian_of_fdefun(f, t, y)
     ForwardDiff.jacobian(y) do y
     du = similar(y)
-    f(du, y, t)
+    f(du, y, nothing, t)
     du
     end
 end

src/FractionalSystems.jl

@@ -2,7 +2,7 @@ module FractionalSystems
 
 using ControlSystems, ControlSystemIdentification
 using LinearAlgebra, SpecialFunctions
-using FractionalDiffEq
+import FractionalDiffEq: FODESystem
 using ForwardDiff
 
 include("foss/foss.jl")

test/FOLE.jl

@@ -3,16 +3,15 @@ using Test
 
 @testset "Test Lyapunov exponents" begin
     function RF(du, u, t)
-        du[1] = u[2]*(u[3]-1+u[1]*u[1])+0.1*u[1];
-        du[2] = u[1]*(3*u[3]+1-u[1]*u[1])+0.1*u[2];
-        du[3] = -2*u[3]*(0.98+u[1]*u[2]);
-        du
+        du[1] = u[2]*(u[3]-1+u[1]*u[1])+0.1*u[1]
+        du[2] = u[1]*(3*u[3]+1-u[1]*u[1])+0.1*u[2]
+        du[3] = -2*u[3]*(0.98+u[1]*u[2])
     end
 
     (LE, tspan)=FOLyapunov(RF, 0.999, 0, 0.02, 300, [0.1; 0.1; 0.1], 0.005, 1000)
 
     # Test the latest computed(final) Lyapunov exponents
-    @test isapprox(LE[end-2:end], [  0.06111650166568285
+    @test isapprox(LE[end-2:end], [0.06111650166568285
     0.0038981396237095034
    -1.8324646820425692]; atol=1e-3)
 end