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]
@@ -56,12 +56,14 @@
  end
  # Evaluation solution in TFINAL when TFINAL is not in the mesh
  if tfinal < t[N+1]
  c = (tfinal - t[N])/dt
  t[N+1] = tfinal
  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]
@@ -58,12 +58,14 @@
  end
  # Evaluation solution in TFINAL when TFINAL is not in the mesh
  if tfinal < t[N+1]
  c = (tfinal - t[N])/dt
  t[N+1] = tfinal
  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