Merge pull request #94 from ErikQQY/qqy/refactor_fode_sol

Refactor FODEProblem solutions

src/FractionalDiffEq.jl

@@ -1,7 +1,8 @@
 module FractionalDiffEq
 
 using LinearAlgebra
-using SciMLBase, DiffEqBase
+using SciMLBase
+using DiffEqBase
 using SpecialFunctions
 using SparseArrays
 using InvertedIndices: Not
@@ -37,7 +38,6 @@ include("fodesystem/bdf.jl")
 include("fodesystem/newton_gregory.jl")
 include("fodesystem/trapezoid.jl")
 include("fodesystem/explicit_pi.jl")
-include("fodesystem/Euler.jl")
 include("fodesystem/GLWithMemory.jl")
 include("fodesystem/NonLinear.jl")
 include("fodesystem/newton_polynomials.jl")
@@ -93,7 +93,6 @@ export FODESystemSolution, FDDESystemSolution, FFMODESystem
 export PIPECE, PIRect, PITrap
 export PECE, FODEMatrixDiscrete, ClosedForm, ClosedFormHankelM, GL
 export AtanganaSedaAB
-#export Euler
 
 # System of FODE solvers
 export NonLinearAlg, BDF, NewtonGregory, Trapezoid, PIEX, NewtonPolynomial

src/fodesystem/GLWithMemory.jl

@@ -1,4 +1,4 @@
-function solve(prob::FODEProblem, ::GL; dt = 0.0)
+function solve(prob::FODEProblem, alg::GL; dt = 0.0)
     dt ≤ 0 ? throw(ArgumentError("dt must be positive")) : nothing
     # GL method is only for same order FODE
     @unpack f, order, u0, tspan, p = prob
@@ -35,5 +35,7 @@ function solve(prob::FODEProblem, ::GL; dt = 0.0)
         f(du, result[:, k-1], p, t0+(k-1)*dt)
         result[:, k] = @. horder*du-summation
     end
-    return FODESystemSolution(t, result)
+    u = collect(Vector{eltype(u0)}, eachcol(result))
+
+    return DiffEqBase.build_solution(prob, alg, t, u)
 end

src/fodesystem/NonLinear.jl

@@ -1,4 +1,4 @@
-function solve(prob::FODEProblem, ::NonLinearAlg; dt = 0.0, L0=1e10)
+function solve(prob::FODEProblem, alg::NonLinearAlg; dt = 0.0, L0=1e10)
     dt ≤ 0 ? throw(ArgumentError("dt must be positive")) : nothing
     @unpack f, order, u0, tspan, p = prob
     t0 = tspan[1]; T = tspan[2]
@@ -32,7 +32,9 @@ function solve(prob::FODEProblem, ::NonLinearAlg; dt = 0.0, L0=1e10)
         z[:, k] = x1 - u0
     end
     result = z + repeat(u0, 1, m)
-    return FODESystemSolution(time, result)
+    u = collect(Vector{eltype(result)}, eachcol(result))
+
+    return DiffEqBase.build_solution(prob, alg, time, u)
 end
 
 """

src/fodesystem/PIPECE.jl

@@ -12,7 +12,7 @@ mutable struct M
     bn_fft
 end
 
-function solve(prob::FODEProblem, ::PECE; dt = 0.0)
+function solve(prob::FODEProblem, alg::PECE; dt = 0.0)
     dt ≤ 0 ? throw(ArgumentError("dt must be positive")) : nothing
     @unpack f, order, u0, tspan, p = prob
     t0 = tspan[1]; T = tspan[2]
@@ -139,7 +139,9 @@ function solve(prob::FODEProblem, ::PECE; dt = 0.0)
     end
 
     t = t[1:N+1]; y = y[:, 1:N+1]
-    return FODESystemSolution(t, y)
+    u = collect(Vector{eltype(u0)}, eachcol(y))
+
+    return DiffEqBase.build_solution(prob, alg, t, u)
 
 end
 

src/fodesystem/atangana_seda.jl

@@ -1,4 +1,4 @@
-function solve(prob::FODEProblem, ::AtanganaSedaAB; dt = 0.0)
+function solve(prob::FODEProblem, alg::AtanganaSedaAB; dt = 0.0)
     dt ≤ 0 ? throw(ArgumentError("dt must be positive")) : nothing
     @unpack f, order, u0, tspan, p = prob
     order = order[1]
@@ -41,7 +41,9 @@ function solve(prob::FODEProblem, ::AtanganaSedaAB; dt = 0.0)
 
         result[:, n+1] = u0+(1-order)/AB*temptemptemp+((dt^order)*order/(AB*gamma(order+1)))*temptemp1 + ((dt^order)*order/(AB*gamma(order+2)))*temptemp2+((dt^order)*order/(2*AB*gamma(order+3)))*temptemp3
     end
-    return FODESystemSolution(t, result)
+    u = collect(Vector{eltype(u0)}, eachcol(result))
+
+    return DiffEqBase.build_solution(prob, alg, t, u)
 end
 
 

src/fodesystem/bdf.jl

@@ -1,4 +1,4 @@
-function solve(prob::FODEProblem, ::BDF; dt=0.0, reltol=1e-6, abstol=1e-6, maxiters = 100)
+function solve(prob::FODEProblem, alg::BDF; dt=0.0, reltol=1e-6, abstol=1e-6, maxiters = 100)
     dt ≤ 0 ? throw(ArgumentError("dt must be positive")) : nothing
     @unpack f, order, u0, tspan, p = prob
     t0 = tspan[1]; tfinal = tspan[2]
@@ -61,7 +61,9 @@ function solve(prob::FODEProblem, ::BDF; dt=0.0, reltol=1e-6, abstol=1e-6, maxit
         y[:, N+1] = (1-c)*y[:, N] + c*y[:, N+1]
     end
     t = t[1:N+1]; y = y[:, 1:N+1]
-    return FODESystemSolution(t, y)
+    y = collect(Vector{eltype(y)}, eachcol(y))
+
+    return DiffEqBase.build_solution(prob, alg, t, y)
 end
 
 

src/fodesystem/explicit_pi.jl

@@ -1,4 +1,4 @@
-function solve(prob::FODEProblem, ::PIEX; dt = 0.0)
+function solve(prob::FODEProblem, alg::PIEX; dt = 0.0)
     dt ≤ 0 ? throw(ArgumentError("dt must be positive")) : nothing
     @unpack f, order, u0, tspan, p = prob
 
@@ -114,7 +114,9 @@ function solve(prob::FODEProblem, ::PIEX; dt = 0.0)
     end
 
     t = t[1:N+1]; y = y[:, 1:N+1]
-    return FODESystemSolution(t, y)
+    u = collect(Vector{eltype(u0)}, eachcol(y))
+
+    return DiffEqBase.build_solution(prob, alg, t, u)
 end
 
 

src/fodesystem/newton_gregory.jl

@@ -1,4 +1,4 @@
-function solve(prob::FODEProblem, ::NewtonGregory; dt = 0.0, reltol=1e-6, abstol=1e-6, maxiters = 100)
+function solve(prob::FODEProblem, alg::NewtonGregory; dt = 0.0, reltol=1e-6, abstol=1e-6, maxiters = 100)
     dt ≤ 0 ? throw(ArgumentError("dt must be positive")) : nothing
     @unpack f, order, u0, tspan, p = prob
     t0 = tspan[1]; tfinal = tspan[2]
@@ -63,7 +63,9 @@ function solve(prob::FODEProblem, ::NewtonGregory; dt = 0.0, reltol=1e-6, abstol
         y[:, N+1] = (1-c)*y[:, N] + c*y[:, N+1]
     end
     t = t[1:N+1]; y = y[:, 1:N+1]
-    return FODESystemSolution(t, y)
+    u = collect(Vector{eltype(u0)}, eachcol(y))
+
+    return DiffEqBase.build_solution(prob, alg, t, u)
 end
 
 

src/fodesystem/newton_polynomials.jl

@@ -1,4 +1,4 @@
-function solve(prob::FODEProblem, ::NewtonPolynomial; dt = 0.0)
+function solve(prob::FODEProblem, alg::NewtonPolynomial; dt = 0.0)
     dt ≤ 0 ? throw(ArgumentError("dt must be positive")) : nothing
     @unpack f, order, u0, tspan, p = prob
     t0 = tspan[1]; tfinal = tspan[2]
@@ -27,5 +27,7 @@ function solve(prob::FODEProblem, ::NewtonPolynomial; dt = 0.0)
         f(temp3, result[:, n-2], p, t[n-2])
         result[:, n+1] = result[:, n] + (1-order)/M*(temp1-temp2)+order.*M.*dt.*(23/12*temp1-4/3*temp2+5/12*temp3)
     end
-    return FODESystemSolution(t, result)
+    u = collect(Vector{eltype(u0)}, eachcol(result))
+
+    return DiffEqBase.build_solution(prob, alg, t, u)
 end

src/fodesystem/trapezoid.jl

@@ -1,4 +1,4 @@
-function solve(prob::FODEProblem, ::Trapezoid; dt = 0.0, reltol=1e-6, abstol=1e-6)
+function solve(prob::FODEProblem, alg::Trapezoid; dt = 0.0, reltol=1e-6, abstol=1e-6)
     dt ≤ 0 ? throw(ArgumentError("dt must be positive")) : nothing
     @unpack f, order, u0, tspan, p = prob
     t0 = tspan[1]; tfinal = tspan[2]
@@ -64,7 +64,9 @@ function solve(prob::FODEProblem, ::Trapezoid; dt = 0.0, reltol=1e-6, abstol=1e-
         y[:, N+1] = (1-c)*y[:, N] + c*y[:, N+1]
     end
     t = t[1:N+1]; y = y[:, 1:N+1]
-    return FODESystemSolution(t, y)
+    u = collect(Vector{eltype(u0)}, eachcol(y))
+
+    return DiffEqBase.build_solution(prob, alg, t, u)
 end
 
 

src/lyapunov.jl

@@ -87,9 +87,8 @@ function noncommensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out
         # Here we directly use the buildin PECE algorithm to solve the extend system, SO FAST!!!
         prob = FODEProblem(extend_fun, q[:], x[:], (t, t+h_norm), p)
         sol = solve(prob, PECE(), dt = h)
-        Y = sol.u
+        Y = mapreduce(permutedims, vcat, sol.u)
         t = t+h_norm
-        Y = Y'
 
         x = Y[size(Y, 1), :]
         for i=1:ne

test/FODETests.jl → test/fode.jl

@@ -1,30 +1,4 @@
-#=
-@testset "Test Diethelm PECE algorithms" begin
-    fun(y, p, t) = 1-y
-    prob = SingleTermFODEProblem(fun, 1.8, [0, 0], (0, 5))
-    sol = solve(prob, 0.5, PECE())
-
-    @test isapprox(sol.u, [0.16153482345602124
-    0.5289988135096678
-    0.848191872231344
-    1.0899386876765813
-    1.2276684080819034
-    1.2684889407947215
-    1.2385617809129195
-    1.1700696107638846
-    1.0918355406175286
-    1.0242165055531318
-    0.977769035221195]; atol=1e-3)
-end
 
-@testset "Test forward Euler method" begin
-    fun(y, p, t) = 1-y
-    prob = SingleTermFODEProblem(fun, 0.5, 0, (0, 5))
-    sol = solve(prob, 0.5, Euler())
-
-    @test isapprox(sol.u, [0.0, 0.7978845608028654, 0.4917593947279313, 0.7259128254126388, 0.6517091834824043, 0.7289344741198424, 0.7179084631988641, 0.7483267686842404, 0.7528156154229637, 0.7681082183772983, 0.7753604356669395]; atol=1e-3)
-end
-=#
 @testset "Test Matrix discrete method" begin
     fun(x, y) = 1-y
     u0 = [0, 0]
@@ -86,7 +60,7 @@ end
 
 @testset "Test GL method for FODEProblem" begin
     alpha = [0.99, 0.99, 0.99]
-    u0 = [1, 0, 1]
+    u0 = [1.0, 0.0, 1.0]
     tspan = (0.0, 1.0)
     function testf!(du, u, p, t)
         a, b, c = 10, 28, 8/3
@@ -96,7 +70,7 @@ end
     end
     prob = FODEProblem(testf!, alpha, u0, tspan)
     sol = solve(prob, GL(), dt=0.5)
-    @test isapprox(sol.u, [1.0  -4.04478    84.8074
+    @test isapprox(test_sol(sol), [1.0  -4.04478    84.8074
     0.0  13.5939    -51.1251
     1.0  -0.352607  -27.5541]; atol=1e-4)
 end
@@ -115,9 +89,9 @@ end
     u0 = [0.2; -0.1; 0.1];
     tspan = (0, 0.5);
     prob = FODEProblem(chua!, α, u0, tspan)
-    result = solve(prob, NonLinearAlg(), dt = 0.1)
+    sol = solve(prob, NonLinearAlg(), dt = 0.1)
 
-    @test isapprox(result.u, [  0.2   0.11749     0.074388  0.0733938  0.117483   0.210073
+    @test isapprox(test_sol(sol), [  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
@@ -357,7 +331,7 @@ end
     prob = FODEProblem(test, alpha, y0, tspan)
     sol = solve(prob, PECE(), dt=0.01)
 
-    @test isapprox(sol.u, [1.2  1.2061   1.21042  1.21421  1.21767  1.22089  1.22392  1.22678  1.2295   1.2321   1.23459
+    @test isapprox(test_sol(sol), [1.2  1.2061   1.21042  1.21421  1.21767  1.22089  1.22392  1.22678  1.2295   1.2321   1.23459
     2.8  2.78118  2.76953  2.7596   2.75064  2.74235  2.73455  2.72715  2.72008  2.71329  2.70673]; atol=1e-3)
 end
 
@@ -374,7 +348,7 @@ end
     prob = FODEProblem(test, alpha, y0, tspan)
     sol = solve(prob, PIEX(), dt = 0.015)
 
-    @test isapprox(sol.u, [1.2  1.20865  1.21458  1.21975  1.22443  1.22874  1.23276  1.23652  1.24006  1.24339  1.24653  1.2495   1.25229  1.25493  1.25575
+    @test isapprox(test_sol(sol), [1.2  1.20865  1.21458  1.21975  1.22443  1.22874  1.23276  1.23652  1.24006  1.24339  1.24653  1.2495   1.25229  1.25493  1.25575
     2.8  2.77486  2.75941  2.74621  2.73429  2.72326  2.71291  2.70309  2.69373  2.68477  2.67616  2.66786  2.65986  2.65213  2.64964]; atol=1e-4)
 end
 
@@ -390,7 +364,7 @@ end
     prob = FODEProblem(brusselator, alpha, u0, tspan)
     sol = solve(prob, NewtonGregory(), dt = 0.1)
 
-    @test isapprox(sol.u, [ 0.2   0.200531  0.201161  0.201809  0.202453  0.203089
+    @test isapprox(test_sol(sol), [ 0.2   0.200531  0.201161  0.201809  0.202453  0.203089
     0.03  0.165554  0.265678  0.355635  0.439544  0.519211]; atol=1e-4)
 end
 
@@ -406,7 +380,7 @@ end
     prob = FODEProblem(brusselator, alpha, u0, tspan)
     sol = solve(prob, BDF(), dt = 0.1)
 
-    @test isapprox(sol.u, [ 0.2   0.200531  0.201161  0.201809  0.202453  0.203089
+    @test isapprox(test_sol(sol), [ 0.2   0.200531  0.201161  0.201809  0.202453  0.203089
     0.03  0.165554  0.265678  0.355635  0.439544  0.519211]; atol=1e-4)
 end
 
@@ -422,25 +396,26 @@ end
     prob = FODEProblem(brusselator, alpha, u0, tspan)
     sol = solve(prob, Trapezoid(), dt = 0.1)
 
-    @test isapprox(sol.u, [0.2   0.200531  0.201161  0.201808  0.202452  0.203088
+    @test isapprox(test_sol(sol), [0.2   0.200531  0.201161  0.201808  0.202452  0.203088
     0.03  0.165554  0.265678  0.355635  0.439545  0.519211]; atol=1e-4)
 end
 
 
 @testset "Test Newton Polynomial" begin
     t0=0; tfinal=5
+    tspan = (0.0, 5.0)
     α = [0.98, 0.98, 0.98]
-    u0 = [-1, 1, 1]
+    u0 = [-1.0, 1.0, 1.0]
     function fun(du, u, p, t)
         b=0.1;
         du[1] = -sin(u[2])-b*u[1]
         du[2] = -sin(u[3])-b*u[2]
         du[3] = -sin(u[1])-b*u[3]
     end
-    prob = FODEProblem(fun, α, u0, (t0, tfinal))
+    prob = FODEProblem(fun, α, u0, tspan)
     sol = solve(prob, NewtonPolynomial(), dt = 0.5)
 
-    @test isapprox(sol.u, [-1.0  -1.37074   -1.46124    -1.21771   -0.884241  -0.49327   -0.0897447   0.338635   0.793532   1.2323    1.55187
+    @test isapprox(test_sol(sol), [-1.0  -1.37074   -1.46124    -1.21771   -0.884241  -0.49327   -0.0897447   0.338635   0.793532   1.2323    1.55187
     1.0   0.529265  -0.0101035  -0.430956  -0.733314  -0.916321  -1.06158    -1.2647    -1.5767    -1.99613  -2.37772
     1.0   1.37074    1.8176      2.17624    2.48899    2.67047    2.6624      2.46322    2.07376    1.55057   1.0049]; atol=1e-4)
 end
@@ -478,7 +453,7 @@ end
     prob = FODEProblem(fun, α, u0, tspan)
     sol = solve(prob, AtanganaSedaAB(), dt = 0.5)
 
-    @test isapprox(sol.u, [ -0.1   0.7    1.18511  -1.41522  -31.0872
+    @test isapprox(test_sol(sol), [ -0.1   0.7    1.18511  -1.41522  -31.0872
     0.1   0.525  1.08088  -4.03383   32.0085
    -0.1  -0.065  1.03463   4.10537   13.3441]; atol=1e-2)
 end

test/runtests.jl

@@ -1,12 +1,16 @@
 using FractionalDiffEq, SpecialFunctions
 using Test
 
+function test_sol(sol)
+    return transpose(mapreduce(permutedims, vcat, sol.u))
+end
+
 @testset "FractionalDiffEq.jl" begin
-    include("FODETests.jl")
-    include("FDDETests.jl")
-    include("FDiscreteTests.jl")
-    include("DODETests.jl")
-    include("FFODETests.jl")
+    include("fode.jl")
+    include("fdde.jl")
+    include("discrete.jl")
+    include("dode.jl")
+    include("ffode.jl")
     include("auxillary.jl")
-    include("FOLEtest.jl")
+    include("fole.jl")
 end