From 8e5db53b8ace2d3c2bd66ee4c1e53ac4601f2e4f Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Fri, 28 Jul 2023 19:01:48 +0800 Subject: [PATCH] Add DiffEqBase and SciMLBase Signed-off-by: ErikQQY <2283984853@qq.com> --- Project.toml | 3 + src/FractionalDiffEq.jl | 21 ++---- src/ffode/atangana_seda.jl | 2 + src/fodesystem/Euler.jl | 12 ++-- src/fodesystem/GLWithMemory.jl | 25 +++---- src/fodesystem/NonLinear.jl | 16 ++--- src/fodesystem/PIPECE.jl | 29 ++++---- src/fodesystem/atangana_seda.jl | 26 ++++---- src/fodesystem/bdf.jl | 10 +-- src/fodesystem/explicit_pi.jl | 29 ++++---- src/fodesystem/newton_gregory.jl | 10 +-- src/fodesystem/newton_polynomials.jl | 12 ++-- src/fodesystem/trapezoid.jl | 10 +-- src/lyapunov.jl | 10 +-- src/multitermsfode/explicit_pi.jl | 1 + src/singletermfode/Euler.jl | 32 --------- src/singletermfode/GL.jl | 43 ------------ src/singletermfode/PECE.jl | 92 -------------------------- src/singletermfode/atangana_seda.jl | 33 --------- src/singletermfode/product_integral.jl | 47 ------------- src/types/problems.jl | 69 +++++++++++++------ src/utils.jl | 18 ----- test/FODETests.jl | 42 ++++++------ test/auxillary.jl | 11 +-- 24 files changed, 183 insertions(+), 420 deletions(-) delete mode 100644 src/singletermfode/Euler.jl delete mode 100644 src/singletermfode/GL.jl delete mode 100644 src/singletermfode/PECE.jl delete mode 100644 src/singletermfode/atangana_seda.jl delete mode 100644 src/singletermfode/product_integral.jl diff --git a/Project.toml b/Project.toml index 394adb56c..dc83a7bf4 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Qingyu Qu "] version = "0.3.1" [deps] +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a" @@ -12,10 +13,12 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" SpecialMatrices = "928aab9d-ef52-54ac-8ca1-acd7ca42c160" ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24" +TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] diff --git a/src/FractionalDiffEq.jl b/src/FractionalDiffEq.jl index 31dfe67e8..b994a6280 100644 --- a/src/FractionalDiffEq.jl +++ b/src/FractionalDiffEq.jl @@ -1,6 +1,7 @@ module FractionalDiffEq using LinearAlgebra +using SciMLBase, DiffEqBase using SpecialFunctions using SparseArrays using InvertedIndices: Not @@ -13,19 +14,12 @@ using ToeplitzMatrices using RecipesBase using ForwardDiff using Polynomials: Polynomial +using TruncatedStacktraces include("types/problems.jl") include("types/algorithms.jl") include("types/solutions.jl") - -# Single-term fractional ordinary differential equations -include("singletermfode/PECE.jl") -include("singletermfode/product_integral.jl") -include("singletermfode/GL.jl") -include("singletermfode/atangana_seda.jl") -include("singletermfode/Euler.jl") - # Multi-terms fractional ordinary differential equations include("multitermsfode/matrix.jl") include("multitermsfode/hankelmatrix.jl") @@ -80,7 +74,7 @@ export solve, FDEProblem export FDDEProblem, FDDESystem, FDDEMatrixProblem # FODE problems -export SingleTermFODEProblem, MultiTermsFODEProblem, FODESystem, DODEProblem, FFPODEProblem, FFEODEProblem, FFMODEProblem +export FODEProblem, MultiTermsFODEProblem, DODEProblem, FFPODEProblem, FFEODEProblem, FFMODEProblem # Fractional Discrete probelms export FractionalDiscreteProblem, FractionalDiscreteSystem @@ -94,14 +88,15 @@ export FODESolution, FDifferenceSolution, DODESolution, FFMODESolution export FODESystemSolution, FDDESystemSolution, FFMODESystem # FODE solvers -export PIEX, PIPECE, PIRect, PITrap +export PIPECE, PIRect, PITrap export PECE, FODEMatrixDiscrete, ClosedForm, ClosedFormHankelM, GL -export AtanganaSeda, AtanganaSedaAB -export Euler +export AtanganaSedaAB +#export Euler # System of FODE solvers export NonLinearAlg, FLMMBDF, FLMMNewtonGregory, FLMMTrap, PIEX, NewtonPolynomial export AtanganaSedaCF +export AtanganaSeda # FDDE solvers export DelayPECE, DelayABM, MatrixForm @@ -129,6 +124,4 @@ export FOLyapunov, FOLE # Distributed order auxiliary SpecialFunctions export DOB, DOF, DORANORT, isFunction -export ourfft, ourifft, rowfft, FastConv - end \ No newline at end of file diff --git a/src/ffode/atangana_seda.jl b/src/ffode/atangana_seda.jl index f93bc06f5..e6f0c9a6e 100644 --- a/src/ffode/atangana_seda.jl +++ b/src/ffode/atangana_seda.jl @@ -1,3 +1,5 @@ +struct AtanganaSeda <: FODESystemAlgorithm end + function solve(prob::Union{FFMODEProblem, FFMODESystem}, h::Float64, ::AtanganaSeda) if typeof(prob.order[2]) <: Function solve_cf_variable_ffodesystem(prob, h) diff --git a/src/fodesystem/Euler.jl b/src/fodesystem/Euler.jl index a9f214e90..f54846412 100644 --- a/src/fodesystem/Euler.jl +++ b/src/fodesystem/Euler.jl @@ -1,7 +1,9 @@ -function solve(prob::FODESystem, h, ::Euler) - @unpack f, α, u0, tspan, p = prob +struct Euler <: FODESystemAlgorithm end + +function solve(prob::FODEProblem, h, ::Euler) + @unpack f, order, u0, tspan, p = prob t0=tspan[1]; tfinal=tspan[2] - α = α[1] + order = order[1] t = collect(Float64, t0:h:tfinal) N::Int = ceil(Int, (tfinal-t0)/h) l = length(u0) @@ -12,9 +14,9 @@ function solve(prob::FODESystem, h, ::Euler) tmp = zeros(Float64, l) for j=1:n f(tmp, result[:, j], p, t[j]) - temp = ((n-j+1)^α-(n-j)^α)*tmp + temp = ((n-j+1)^order-(n-j)^order)*tmp end - result[:, n+1] = result[:, 1] + h^α/gamma(α+1)*temp + result[:, n+1] = result[:, 1] + h^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 8494bccc8..9c0a8e2ed 100644 --- a/src/fodesystem/GLWithMemory.jl +++ b/src/fodesystem/GLWithMemory.jl @@ -1,7 +1,7 @@ """ # Usage - solve(prob::FODESystem, h, GL()) + solve(prob::FODEProblem, h, GL()) Use Grunwald Letnikov difference method to solve system of system of FODE. @@ -18,16 +18,17 @@ doi={10.1109/MOCAST.2019.8742063}} Python version by https://github.com/DClementeL/Grunwald_Letnikov """ -# Grunwald Letnikov discretization method dispatch for FODESystem +# Grunwald Letnikov discretization method dispatch for FODEProblem # struct GLWithMemory <: FractionalDiffEqAlgorithm end +struct GL <: FODESystemAlgorithm end -function solve(prob::FODESystem, h, ::GL) +function solve(prob::FODEProblem, h, ::GL) # GL method is only for same order FODE - @unpack f, α, u0, tspan, p = prob + @unpack f, order, u0, tspan, p = prob t0 = tspan[1]; T = tspan[2] t = collect(Float64, t0:h:T) - α = α[1] - hα = h^α[1] + order = order[1] + horder = h^order[1] n::Int64 = floor(Int64, (T-t0)/h)+1 l = length(u0) @@ -35,12 +36,12 @@ function solve(prob::FODESystem, h, ::GL) result = zeros(Float64, length(u0), n) result[:, 1] = u0 - # generating generalized binomial Cα - Cα = zeros(Float64, n) - Cα[1] = 1 + # generating generalized binomial Corder + Corder = zeros(Float64, n) + Corder[1] = 1 @fastmath @inbounds @simd for j in range(2, n, step=1) - Cα[j] = (1-(1+α)/(j-1))*Cα[j-1] + Corder[j] = (1-(1+order)/(j-1))*Corder[j-1] end du = zeros(Float64, l) @@ -50,12 +51,12 @@ function solve(prob::FODESystem, h, ::GL) @fastmath @inbounds @simd for j in range(1, k-1, step=1) for i in eachindex(summation) - summation[i] += Cα[j+1]*result[i, k-j] + summation[i] += Corder[j+1]*result[i, k-j] end end f(du, result[:, k-1], p, t0+(k-1)*h) - result[:, k] = @. hα*du-summation + result[:, k] = @. horder*du-summation end return FODESystemSolution(t, result) end \ No newline at end of file diff --git a/src/fodesystem/NonLinear.jl b/src/fodesystem/NonLinear.jl index 38e1bb4f2..2c23c095e 100644 --- a/src/fodesystem/NonLinear.jl +++ b/src/fodesystem/NonLinear.jl @@ -1,7 +1,7 @@ """ # Usage - solve(prob::FODESystem, h, NonLinearAlg()) + solve(prob::FODEProblem, h, NonLinearAlg()) Nonlinear algorithm for nonlinear fractional differential equations. @@ -11,8 +11,8 @@ Dingyu Xue, Northeastern University, China ISBN:9787030543981 """ struct NonLinearAlg <: FODESystemAlgorithm end -function solve(prob::FODESystem, h, ::NonLinearAlg, L0=1e10) - @unpack f, α, u0, tspan, p = prob +function solve(prob::FODEProblem, h, ::NonLinearAlg, L0=1e10) + @unpack f, order, u0, tspan, p = prob t0 = tspan[1]; T = tspan[2] time = collect(t0:h:T) n = length(u0) @@ -20,7 +20,7 @@ function solve(prob::FODESystem, h, ::NonLinearAlg, L0=1e10) g = genfun(1) g = g[:] u0 = u0[:] - ha = h.^α + ha = h.^order z = zeros(Float64, n, m) x1::Vector{Float64} = copy(u0) @@ -29,7 +29,7 @@ function solve(prob::FODESystem, h, ::NonLinearAlg, L0=1e10) W = zeros(n, SetMemoryEffect) #Initializing W a n*m matrix @fastmath @inbounds @simd for i = 1:n - W[i, :] = getvec(α[i], SetMemoryEffect, g) + W[i, :] = getvec(order[i], SetMemoryEffect, g) end du = zeros(n) @@ -60,12 +60,12 @@ function genfun(p) return (1 .-a')*inv(A') end -function getvec(α, n, g) +function getvec(order, n, g) p = length(g)-1 - b = 1 + α + b = 1 + order g0 = g[1] w = Float64[] - push!(w, g[1]^α) + push!(w, g[1]^order) @fastmath @inbounds @simd for m = 2:p M = m-1 diff --git a/src/fodesystem/PIPECE.jl b/src/fodesystem/PIPECE.jl index 63a6ac69c..794c80c53 100644 --- a/src/fodesystem/PIPECE.jl +++ b/src/fodesystem/PIPECE.jl @@ -1,6 +1,6 @@ #= """ - solve(prob::FODESystem, h, PECE()) + solve(prob::FODEProblem, h, PECE()) Use the Adams-Bashforth-Moulton method to solve the system of FODEs. @@ -29,27 +29,30 @@ mutable struct M an_fft bn_fft end + +struct PECE <: FODESystemAlgorithm end + #TODO: Rename as PIPECE -function solve(prob::FODESystem, h, ::PECE) - @unpack f, α, u0, tspan, p = prob +function solve(prob::FODEProblem, h, ::PECE) + @unpack f, order, u0, tspan, p = prob t0 = tspan[1]; T = tspan[2] METH = M(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)# Initialization mu_tol = 1.0e-6 mu = 1 - α = α[:] + order = order[:] # issue [#64](https://github.com/SciFracX/FractionalDiffEq.jl/issues/64) - max_order = findmax(α)[1] + max_order = findmax(order)[1] if max_order > 1 @error "This method doesn't support high order FDEs" end # Check compatibility size of the problem with number of fractional orders - alpha_length = length(α) + alpha_length = length(order) problem_size = size(u0, 1) - m_alpha = ceil.(Int, α) + m_alpha = ceil.(Int, order) m_alpha_factorial = zeros(alpha_length, maximum(m_alpha)) for i = 1 : alpha_length for j = 0 : m_alpha[i]-1 @@ -78,7 +81,7 @@ function solve(prob::FODESystem, h, ::PECE) bn = zeros(alpha_length, NNr+1); an = copy(bn); a0 = copy(bn) for i_alpha = 1:alpha_length find_alpha = Float64[] - if α[i_alpha] == α[1:i_alpha-1] + if order[i_alpha] == order[1:i_alpha-1] push!(find_alpha, i_alpha) end @@ -87,16 +90,16 @@ function solve(prob::FODESystem, h, ::PECE) an[i_alpha, :] = an[find_alpha[1], :] a0[i_alpha, :] = a0[find_alpha[1], :] else - nalpha = nvett.^α[i_alpha] + nalpha = nvett.^order[i_alpha] nalpha1 = nalpha.*nvett bn[i_alpha, :] = nalpha[2:end] - nalpha[1:end-1] an[i_alpha, :] = [1; (nalpha1[1:end-2] - 2*nalpha1[2:end-1] + nalpha1[3:end]) ] - a0[i_alpha, :] = [0; nalpha1[1:end-2]-nalpha[2:end-1].*(nvett[2:end-1].-α[i_alpha].-1)] + a0[i_alpha, :] = [0; nalpha1[1:end-2]-nalpha[2:end-1].*(nvett[2:end-1].-order[i_alpha].-1)] end end METH.bn = bn; METH.an = an; METH.a0 = a0 - METH.halpha1 = h.^α./gamma.(α.+1) - METH.halpha2 = h.^α./gamma.(α.+2) + METH.halpha1 = h.^order./gamma.(order.+1) + METH.halpha2 = h.^order./gamma.(order.+2) METH.mu = mu; METH.mu_tol = mu_tol # Evaluation of FFT of coefficients of the PECE method @@ -116,7 +119,7 @@ function solve(prob::FODESystem, h, ::PECE) coef_end = 2^l*r for i_alpha = 1 : alpha_length find_alpha = Float64[] - if α[i_alpha] == α[1:i_alpha-1] + if order[i_alpha] == order[1:i_alpha-1] push!(find_alpha, i_alpha) end if isempty(find_alpha) == false diff --git a/src/fodesystem/atangana_seda.jl b/src/fodesystem/atangana_seda.jl index 92005e4ff..fd565781b 100644 --- a/src/fodesystem/atangana_seda.jl +++ b/src/fodesystem/atangana_seda.jl @@ -5,13 +5,13 @@ Solve Atangana-Baleanu fractional order differential equations using Newton Poly """ struct AtanganaSedaAB <: FODESystemAlgorithm end -function solve(prob::FODESystem, h, ::AtanganaSedaAB) - @unpack f, α, u0, tspan, p = prob - α = α[1] +function solve(prob::FODEProblem, h, ::AtanganaSedaAB) + @unpack f, order, u0, tspan, p = prob + order = order[1] t0 = tspan[1]; tfinal = tspan[2] t = collect(t0:h:tfinal) N::Int = ceil(Int, (tfinal-t0)/h) - AB = 1-α+α/gamma(α) + AB = 1-order+order/gamma(order) t = collect(Float64, t0:h:tfinal) l = length(u0) @@ -36,16 +36,16 @@ function solve(prob::FODESystem, h, ::AtanganaSedaAB) f(temp1, result[:, n-1], p, t[n-1]) for j=3:n f(temp1, result[:, j-2], p, t[j-2]) - temptemp1 += ((n+1-j)^α-(n-j)^α)*temp1 + temptemp1 += ((n+1-j)^order-(n-j)^order)*temp1 f(temp2, result[:, j-1], p, t[j-1]) - temptemp2 += (temp2-temp1)*((n+1-j)^α*(n-j+3+2*α)-(n-j)^α*(n-j+3+3*α)) + temptemp2 += (temp2-temp1)*((n+1-j)^order*(n-j+3+2*order)-(n-j)^order*(n-j+3+3*order)) f(temp3, result[:, j], p, t[j]) - temptemp3 += (temp3-2*temp2+temp1)*((n-j+1)^α*(2*(n-j)^2+(3*α+10)*(n-j)+2*(α)^2+9*α+12)-(n-j)^α*(2*(n-j)^2+(5*α+10)*(n-j)+6*α^2+18*α+12)) + temptemp3 += (temp3-2*temp2+temp1)*((n-j+1)^order*(2*(n-j)^2+(3*order+10)*(n-j)+2*(order)^2+9*order+12)-(n-j)^order*(2*(n-j)^2+(5*order+10)*(n-j)+6*order^2+18*order+12)) end temptemptemp = zeros(l) f(temptemptemp, result[:, n], p, t[n]) - result[:, n+1] = u0+(1-α)/AB*temptemptemp+((h^α)*α/(AB*gamma(α+1)))*temptemp1 + ((h^α)*α/(AB*gamma(α+2)))*temptemp2+((h^α)*α/(2*AB*gamma(α+3)))*temptemp3 + result[:, n+1] = u0+(1-order)/AB*temptemptemp+((h^order)*order/(AB*gamma(order+1)))*temptemp1 + ((h^order)*order/(AB*gamma(order+2)))*temptemp2+((h^order)*order/(2*AB*gamma(order+3)))*temptemp3 end return FODESystemSolution(t, result) end @@ -65,12 +65,12 @@ https://doi.org/10.1016/c2020-0-02711-8 """ struct AtanganaSedaCF <: FODESystemAlgorithm end #FIXME: Tests -function solve(prob::FODESystem, h, ::AtanganaSedaCF) - @unpack f, α, u0, tspan, p = prob +function solve(prob::FODEProblem, h, ::AtanganaSedaCF) + @unpack f, order, u0, tspan, p = prob t0 = tspan[1]; tfinal = tspan[2] t=collect(Float64, t0:h:tfinal) - α=α[1] - M=1-α+α/gamma(α) + order=order[1] + M=1-order+order/gamma(order) N=ceil(Int, (tfinal-t0)/h) l=length(u0) result = zeros(l, N+1) @@ -85,7 +85,7 @@ function solve(prob::FODESystem, h, ::AtanganaSedaCF) f(tempn1, result[:, n-1], p, t[n-1]) f(tempn2, result[:, n-2], p, t[n-2]) - result[:, n+1] = result[:, n] + (1-α)/M*(tempn-tempn1)+α*M*h*(23/12*tempn-4/3*tempn1+5/12*tempn2) + result[:, n+1] = result[:, n] + (1-order)/M*(tempn-tempn1)+order*M*h*(23/12*tempn-4/3*tempn1+5/12*tempn2) end return FODESystemSolution(t, result) end \ No newline at end of file diff --git a/src/fodesystem/bdf.jl b/src/fodesystem/bdf.jl index 4396ae250..9eedbfd1e 100644 --- a/src/fodesystem/bdf.jl +++ b/src/fodesystem/bdf.jl @@ -1,5 +1,5 @@ """ - solve(prob::FODESystem, h, FLMMBDF()) + solve(prob::FODEProblem, h, FLMMBDF()) Use [Backward Differentiation Formula](https://en.wikipedia.org/wiki/Backward_differentiation_formula) generated weights fractional linear multi-steps method to solve system of FODE. @@ -17,13 +17,13 @@ Use [Backward Differentiation Formula](https://en.wikipedia.org/wiki/Backward_di """ struct FLMMBDF <: FODESystemAlgorithm end -function solve(prob::FODESystem, h, ::FLMMBDF; reltol=1e-6, abstol=1e-6, maxiters = 100) - @unpack f, α, u0, tspan, p = prob +function solve(prob::FODEProblem, h, ::FLMMBDF; reltol=1e-6, abstol=1e-6, maxiters = 100) + @unpack f, order, u0, tspan, p = prob t0 = tspan[1]; tfinal = tspan[2] - alpha = α[1] + alpha = order[1] # issue [#64](https://github.com/SciFracX/FractionalDiffEq.jl/issues/64) - max_order = findmax(α)[1] + max_order = findmax(order)[1] if max_order > 1 @error "This method doesn't support high order FDEs" end diff --git a/src/fodesystem/explicit_pi.jl b/src/fodesystem/explicit_pi.jl index 2b7b607a8..10a80c492 100644 --- a/src/fodesystem/explicit_pi.jl +++ b/src/fodesystem/explicit_pi.jl @@ -1,6 +1,7 @@ + """ - solve(prob::FODESystem, h, PIEX()) + solve(prob::FODEProblem, h, PIEX()) Use explicit Product integration method to solve system of FODE. """ @@ -19,24 +20,24 @@ mutable struct M bn_fft end -function solve(prob::FODESystem, h, ::PIEX) - @unpack f, α, u0, tspan, p = prob +function solve(prob::FODEProblem, h, ::PIEX) + @unpack f, order, u0, tspan, p = prob t0 = tspan[1]; T = tspan[2] # issue [#64](https://github.com/SciFracX/FractionalDiffEq.jl/issues/64) - max_order = findmax(α)[1] + max_order = findmax(order)[1] if max_order > 1 @error "This method doesn't support high order FDEs" end METH = M(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)# Initialization - α = α[:] + order = order[:] # Check compatibility size of the problem with number of fractional orders - alpha_length = length(α) + alpha_length = length(order) problem_size = size(u0, 1) - m_alpha = ceil.(Int, α) + m_alpha = ceil.(Int, order) m_alpha_factorial = zeros(alpha_length, maximum(m_alpha)) for i = 1 : alpha_length for j = 0 : m_alpha[i]-1 @@ -63,21 +64,21 @@ function solve(prob::FODESystem, h, ::PIEX) bn = zeros(alpha_length, NNr+1)#; an = copy(bn); a0 = copy(bn) for i_alpha = 1:alpha_length find_alpha = Float64[] - if α[i_alpha] == α[1:i_alpha-1] + if order[i_alpha] == order[1:i_alpha-1] push!(find_alpha, i_alpha) end if isempty(find_alpha) == false bn[i_alpha, :] = bn[find_alpha[1], :] - elseif abs(α[i_alpha]-1) < 1e-14 + elseif abs(order[i_alpha]-1) < 1e-14 bn[i_alpha, :] = [1; zeros(NNr)] else - nalpha = nvett.^α[i_alpha] + nalpha = nvett.^order[i_alpha] bn[i_alpha, :] = nalpha[2:end] - nalpha[1:end-1] end end METH.bn = bn - METH.halpha1 = h.^α./gamma.(α.+1) + METH.halpha1 = h.^order./gamma.(order.+1) # Evaluation of FFT of coefficients of the PECE method METH.r = r @@ -96,7 +97,7 @@ function solve(prob::FODESystem, h, ::PIEX) coef_end = 2^l*r for i_alpha = 1 : alpha_length find_alpha = Float64[] - if α[i_alpha] == α[1:i_alpha-1] + if order[i_alpha] == order[1:i_alpha-1] push!(find_alpha, i_alpha) end if isempty(find_alpha) == false @@ -113,14 +114,14 @@ function solve(prob::FODESystem, h, ::PIEX) t = t0 .+ collect(0:N)*h y[:, 1] = u0[:, 1] fy[:, 1] = f_temp - (y, fy) = PIEX_system_triangolo(1, r-1, t, y, fy, zn, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f, α, p) + (y, fy) = PIEX_system_triangolo(1, r-1, t, y, fy, zn, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f, order, p) # Main process of computation by means of the FFT algorithm ff = zeros(1, 2^(Qr+2)); ff[1:2] = [0; 2] ; card_ff = 2 nx0::Int = 0; ny0::Int = 0 for qr = 0 : Qr L = 2^qr - (y, fy) = PIEX_system_disegna_blocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, zn, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f, α, p) + (y, fy) = PIEX_system_disegna_blocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, zn, N, METH, problem_size, alpha_length, m_alpha, m_alpha_factorial, u0, t0, f, order, p) ff[1:2*card_ff] = [ff[1:card_ff]; ff[1:card_ff]] card_ff = 2*card_ff ff[card_ff] = 4*L diff --git a/src/fodesystem/newton_gregory.jl b/src/fodesystem/newton_gregory.jl index 6341e3031..effc036d4 100644 --- a/src/fodesystem/newton_gregory.jl +++ b/src/fodesystem/newton_gregory.jl @@ -1,5 +1,5 @@ """ - solve(prob::FODESystem, h, FLMMNewtonGregory()) + solve(prob::FODEProblem, h, FLMMNewtonGregory()) Use [Newton Gregory](https://www.geeksforgeeks.org/newton-forward-backward-interpolation/) generated weights fractional linear multiple steps method to solve system of FODE. @@ -17,13 +17,13 @@ Use [Newton Gregory](https://www.geeksforgeeks.org/newton-forward-backward-inter """ struct FLMMNewtonGregory <: FODESystemAlgorithm end -function solve(prob::FODESystem, h, ::FLMMNewtonGregory; reltol=1e-6, abstol=1e-6, maxiters = 100) - @unpack f, α, u0, tspan, p = prob +function solve(prob::FODEProblem, h, ::FLMMNewtonGregory; reltol=1e-6, abstol=1e-6, maxiters = 100) + @unpack f, order, u0, tspan, p = prob t0 = tspan[1]; tfinal = tspan[2] - alpha = α[1] + alpha = order[1] # issue [#64](https://github.com/SciFracX/FractionalDiffEq.jl/issues/64) - max_order = findmax(α)[1] + max_order = findmax(order)[1] if max_order > 1 @error "This method doesn't support high order FDEs" end diff --git a/src/fodesystem/newton_polynomials.jl b/src/fodesystem/newton_polynomials.jl index 0c770ff8d..e658f5822 100644 --- a/src/fodesystem/newton_polynomials.jl +++ b/src/fodesystem/newton_polynomials.jl @@ -1,5 +1,5 @@ """ - solve(prob::FODESystem, h, NewtonPolynomial()) + solve(prob::FODEProblem, h, NewtonPolynomial()) Solve FODE system using Newton Polynomials methods. @@ -12,12 +12,12 @@ https://doi.org/10.1016/c2020-0-02711-8 """ struct NewtonPolynomial <: FODESystemAlgorithm end -function solve(prob::FODESystem, h, ::NewtonPolynomial) - @unpack f, α, u0, tspan, p = prob +function solve(prob::FODEProblem, h, ::NewtonPolynomial) + @unpack f, order, u0, tspan, p = prob t0 = tspan[1]; tfinal = tspan[2] - α = α[1] + order = order[1] t = collect(Float64, t0:h:tfinal) - M = 1-α+α/gamma(α) + M = 1-order+order/gamma(order) N::Int = ceil(Int, (tfinal-t0)/h) l = length(u0) result = zeros(Float64, l, N+1) @@ -38,7 +38,7 @@ function solve(prob::FODESystem, h, ::NewtonPolynomial) f(temp1, result[:, n], p, t[n]) f(temp2, result[:, n-1], p, t[n-1]) f(temp3, result[:, n-2], p, t[n-2]) - result[:, n+1] = result[:, n] + (1-α)/M*(temp1-temp2)+α.*M.*h.*(23/12*temp1-4/3*temp2+5/12*temp3) + result[:, n+1] = result[:, n] + (1-order)/M*(temp1-temp2)+order.*M.*h.*(23/12*temp1-4/3*temp2+5/12*temp3) end return FODESystemSolution(t, result) end \ No newline at end of file diff --git a/src/fodesystem/trapezoid.jl b/src/fodesystem/trapezoid.jl index aea809e2d..a45f6e594 100644 --- a/src/fodesystem/trapezoid.jl +++ b/src/fodesystem/trapezoid.jl @@ -1,5 +1,5 @@ """ - solve(prob::FODEsystem, FLMMTrap()) + solve(prob::FODEProblem, FLMMTrap()) Use [Trapezoidal](https://en.wikipedia.org/wiki/Trapezoidal_rule_(differential_equations)) with generating function ``f(x)=\\frac{1+x}{2(1-x)^\\alpha}`` generated weights fractional linear multiple steps method to solve system of FODE. @@ -17,15 +17,15 @@ Use [Trapezoidal](https://en.wikipedia.org/wiki/Trapezoidal_rule_(differential_e """ struct FLMMTrap <: FODESystemAlgorithm end -function solve(prob::FODESystem, h, ::FLMMTrap; reltol=1e-6, abstol=1e-6) - @unpack f, α, u0, tspan, p = prob +function solve(prob::FODEProblem, h, ::FLMMTrap; reltol=1e-6, abstol=1e-6) + @unpack f, order, u0, tspan, p = prob t0 = tspan[1]; tfinal = tspan[2] u0 = u0 - alpha = α[1] + alpha = order[1] maxiters = 100 # issue [#64](https://github.com/SciFracX/FractionalDiffEq.jl/issues/64) - max_order = findmax(α)[1] + max_order = findmax(order)[1] if max_order > 1 @error "This method doesn't support high order FDEs" end diff --git a/src/lyapunov.jl b/src/lyapunov.jl index 721c62063..163fd98b2 100644 --- a/src/lyapunov.jl +++ b/src/lyapunov.jl @@ -42,9 +42,9 @@ end # Dispatch for parameters not specified FOLyapunov(fun, order, t_start, h_norm, t_end, u0, h, out) = FOLyapunov(fun, order, t_start, h_norm, t_end, u0, h, out, nothing) -#FOLyapunov(sys::FODESystem, h_norm, h, out) = FOLyapunov(sys.f, sys.α, sys.tspan[1], h_norm, sys.tspan[2], sys.u0, h, out) -FOLyapunov(sys::FODESystem, h_norm, h, out, p) = FOLyapunov(sys.f, sys.α, sys.tspan[1], h_norm, sys.tspan[2], sys.u0, h, out, p) -FOLyapunov(sys::FODESystem, h_norm, h, out) = FOLyapunov(sys.f, sys.α, sys.tspan[1], h_norm, sys.tspan[2], sys.u0, h, out, nothing) +#FOLyapunov(sys::FODEProblem, h_norm, h, out) = FOLyapunov(sys.f, sys.α, sys.tspan[1], h_norm, sys.tspan[2], sys.u0, h, out) +FOLyapunov(sys::FODEProblem, h_norm, h, out, p) = FOLyapunov(sys.f, sys.α, sys.tspan[1], h_norm, sys.tspan[2], sys.u0, h, out, p) +FOLyapunov(sys::FODEProblem, h_norm, h, out) = FOLyapunov(sys.f, sys.α, sys.tspan[1], h_norm, sys.tspan[2], sys.u0, h, out, nothing) @inline function is_all_equal(order::AbstractArray) length(order) < 2 && return true @@ -85,7 +85,7 @@ function noncommensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out LExp = zeros(ne) for it=1:n_it # Here we directly use the buildin PECE algorithm to solve the extend system, SO FAST!!! - prob = FODESystem(extend_fun, q[:], x[:], (t, t+h_norm), p) + prob = FODEProblem(extend_fun, q[:], x[:], (t, t+h_norm), p) sol = solve(prob, h, PECE()) Y = sol.u t = t+h_norm @@ -257,7 +257,7 @@ mutable struct M an_fft bn_fft end -#TODO: Decouple ABM methods for FODESystems from FractionalSystems.jl to FractionalDiffEq.jl +#TODO: Decouple ABM methods for FODEProblems from FractionalSystems.jl to FractionalDiffEq.jl function pc(alpha, f_fun, t0, T, y0, h, p) METH = M(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) mu_tol = 1.0e-6 diff --git a/src/multitermsfode/explicit_pi.jl b/src/multitermsfode/explicit_pi.jl index 51186284d..af60a54e5 100644 --- a/src/multitermsfode/explicit_pi.jl +++ b/src/multitermsfode/explicit_pi.jl @@ -1,6 +1,7 @@ #= Multiple dispatch of explicit product integral methods for multi term FODEProblem =# +struct PIEX <: FODESystemAlgorithm end function solve(prob::MultiTermsFODEProblem, h, ::PIEX) @unpack parameters, orders, rightfun, u0, tspan = prob diff --git a/src/singletermfode/Euler.jl b/src/singletermfode/Euler.jl deleted file mode 100644 index 2ccbce4fe..000000000 --- a/src/singletermfode/Euler.jl +++ /dev/null @@ -1,32 +0,0 @@ -""" - solve(prob, h, Euler()) - -The basic forward Euler method for fractional ordinary differential equations. - -```tex -@inproceedings{Li2015NumericalMF, - title={Numerical Methods for Fractional Calculus}, - author={Changpin Li and Fanhai Zeng}, - year={2015} -} -``` -""" -struct Euler <: SingleTermFODEAlgorithm end - -function solve(prob::SingleTermFODEProblem, h::Float64, ::Euler) - @unpack f, α, u0, tspan, p = prob - fun(t, u) = f(u, p, t) - t0 = tspan[1]; tfinal = tspan[2] - t = collect(Float64, t0:h:tfinal) - N::Int = ceil(Int, (tfinal-t0)/h) - y = zeros(Float64, N+1) - y[1] = sum(u0) - for n=1:N - temp = zero(Float64) - for j=1:n - temp += ((n-j+1)^α-(n-j)^α)*fun(t[j], y[j]) - end - y[n+1] = y[1] + h^α/gamma(α+1)*temp - end - return FODESolution(t, y) -end \ No newline at end of file diff --git a/src/singletermfode/GL.jl b/src/singletermfode/GL.jl deleted file mode 100644 index 9889b3dc7..000000000 --- a/src/singletermfode/GL.jl +++ /dev/null @@ -1,43 +0,0 @@ -""" -# Usage - - solve(prob::SingleTermFODEProblem, h, GL()) -Grunwald Letnikov method for fractional ordinary differential equations - -```tex -@INPROCEEDINGS{8742063, -author={Clemente-López, D. and Muñoz-Pacheco, J. M. and Félix-Beltrán, O. G. and Volos, C.}, -booktitle={2019 8th International Conference on Modern Circuits and Systems Technologies (MOCAST)}, -title={Efficient Computation of the Grünwald-Letnikov Method for ARM-Based Implementations of Fractional-Order Chaotic Systems}, -year={2019}, -doi={10.1109/MOCAST.2019.8742063}} -``` -""" -struct GL <: SingleTermFODEAlgorithm end - -function solve(FODE::SingleTermFODEProblem, h::Float64, ::GL) - @unpack f, α, u0, tspan, p = FODE - fun(t, u) = f(u, p, t) - t0 = tspan[1]; T = tspan[2] - N::Int = floor(Int, (T-t0)/h)+1 - c = zeros(Float64, N) - - cp::Float64 = 1.0 - for j = 1:N - c[j] = (1-(1+α)/j)*cp - cp = c[j] - end - - # Initialization - y = zeros(Float64, N) - y[1] = sum(u0) - - @fastmath @inbounds @simd for i = 2:N - right = 0 - @fastmath @inbounds @simd for j=1:i-1 - right += c[j]*y[i-j] - end - y[i] = fun(t0+(i-1)*h, y[i-1])*h^α - right - end - return FODESolution(collect(t0:h:T), y) -end \ No newline at end of file diff --git a/src/singletermfode/PECE.jl b/src/singletermfode/PECE.jl deleted file mode 100644 index af49cf19e..000000000 --- a/src/singletermfode/PECE.jl +++ /dev/null @@ -1,92 +0,0 @@ -""" -# Usage - - solve(prob::SingleTermFODEProblem, h, PECE()) - -Predict-Evaluate-Correct-Evaluate algorithm. - -For more details, please refer to [Predictor-Corrector algorithms](https://en.wikipedia.org/wiki/Predictor%E2%80%93corrector_method) - -This PECE algorithm is taken from Diethelm's paper. - -### References - -```tex -@article{ -title={A predictor-corrector approach for the numerical solution of fractional differential equations}, -author={Diethelm, Kai and Ford, Neville J. and Freed, Alan D.} -doi={https://doi.org/10.1023/A:1016592219341} -} -``` -""" -struct PECE <: SingleTermFODEAlgorithm end -#TODO: Use Richardson extrapolation to refine the PECE algorithms -#TODO: Rename as ABM - -""" - solve(problem::SingleTermFODEProblem, h, PECE()) - -Generalproblem solving API for solving FDE problems. -""" -function solve(FODE::SingleTermFODEProblem, h::Float64, ::PECE) - @unpack f, α, u0, tspan, p = FODE - fun(t, u) = f(u, p, t) - t0 = tspan[1]; T = tspan[2] - N::Int64 = round(Int, (T-t0)/h) - y = zeros(Float64, N+1) # Initialize the solution array - - - @fastmath @inbounds @simd for n ∈ 0:N - # Handling initial value - temp = zero(Float64) - for k in 0:ceil(Int, α)-1 - temp += u0[k+1]*(t0+(n+1)*h)^k/factorial(k) - end - y[n+1] = temp + h^α/gamma(α+2)*(fun(t0+(n+1)*h, predictor(fun, y, α, n, h, u0, t0)) + right(fun, y, α, n, h)) - end - - tspan = collect(Float64, t0:h:T) - return FODESolution(tspan, y) -end - -function right(fun::Function, y, α::Float64, n::Int64, h::Float64) - temp = zero(Float64) - - @fastmath @inbounds @simd for j = 0:n - temp += A(j, n, α)*fun(j*h, y[j+1]) - end - - return temp -end - -function predictor(fun::Function, y, α::Float64, n::Integer, h::Float64, u0::Union{Number, AbstractArray}, t0::Number) - predict = zero(Float64) - leftsum = zero(Float64) - - l::Int = ceil(Int, α) - - # Handling initial value - @fastmath @inbounds @simd for k in 0:l-1 - leftsum += u0[k+1]*(t0+(n+1)*h)^k/factorial(k) - end - - @fastmath @inbounds @simd for j in 0:n - predict += B(j, n, α)*fun(j*h, y[j+1]) - end - - return leftsum + h^α/α*predict -end - - -function A(j::Int64, n::Int64, α::Float64) - if j == 0 - return n^(α+1) - (n-α)*(n+1)^α - elseif 1 ≤ j ≤ n - return (n-j+2)^(α+1) + (n-j)^(α+1) - 2*(n-j+1)^(α+1) - elseif j == n+1 - return 1 - end -end - -# Generalized binomials -B(j::Int64, n::Int64, α::Float64) = ((n + 1 - j)^α - (n - j)^α) # Moved the h^α/α to the end of predictor: return leftsum + h^α/α*predict diff --git a/src/singletermfode/atangana_seda.jl b/src/singletermfode/atangana_seda.jl deleted file mode 100644 index 338ca6927..000000000 --- a/src/singletermfode/atangana_seda.jl +++ /dev/null @@ -1,33 +0,0 @@ -""" - solve(prob::SingleTermFODEProblem, h, AS()) - -Atangana-Seda method for Caputo single term FODE. -""" -struct AtanganaSeda <: SingleTermFODEAlgorithm end - -function solve(prob::SingleTermFODEProblem, h::Float64, ::AtanganaSeda) - @unpack f, α, u0, tspan, p = prob - fun(t, u) = f(u, p, t) - t0 = tspan[1]; tfinal = tspan[2] - N::Int = ceil(Int, (tfinal-t0)/h) - t = collect(Float64, t0:h:tfinal) - u = zeros(Float64, N+1) - - u[1] = u0 - u[2] = u[1]+h*fun(t[1], u[1]) # One step Euler method to first evaluate the second value - - u[3] = u[2]+(h/2)*(3*fun(t[2], u[2])-fun(t[1], u[1])) # Two-step Adams-Bashforth method to evaluate the third value - - for n=3:N - temp1=0 - temp2=0 - temp3=0 - for j=3:n - temp1 += fun(t[j-2], u[j-2])*((n-j+1)^α-(n-j)^α) - temp2 += (fun(t[j-1], u[j-1])-fun(t[j-2], u[j-2]))*((n-j+1)^α*(n.-j.+3 .+2*α).-(n.-j).^α.*(n.-j.+3 .+3*α)) - temp3 += ((fun(t[j], u[j])-2*fun(t[j-1], u[j-1])+fun(t[j-2], u[j-2])))*(((n-j+1)^α*(2*(n-j)^2+(3*α+10)*(n-j)+2*α^2+9*α+12)) - (n-j)^α*(2*(n-j)^2+(5*α+10)*(n-j)+6*α^2 +18*α+12)) - end - u[n+1] = @. u[1]+(h^α/gamma(α+1))*temp1+(h^(α)/gamma(α+2))*temp2 + (h^α/(2*gamma(α+3)))*temp3 - end - return FODESolution(t, u) -end \ No newline at end of file diff --git a/src/singletermfode/product_integral.jl b/src/singletermfode/product_integral.jl deleted file mode 100644 index 955d9a6b4..000000000 --- a/src/singletermfode/product_integral.jl +++ /dev/null @@ -1,47 +0,0 @@ -""" -# Usage - - solve(prob::SingleTermFODEProblem, h, PIEX()) - -### References - -```tex -@inproceedings{Garrappa2018NumericalSO, - title={Numerical Solution of Fractional Differential Equations: A Survey and a Software Tutorial}, - author={Roberto Garrappa}, - year={2018} -} -``` -""" -struct PIEX <: SingleTermFODEAlgorithm end - -function solve(FODE::SingleTermFODEProblem, h::Float64, ::PIEX) - @unpack f, α, u0, tspan, p = FODE - fun(t, u) = f(u, p, t) - t0 = tspan[1]; T = tspan[2] - N::Int64 = round(Int, (T-t0)/h) - y = zeros(Float64, N+1) - - y[1]=sum(u0) - for j in range(2, N+1, step=1) - middle=0 - @turbo for i=0:j-1 - middle += bcoefficients(j-i, α)*fun(t0+i*h, y[i+1]) - end - middle = middle/gamma(α+1) - y[j] = sum(u0) + middle*h^α - end - t = collect(Float64, t0:h:T) - return FODESolution(t, y) -end -#= -function acoefficients(n, α) - if n == 0 - return 1 - else - return ((n-1)^(α+1)-2*n^(α+1)+(n+1)^(α+1)) - end -end -=# - -bcoefficients(n, α) = ((n+1)^α-n^α) \ No newline at end of file diff --git a/src/types/problems.jl b/src/types/problems.jl index 5669cfc08..8724e77a2 100644 --- a/src/types/problems.jl +++ b/src/types/problems.jl @@ -1,18 +1,8 @@ -""" - FDEProblem - -General type for all kinds of problems in FractionalDiffEq.jl. -""" +using SciMLBase, DiffEqBase +abstract type AbstractFDEProblem <: SciMLBase.AbstractDEProblem end +abstract type AbstractFODEProblem{uType, tType, oType, isinplace} <: AbstractFDEProblem end +abstract type AbstractFDDEProblem{uType, tType, oType, lType, isinplace} <: AbstractFDEProblem end abstract type FDEProblem end - - -""" - MultiTermsFODEProblem(parameters, orders, rightfun, u0, T) - - MultiTermsFODEProblem(parameters, orders, rightfun, rparameters, rorders, u0, T) - -Define a multi-terms fractional ordinary differential equation. -""" struct MultiTermsFODEProblem <: FDEProblem parameters orders @@ -34,15 +24,52 @@ MultiTermsFODEProblem(parameters, orders, rightfun, u0, T) = MultiTermsFODEProbl Define a single term fractional ordinary differential equation, there are only one fractional differential operator in this problem. """ -struct SingleTermFODEProblem <: FDEProblem - f::Function - α::Float64 - u0::Union{AbstractArray, Number} - tspan::Union{Tuple, Number} - p::Union{AbstractArray, Number, Nothing} +#= +abstract type AbstractTestProblem{uType, tType, isinplace} <: SciMLBase.AbstractODEProblem{uType, tType, isinplace} end +=# + + +SciMLBase.isinplace(prob::AbstractFODEProblem{uType, tType, oType, iip}) where {uType, tType, oType, iip} = iip +SciMLBase.isinplace(prob::AbstractFDDEProblem{uType, tType, oType, lType, isinplace}) where {uType, tType, oType, lType, isinplace} = iip + +struct StandardFODEProblem end + +struct FODEProblem{uType, tType, oType, isinplace, P, F, bF, PT, K} <: + AbstractFODEProblem{uType, tType, oType, isinplace} + f::F + order::oType + u0::uType + tspan::tType + p::P + problem_type::PT + kwargs::K + SciMLBase.@add_kwonly function FODEProblem{iip}(f::SciMLBase.AbstractODEFunction, order, u0, tspan, + p = SciMLBase.NullParameters(), + problem_type = StandardFODEProblem(); + kwargs...) where {iip} + _tspan = SciMLBase.promote_tspan(tspan) + #warn_paramtype(p) + new{typeof(u0), typeof(_tspan), typeof(order), iip, typeof(p), + typeof(f), typeof(order), + typeof(problem_type), typeof(kwargs)}(f, order, u0, _tspan, p, + problem_type, kwargs) + end + + function FODEProblem{iip}(f, order, u0, tspan, p = SciMLBase.NullParameters(); kwargs...) where {iip} + FODEProblem(ODEFunction{iip}(f), order, u0, tspan, p; kwargs...) + end +end + +TruncatedStacktraces.@truncate_stacktrace SingleTermFODEProblem 3 1 2 + +function FODEProblem(f::SciMLBase.AbstractODEFunction, order, u0, tspan, args...; kwargs...) + FODEProblem{SciMLBase.isinplace(f, 4)}(f, order, u0, tspan, args...; kwargs...) +end + +function FODEProblem(f, order, u0, tspan, p = SciMLBase.NullParameters(); kwargs...) + FODEProblem(ODEFunction(f), order, u0, tspan, p; kwargs...) end -SingleTermFODEProblem(f::Function, α::Float64, u0::Union{AbstractArray, Number}, tspan::Union{Tuple, Number}) = SingleTermFODEProblem(f, α, u0, tspan, nothing) """ FDDEProblem(f, ϕ, α, τ, tspan) diff --git a/src/utils.jl b/src/utils.jl index dbdce3744..2ca5e1ea6 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -5,15 +5,6 @@ function Base.show(io::IO, sol::FODESolution) println(sol.u) end -function Base.show(io::IO, prob::SingleTermFODEProblem) - printstyled(typeof(prob), color=:light_blue) - printstyled(" with order ") - printstyled("$(prob.α)", color=:red) - println() - println("timespan: $(prob.tspan)") - println("u0: $(prob.u0)") -end - function Base.show(io::IO, prob::MultiTermsFODEProblem) printstyled(typeof(prob), color=:light_blue) printstyled(" with order ") @@ -23,15 +14,6 @@ function Base.show(io::IO, prob::MultiTermsFODEProblem) println("u0: $(prob.u0)") end -function Base.show(io::IO, prob::FODESystem) - printstyled(typeof(prob), color=:light_blue) - printstyled(" with order ") - printstyled("$(prob.α)", color=:red) - println() - println("timespan: $(prob.tspan)") - println("u0: $(prob.u0)") -end - function Base.show(io::IO, prob::FDDEProblem) printstyled(typeof(prob), color=:light_blue) printstyled(" with order ") diff --git a/test/FODETests.jl b/test/FODETests.jl index 87e2b8b9d..ccaa243db 100644 --- a/test/FODETests.jl +++ b/test/FODETests.jl @@ -1,3 +1,4 @@ +#= @testset "Test Diethelm PECE algorithms" begin fun(y, p, t) = 1-y prob = SingleTermFODEProblem(fun, 1.8, [0, 0], (0, 5)) @@ -23,7 +24,7 @@ end @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] @@ -83,7 +84,7 @@ end @test isapprox(result.u, [0.0; 0.08402140107687359; 0.3754974742112727]; atol=1e-3) end -@testset "Test GL method for FODESystem" begin +@testset "Test GL method for FODEProblem" begin h=0.5; tf=1 alpha = [0.99, 0.99, 0.99] x0 = [1, 0, 1] @@ -93,13 +94,13 @@ end du[2] = u[1]*(b-u[3])-u[2] du[3] = u[1]*u[2]-c*u[3] end - prob = FODESystem(testf!, alpha, x0, (0, tf)) + prob = FODEProblem(testf!, alpha, x0, (0, tf)) sol = solve(prob, h, GL()) @test isapprox(sol.u, [1.0 -4.04478 84.8074 0.0 13.5939 -51.1251 1.0 -0.352607 -27.5541]; atol=1e-4) end - +#= @testset "Test GL method" begin fun(y, p, t) = 1-y prob = SingleTermFODEProblem(fun, 0.5, 0, (0, 1)) @@ -117,7 +118,7 @@ end 0.5670394859746068 0.5808448788127054]; atol=1e-3) end - +=# @testset "Test Nonlinear method" begin function chua!(du, x, p, t) a = 10.725; b = 10.593 @@ -131,14 +132,14 @@ end α = [0.93, 0.99, 0.92]; x0 = [0.2; -0.1; 0.1]; h = 0.1; tspan = (0, 0.5); - prob = FODESystem(chua!, α, x0, tspan) + prob = FODEProblem(chua!, α, x0, tspan) result = solve(prob, h, NonLinearAlg()) @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 - +#= @testset "Test Product Integral Explicit method" begin fun(y, p, t)=1-y prob = SingleTermFODEProblem(fun, 0.5, 0, (0, 5)) @@ -156,6 +157,7 @@ end 0.8181895332990293 0.8276979913681598]; atol=1e-4) end +=# #= @testset "Test Product Integral Implicit method" begin fun(t, y)=1-y @@ -394,7 +396,7 @@ end 0.8406427306579508]; atol=1e-4) end -@testset "Test product integration PECE method for FODESystem" begin +@testset "Test product integration PECE method for FODEProblem" begin function test(du, u, p, t) du[1] = 1-4*u[1]+u[1]^2*u[2] du[2] = 3*u[1]-u[1]^2*u[2] @@ -404,14 +406,14 @@ end tspan=(0, 0.1); h=0.01 y0=[1.2; 2.8] param=[1 3] - prob = FODESystem(test, alpha, y0, tspan) + prob = FODEProblem(test, alpha, y0, tspan) sol = solve(prob, h, PECE()) @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 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 -@testset "Test explicit product integration method for FODESystem" begin +@testset "Test explicit product integration method for FODEProblem" begin function test(du, u, p, t) du[1] = 1-4*u[1]+u[1]^2*u[2] du[2] = 3*u[1]-u[1]^2*u[2] @@ -421,7 +423,7 @@ end tspan = (0, 0.2); h=0.015 y0=[1.2; 2.8] param=[1 3] - prob = FODESystem(test, alpha, y0, tspan) + prob = FODEProblem(test, alpha, y0, tspan) sol = solve(prob, h, PIEX()) @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 @@ -437,7 +439,7 @@ end alpha=[0.8; 0.8] y0=[0.2; 0.03] h=0.1; tspan=(0, 0.5) - prob = FODESystem(Brusselator, alpha, y0, tspan) + prob = FODEProblem(Brusselator, alpha, y0, tspan) sol = solve(prob, h, FLMMNewtonGregory()) @test isapprox(sol.u, [ 0.2 0.200531 0.201161 0.201809 0.202453 0.203089 @@ -453,7 +455,7 @@ end alpha=[0.8; 0.8] y0=[0.2; 0.03] h=0.1; tspan=(0, 0.5) - prob = FODESystem(Brusselator, alpha, y0, tspan) + prob = FODEProblem(Brusselator, alpha, y0, tspan) sol = solve(prob, h, FLMMBDF()) @test isapprox(sol.u, [ 0.2 0.200531 0.201161 0.201809 0.202453 0.203089 @@ -469,7 +471,7 @@ end alpha=[0.8; 0.8] y0=[0.2; 0.03] h=0.1; tspan=(0, 0.5) - prob = FODESystem(Brusselator, alpha, y0, tspan) + prob = FODEProblem(Brusselator, alpha, y0, tspan) sol = solve(prob, h, FLMMTrap()) @test isapprox(sol.u, [0.2 0.200531 0.201161 0.201808 0.202452 0.203088 @@ -487,7 +489,7 @@ end du[2] = -sin(u[3])-b*u[2] du[3] = -sin(u[1])-b*u[3] end - prob = FODESystem(fun, α, u0, (t0, tfinal)) + prob = FODEProblem(fun, α, u0, (t0, tfinal)) sol = solve(prob, h, NewtonPolynomial()) @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 @@ -495,7 +497,7 @@ end 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 - +#= @testset "Test Atangana Seda method" begin fun(y, p, t) = 1-y prob = SingleTermFODEProblem(fun, 0.5, 0, (0, 5)) @@ -513,8 +515,8 @@ end -18.488162171353544 44.23628311938973]; atol=1e-3) end - -@testset "Test Atangana Seda method for Atangana-Baleanu FodeSystem" begin +=# +@testset "Test Atangana Seda method for Atangana-Baleanu FODEProblem" begin t0=0;tfinal=2;h=0.5 α = [0.99, 0.99, 0.99] u0 = [-0.1; 0.1; -0.1] @@ -524,7 +526,7 @@ end du[2] = u[1]*u[2]-b*u[1]*u[3]-u[2]+δ du[3] = b*u[1]*u[2]+u[1]*u[3]-u[3] end - prob = FODESystem(fun, α, u0, (t0, tfinal)) + prob = FODEProblem(fun, α, u0, (t0, tfinal)) sol = solve(prob, h, AtanganaSedaAB()) @test isapprox(sol.u, [ -0.1 0.7 1.18511 -1.41522 -31.0872 @@ -542,7 +544,7 @@ end u0 = [0.5, 0.9, 0.1] tspan = (0, 2) α = ones(3)*0.98; - prob = FODESystem(LotkaVolterra, α, u0, tspan) + prob = FODEProblem(LotkaVolterra, α, u0, tspan) sol = solve(prob, 0.5, AtanganaSedaCF()) isapprox(sol.u, [ 0.0 1.375 3.37602 -5.46045 449.712 diff --git a/test/auxillary.jl b/test/auxillary.jl index 197f69833..b6d9092bf 100644 --- a/test/auxillary.jl +++ b/test/auxillary.jl @@ -51,13 +51,6 @@ end @test isFunction("Hello")==false end -@testset "Test singleterm show method" begin - # SingleTermFODEProblem - singlefun(x, y) = 1-y - singletermprob = SingleTermFODEProblem(singlefun, 1.8, 0, 5) - @test_nowarn show(singletermprob) -end - @testset "Test Multiterms show method" begin # MultiTermsFODEProblem rightfun(x, y) = 172/125*cos(4/5*x) @@ -68,7 +61,7 @@ end @test_nowarn show(multitermssol) end -@testset "Test FODESystem show method" begin +@testset "Test FODEProblem show method" begin # FODESystem h=0.5; tspan=(0, 1) alpha = [0.99, 0.99, 0.99] @@ -79,7 +72,7 @@ end du[2] = u[1]*(b-u[3])-u[2] du[3] = u[1]*u[2]-c*u[3] end - fodesystemprob = FODESystem(testf!, alpha, x0, tspan) + fodesystemprob = FODEProblem(testf!, alpha, x0, tspan) fodesystemsol = solve(fodesystemprob, h, GL()) @test_nowarn show(fodesystemprob)