diff --git a/src/FractionalDiffEq.jl b/src/FractionalDiffEq.jl index 34480217..d237c651 100644 --- a/src/FractionalDiffEq.jl +++ b/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 diff --git a/src/fodesystem/Euler.jl b/src/fodesystem/Euler.jl deleted file mode 100644 index 286c740c..00000000 --- a/src/fodesystem/Euler.jl +++ /dev/null @@ -1,21 +0,0 @@ -function solve(prob::FODEProblem, ::Euler; 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] - order = order[1] - t = collect(Float64, t0:dt:tfinal) - N::Int = ceil(Int, (tfinal-t0)/dt) - l = length(u0) - result = zeros(Float64, l, N+1) - result[:, 1] = u0 - for n=1:N - temp = zeros(Float64, l) - tmp = zeros(Float64, l) - for j=1:n - f(tmp, result[:, j], p, t[j]) - temp = ((n-j+1)^order-(n-j)^order)*tmp - end - result[:, n+1] = result[:, 1] + dt^order/gamma(order+1)*temp - end - return result -end \ No newline at end of file diff --git a/src/fodesystem/GLWithMemory.jl b/src/fodesystem/GLWithMemory.jl index 37dd203f..0c6c27b1 100644 --- a/src/fodesystem/GLWithMemory.jl +++ b/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 \ No newline at end of file diff --git a/src/fodesystem/NonLinear.jl b/src/fodesystem/NonLinear.jl index 42fa88bb..64fd11fe 100644 --- a/src/fodesystem/NonLinear.jl +++ b/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 """ diff --git a/src/fodesystem/PIPECE.jl b/src/fodesystem/PIPECE.jl index 668bd1f7..bc3fc37d 100644 --- a/src/fodesystem/PIPECE.jl +++ b/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 diff --git a/src/fodesystem/atangana_seda.jl b/src/fodesystem/atangana_seda.jl index 30c117f5..5c6015e7 100644 --- a/src/fodesystem/atangana_seda.jl +++ b/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 diff --git a/src/fodesystem/bdf.jl b/src/fodesystem/bdf.jl index 19ea12af..117f81a7 100644 --- a/src/fodesystem/bdf.jl +++ b/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 diff --git a/src/fodesystem/explicit_pi.jl b/src/fodesystem/explicit_pi.jl index 332518e9..4a52c79d 100644 --- a/src/fodesystem/explicit_pi.jl +++ b/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 diff --git a/src/fodesystem/newton_gregory.jl b/src/fodesystem/newton_gregory.jl index a68f1a33..8a679092 100644 --- a/src/fodesystem/newton_gregory.jl +++ b/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 diff --git a/src/fodesystem/newton_polynomials.jl b/src/fodesystem/newton_polynomials.jl index 02819732..f781546e 100644 --- a/src/fodesystem/newton_polynomials.jl +++ b/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 \ No newline at end of file diff --git a/src/fodesystem/trapezoid.jl b/src/fodesystem/trapezoid.jl index 50bfc179..6bf75072 100644 --- a/src/fodesystem/trapezoid.jl +++ b/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 diff --git a/src/lyapunov.jl b/src/lyapunov.jl index 6f750242..469f0c97 100644 --- a/src/lyapunov.jl +++ b/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 diff --git a/test/FDiscreteTests.jl b/test/discrete.jl similarity index 100% rename from test/FDiscreteTests.jl rename to test/discrete.jl diff --git a/test/DODETests.jl b/test/dode.jl similarity index 100% rename from test/DODETests.jl rename to test/dode.jl diff --git a/test/FDDETests.jl b/test/fdde.jl similarity index 100% rename from test/FDDETests.jl rename to test/fdde.jl diff --git a/test/FFODETests.jl b/test/ffode.jl similarity index 100% rename from test/FFODETests.jl rename to test/ffode.jl diff --git a/test/FODETests.jl b/test/fode.jl similarity index 87% rename from test/FODETests.jl rename to test/fode.jl index 0205b5c9..5650ce04 100644 --- a/test/FODETests.jl +++ b/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 diff --git a/test/FOLEtest.jl b/test/fole.jl similarity index 100% rename from test/FOLEtest.jl rename to test/fole.jl diff --git a/test/runtests.jl b/test/runtests.jl index b96d2df7..48a96e19 100644 --- a/test/runtests.jl +++ b/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