From cbb523c1c4a807f31ce736a98da5f70e9bfddcf7 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Tue, 28 Nov 2023 22:25:10 +0800 Subject: [PATCH 01/14] Refactor FDDEProblem and solvers Signed-off-by: ErikQQY <2283984853@qq.com> --- Project.toml | 5 +- src/FractionalDiffEq.jl | 2 + src/delay/DelayABM.jl | 73 ++++---- src/delay/DelayPECE.jl | 185 ++++++++++++-------- src/delay/product_integral.jl | 34 ++-- src/fodesystem/explicit_pi.jl | 18 -- src/fodesystem/newton_gregory.jl | 6 - src/fodesystem/trapezoid.jl | 3 - src/lyapunov.jl | 21 --- src/multitermsfode/implicit_pi_trapezoid.jl | 2 - src/types/problem_utils.jl | 28 +++ src/types/problems.jl | 71 ++++++-- test/FDDETests.jl | 80 +++++---- 13 files changed, 299 insertions(+), 229 deletions(-) create mode 100644 src/types/problem_utils.jl diff --git a/Project.toml b/Project.toml index 07611ad0..6495a195 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,7 @@ LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" @@ -25,8 +26,8 @@ TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] -ApproxFun = "0.13" -ArrayInterface = "7" +ApproxFun = "0.13.25" +ArrayInterface = "7.6" DiffEqBase = "6" FFTW = "1.0.0" ForwardDiff = "0.10" diff --git a/src/FractionalDiffEq.jl b/src/FractionalDiffEq.jl index b994a628..ef6f8315 100644 --- a/src/FractionalDiffEq.jl +++ b/src/FractionalDiffEq.jl @@ -20,6 +20,8 @@ include("types/problems.jl") include("types/algorithms.jl") include("types/solutions.jl") +include("types/problem_utils.jl") + # Multi-terms fractional ordinary differential equations include("multitermsfode/matrix.jl") include("multitermsfode/hankelmatrix.jl") diff --git a/src/delay/DelayABM.jl b/src/delay/DelayABM.jl index ea449bdb..c50ed5a4 100644 --- a/src/delay/DelayABM.jl +++ b/src/delay/DelayABM.jl @@ -1,5 +1,5 @@ """ - solve(FDDE::FDDEProblem, h, DelayABM()) + solve(FDDE::FDDEProblem, dt, DelayABM()) Use the Adams-Bashforth-Moulton method to solve fractional delayed differential equations. @@ -16,11 +16,12 @@ Use the Adams-Bashforth-Moulton method to solve fractional delayed differential struct DelayABM <: FDDEAlgorithm end #FIXME: There are still some improvments about initial condition #FIXME: Fix DelayABM method for FDDESystem : https://www.researchgate.net/publication/245538900_A_Predictor-Corrector_Scheme_For_Solving_Nonlinear_Delay_Differential_Equations_Of_Fractional_Order -#FIXME: Also the problem definition f(t, ϕ, y) or f(t, y, ϕ)? -function solve(FDDE::FDDEProblem, h, ::DelayABM) - @unpack f, ϕ, α, τ, tspan = FDDE - N::Int = round(Int, tspan/h) - Ndelay::Int = round(Int, τ/h) +#FIXME: Also the problem definition f(t, h, y) or f(t, y, h)? +function solve(FDDE::FDDEProblem, dt, ::DelayABM) + @unpack f, h, order, constant_lags, tspan, p = FDDE + τ = constant_lags[1] + N::Int = round(Int, tspan[2]/dt) + Ndelay::Int = round(Int, τ/dt) x1 = zeros(Float64, Ndelay+N+1) x = zeros(Float64, Ndelay+N+1) #x1[Ndelay+N+1] = 0 @@ -28,32 +29,32 @@ function solve(FDDE::FDDEProblem, h, ::DelayABM) #x[Ndelay+N+1] = 0 # History function handling - #FIXME: When the value of history function ϕ is different with the initial value? - if typeof(ϕ) <: Number - x[1:Ndelay] = ϕ*ones(Ndelay) - elseif typeof(ϕ) <: Function - x[Ndelay] = ϕ(0) - x[1:Ndelay-1] .= ϕ(collect(Float64, -h*(Ndelay-1):h:(-h))) + #FIXME: When the value of history function h is different with the initial value? + if typeof(h) <: Number + x[1:Ndelay] = h*ones(Ndelay) + elseif typeof(h) <: Function + x[Ndelay] = h(p, 0) + x[1:Ndelay-1] .= h(p, collect(Float64, -dt*(Ndelay-1):dt:(-dt))) end x0 = copy(x[Ndelay]) - x1[Ndelay+1] = x0 + h^α*f(0, x[1], x0)/(gamma(α)*α) + x1[Ndelay+1] = x0 + dt^order*f(x[1], x0, p, 0)/(gamma(order)*order) - x[Ndelay+1] = x0 + h^α*(f(0, x[1], x[Ndelay+1]) + α*f(0, x[1], x0))/gamma(α+2) + x[Ndelay+1] = x0 + dt^order*(f(x[1], x[Ndelay+1], p, 0) + order*f(x[1], x0, p, 0))/gamma(order+2) @fastmath @inbounds @simd for n=1:N - M1=(n^(α+1)-(n-α)*(n+1)^α)*f(0, x[1], x0) + M1=(n^(order+1)-(n-order)*(n+1)^order)*f(x[1], x0, p, 0) - N1=((n+1)^α-n^α)*f(0, x[1], x0) + N1=((n+1)^order-n^order)*f(x[1], x0, p, 0) @fastmath @inbounds @simd for j=1:n - M1 = M1+((n-j+2)^(α+1)+(n-j)^(α+1)-2*(n-j+1)^(α+1))*f(0, x[j], x[Ndelay+j]) - N1 = N1+((n-j+1)^α-(n-j)^α)*f(0, x[j], x[Ndelay+j]) + M1 = M1+((n-j+2)^(order+1)+(n-j)^(order+1)-2*(n-j+1)^(order+1))*f(x[j], x[Ndelay+j], p, 0) + N1 = N1+((n-j+1)^order-(n-j)^order)*f(x[j], x[Ndelay+j], p, 0) end - x1[Ndelay+n+1] = x0+h^α*N1/(gamma(α)*α) - x[Ndelay+n+1] = x0+h^α*(f(0, x[n+1], x[Ndelay+n+1])+M1)/gamma(α+2) + x1[Ndelay+n+1] = x0+dt^order*N1/(gamma(order)*order) + x[Ndelay+n+1] = x0+dt^order*(f(x[n+1], x[Ndelay+n+1], p, 0)+M1)/gamma(order+2) end xresult = zeros(N-Ndelay+1) @@ -83,41 +84,41 @@ DelayABM method for system of fractional delay differential equations. ``` =# -function solve(FDDESys::FDDESystem, h, ::DelayABM) - @unpack f, ϕ, α, τ, T = FDDESys - len = length(ϕ) - N::Int = round(Int, T/h) - Ndelay = round(Int, τ/h) - t = collect(Float64, 0:h:(T-τ)) +function solve(FDDESys::FDDESystem, dt, ::DelayABM) + @unpack f, h, order, τ, T = FDDESys + len = length(h) + N::Int = round(Int, T/dt) + Ndelay = round(Int, τ/dt) + t = collect(Float64, 0:dt:(T-τ)) x1 = zeros(Ndelay+N+1, len) x = zeros(Ndelay+N+1, len) du = zeros(len) # Put the delay term in the array for i=1:len - x[1:Ndelay, i] = ϕ[i]*ones(Ndelay) + x[1:Ndelay, i] = h[i]*ones(Ndelay) end x0 = copy(x[Ndelay, :]) for i=1:len - αi = α[i] + αi = order[i] f(du, x0, x[1, :], 0) - x1[Ndelay+1, i] = x0[i] + h^αi*du[i]/(gamma(αi)*αi) + x1[Ndelay+1, i] = x0[i] + dt^αi*du[i]/(gamma(αi)*αi) end for i=1:len - αi = α[i] + αi = order[i] f(du, x[Ndelay+1, :], x[1, :], 0) du1 = copy(du) f(du, x0, x[1, :], 0) - x[Ndelay+1, i] = x0[i] + h^αi*(du1[i] + αi*du[i])/gamma(αi+2) + x[Ndelay+1, i] = x0[i] + dt^αi*(du1[i] + αi*du[i])/gamma(αi+2) end M1=zeros(len); N1=zeros(len) @fastmath @inbounds @simd for n=1:N for i=1:len - αi = α[i] + αi = order[i] f(du, x0, x[1, :], 0) M1[i]=(n^(αi+1)-(n-αi)*(n+1)^αi)*du[i] N1[i]=((n+1)^αi-n^αi)*du[i] @@ -125,7 +126,7 @@ function solve(FDDESys::FDDESystem, h, ::DelayABM) @fastmath @inbounds @simd for j=1:n for i=1:len - αi = α[i] + αi = order[i] f(du, x[Ndelay+j, :], x[j, :], 0) M1[i] = M1[i]+((n-j+2)^(αi+1)+(n-j)^(αi+1)-2*(n-j+1)^(αi+1))*du[i] N1[i] = N1[i]+((n-j+1)^αi-(n-j)^αi)*du[i] @@ -133,10 +134,10 @@ function solve(FDDESys::FDDESystem, h, ::DelayABM) end for i=1:len - αi = α[i] - x1[Ndelay+n+1, i] = x0[i]+h^αi*N1[i]/(gamma(αi)*αi) + αi = order[i] + x1[Ndelay+n+1, i] = x0[i]+dt^αi*N1[i]/(gamma(αi)*αi) f(du, x[Ndelay+n+1, :], x[n+1, :], 0) - x[Ndelay+n+1, i] = x0[i]+h^αi*(du[i]+M1[i])/gamma(αi+2) + x[Ndelay+n+1, i] = x0[i]+dt^αi*(du[i]+M1[i])/gamma(αi+2) end end diff --git a/src/delay/DelayPECE.jl b/src/delay/DelayPECE.jl index e8374c9b..6299337d 100644 --- a/src/delay/DelayPECE.jl +++ b/src/delay/DelayPECE.jl @@ -1,7 +1,7 @@ """ # Usage - solve(FDDE::FDDEProblem, h, DelayPECE()) + solve(FDDE::FDDEProblem, step_size, DelayPECE()) Using the delayed predictor-corrector method to solve the delayed fractional differential equation problem in the Caputo sense. @@ -34,87 +34,107 @@ Capable of solving both single term FDDE and multiple FDDE, support time varying """ struct DelayPECE <: FDDEAlgorithm end #FIXME: What if we have an FDDE with both variable order and time varying lag?? -function solve(FDDE::FDDEProblem, h, ::DelayPECE) +function solve(FDDE::FDDEProblem, step_size, ::DelayPECE) # If the delays are time varying, we need to specify single delay and multiple delay - if FDDE.τ isa Function + if FDDE.constant_lags[1] isa Function # Here is the PECE solver for single time varying lag - solve_fdde_with_single_lag(FDDE, h) - elseif FDDE.τ isa AbstractArray{Function} + solve_fdde_with_single_lag(FDDE, step_size) + elseif FDDE.constant_lags[1] isa AbstractArray{Function} # Here is the PECE solver for multiple time varying lags - solve_fdde_with_multiple_lags(FDDE, h) #TODO: implement this + solve_fdde_with_multiple_lags(FDDE, step_size) #TODO: implement this # Varying order fractional delay differential equations - elseif FDDE.α isa Function - if length(FDDE.τ) == 1 + elseif FDDE.order isa Function + if length(FDDE.cpnstant_lags[1]) == 1 # Here is the PECE solver for single lag with variable order - solve_fdde_with_single_lag_and_variable_order(FDDE, h) + solve_fdde_with_single_lag_and_variable_order(FDDE, step_size) else # Here is the PECE solver for multiple lags with variable order - solve_fdde_with_multiple_lags_and_variable_order(FDDE, h) + solve_fdde_with_multiple_lags_and_variable_order(FDDE, step_size) end else # If the delays are constant - if length(FDDE.τ) == 1 + if length(FDDE.constant_lags[1]) == 1 # Call the DelayPECE solver for single lag FDDE - solve_fdde_with_single_lag(FDDE, h) + solve_fdde_with_single_lag(FDDE, step_size) else # Call the DelayPECE solver for multiple lags FDDE - solve_fdde_with_multiple_lags(FDDE, h) + solve_fdde_with_multiple_lags(FDDE, step_size) end end end -function solve_fdde_with_single_lag(FDDE::FDDEProblem, h) - @unpack f, ϕ, α, τ, tspan = FDDE - T = tspan - t = collect(0:h:T) +function solve_fdde_with_single_lag(FDDE::FDDEProblem, step_size) + @unpack f, h, order, constant_lags, p, tspan = FDDE + τ = constant_lags[1] + iip = SciMLBase.isinplace(FDDE) + T = tspan[2] + t = collect(0:step_size:T) maxn = size(t, 1) - yp = collect(Float64, 0:h:T+h) + yp = collect(Float64, 0:step_size:T+step_size) y = copy(t) - y[1] = ϕ(0) + y[1] = h(p, 0) for n in 1:maxn-1 yp[n+1] = 0 for j = 1:n - yp[n+1] = yp[n+1]+generalized_binomials(j-1, n-1, α, h)*f(t[j], y[j], v(ϕ, j, τ, h, y, yp, t)) + if iip + tmp = zeros(length(yp[1])) + f(tmp, y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j]) + yp[n+1] = yp[n+1]+generalized_binomials(j-1, n-1, order, step_size)*tmp + else + yp[n+1] = yp[n+1]+generalized_binomials(j-1, n-1, order, step_size)*f(y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j]) + end end - yp[n+1] = yp[n+1]/gamma(α)+ϕ(0) + yp[n+1] = yp[n+1]/gamma(order)+h(p, 0) y[n+1] = 0 @fastmath @inbounds @simd for j=1:n - y[n+1] = y[n+1]+a(j-1, n-1, α, h)*f(t[j], y[j], v(ϕ, j, τ, h, y, yp, t)) + if iip + tmp = zeros(length(y[1])) + f(tmp, y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j]) + y[n+1] = y[n+1]+a(j-1, n-1, order, step_size)*tmp + else + y[n+1] = y[n+1]+a(j-1, n-1, order, step_size)*f(y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j]) + end end - y[n+1] = y[n+1]/gamma(α)+h^α*f(t[n+1], yp[n+1], v(ϕ, n+1, τ, h, y, yp, t))/gamma(α+2)+ϕ(0) + if iip + tmp = zeros(y[1]) + f(tmp, yp[n+1], v(h, n+1, τ, step_size, y, yp, t, p), p, t[n+1]) + y[n+1] = y[n+1]/gamma(order)+step_size^order*tmp/gamma(order+2)+h(p, 0) + else + y[n+1] = y[n+1]/gamma(order)+step_size^order*f(yp[n+1], v(h, n+1, τ, step_size, y, yp, t, p), p, t[n+1])/gamma(order+2)+h(p, 0) + end end V = copy(t) @fastmath @inbounds @simd for n = 1:maxn-1 # The delay term - V[n] = v(ϕ, n, τ, h, y, yp, t) + V[n] = v(h, n, τ, step_size, y, yp, t, p) end return V, y end -function a(j::Int, n::Int, α, h) +function a(j::Int, n::Int, order, step_size) if j == n+1 result = 1 elseif j == 0 - result = n^(α+1)-(n-α)*(n+1)^α + result = n^(order+1)-(n-order)*(n+1)^order else - result = (n-j+2)^(α+1) + (n-j)^(α+1) - 2*(n-j+1)^(α+1) + result = (n-j+2)^(order+1) + (n-j)^(order+1) - 2*(n-j+1)^(order+1) end - return result*h^α / (α*(α + 1)) + return result*step_size^order / (order*(order + 1)) end -generalized_binomials(j, n, α, h) = h^α/α*((n-j+1)^α - (n-j)^α) +generalized_binomials(j, n, order, step_size) = @. step_size^order/order*((n-j+1)^order - (n-j)^order) -function v(ϕ, n, τ, h, y, yp, t) +function v(h, n, τ, step_size, y, yp, t, p) if typeof(τ) <: Function - m = floor.(Int, τ.(t)/h) - δ = m.-τ.(t)./h - if τ(t[n]) >= (n-1)*h - return ϕ((n-1)*h-τ(t[n])) + m = floor.(Int, τ.(t)/step_size) + δ = m.-τ.(t)./step_size + if τ(t[n]) >= (n-1)*step_size + return h(p, (n-1)*step_size-τ(t[n])) else if m[n] > 1 return δ[n]*y[n-m[n]+2] + (1-δ[n])*y[n-m[n]+1] @@ -123,11 +143,11 @@ function v(ϕ, n, τ, h, y, yp, t) end end else - if τ >= (n-1)*h - return ϕ((n-1)*h-τ) + if τ >= (n-1)*step_size + return h(p, (n-1)*step_size-τ) else - m = floor(Int, τ/h) - δ = m-τ/h + m = floor(Int, τ/step_size) + δ = m-τ/step_size if m>1 return δ*y[n-m+2] + (1-δ)*y[n-m+1] @@ -139,33 +159,34 @@ function v(ϕ, n, τ, h, y, yp, t) end -function solve_fdde_with_multiple_lags(FDDE::FDDEProblem, h) - @unpack f, ϕ, α, τ, tspan = FDDE - t = collect(0:h:tspan) +function solve_fdde_with_multiple_lags(FDDE::FDDEProblem, step_size) + @unpack f, h, order, constant_lags, p, tspan = FDDE + τ = constant_lags[1] + t = collect(0:step_size:tspan[2]) maxn = length(t) yp = zeros(Float64, maxn) y = copy(t) - y[1] = ϕ(0) + y[1] = h(p, 0) for n in 1:maxn-1 yp[n+1] = 0 for j = 1:n - yp[n+1] = yp[n+1]+generalized_binomials(j-1, n-1, α, h)*f(t[j], y[j], multiv(ϕ, j, τ, h, y, yp)...) + yp[n+1] = yp[n+1]+generalized_binomials(j-1, n-1, order, step_size)*f(t[j], y[j], multiv(h, j, τ, step_size, y, yp, p)...) end - yp[n+1] = yp[n+1]/gamma(α)+ϕ(0) + yp[n+1] = yp[n+1]/gamma(order)+h(p, 0) y[n+1] = 0 for j=1:n - y[n+1] = y[n+1]+multia(j-1, n-1, α, h)*f(t[j], y[j], multiv(ϕ, j, τ, h, y, yp)...) + y[n+1] = y[n+1]+multia(j-1, n-1, order, step_size)*f(t[j], y[j], multiv(h, j, τ, step_size, y, yp, p)...) end - y[n+1] = y[n+1]/gamma(α)+h^α*f(t[n+1], yp[n+1], multiv(ϕ, n+1, τ, h, y, yp)...)/gamma(α+2) + ϕ(0) + y[n+1] = y[n+1]/gamma(order)+step_size^order*f(t[n+1], yp[n+1], multiv(h, n+1, τ, step_size, y, yp, p)...)/gamma(order+2) + h(p, 0) end V = [] for n = 1:maxn - push!(V, multiv(ϕ, n, τ, h, y, yp)) + push!(V, multiv(h, n, τ, step_size, y, yp, p)) end delayed = zeros(length(τ), length(V)) @@ -175,25 +196,25 @@ function solve_fdde_with_multiple_lags(FDDE::FDDEProblem, h) return delayed, y end -function multia(j::Int, n::Int, α, h) +function multia(j::Int, n::Int, order, step_size) if j == n+1 result = 1 elseif j == 0 - result = n^(α+1)-(n-α)*(n+1)^α + result = n^(order+1)-(n-order)*(n+1)^order elseif j == n - result = 2*(2^(α+1)-1) + result = 2*(2^(order+1)-1) else - result = (n-j+2)^(α+1) + (n-j)^(α+1) - 2*(n-j+1)^(α+1) + result = (n-j+2)^(order+1) + (n-j)^(order+1) - 2*(n-j+1)^(order+1) end - return result*h^α/(α*(α + 1)) + return result*step_size^order/(order*(order + 1)) end -function multiv(ϕ, n, τ, h, y, yp) - if maximum(τ) > n*h - return ϕ.((n-1)*h.-τ) +function multiv(h, n, τ, step_size, y, yp, p) + if maximum(τ) > n*step_size + return h.(p, (n-1)*step_size.-τ) else - m = floor.(Int, τ./h) - δ = m.-τ./h + m = floor.(Int, τ./step_size) + δ = m.-τ./step_size function judge(m) temp1 = findall(x->x>1, m) @@ -213,66 +234,76 @@ end #########################For variable order FDDE########################### -function solve_fdde_with_single_lag_and_variable_order(FDDE::FDDEProblem, h) - @unpack f, ϕ, α, τ, tspan = FDDE - T = tspan - t = collect(0:h:T) +function solve_fdde_with_single_lag_and_variable_order(FDDE::FDDEProblem, step_size) + @unpack f, order, h, constant_lags, p, tspan = FDDE + iip = SciMLBase.isinplace(FDDE) + order = order[1] + τ = constant_lags[1] + T = tspan[2] + t = collect(0:step_size:T) maxn = size(t, 1) - yp = collect(0:h:T+h) + yp = collect(0:step_size:T+step_size) y = copy(t) - y[1] = ϕ(0) + y[1] = h(p, 0) for n in 1:maxn-1 yp[n+1] = 0 for j = 1:n - yp[n+1] = yp[n+1] + generalized_binomials(j-1, n-1, α(t[n+1]), h)*f(t[j], y[j], v(ϕ, j, τ, h, y, yp, t)) + if iip + tmp = zeros(length(yp[1])) + f(tmp, y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j]) + yp[n+1] = yp[n+1]+generalized_binomials(j-1, n-1, order(t[n+1]), step_size)*tmp + else + yp[n+1] = yp[n+1] + generalized_binomials(j-1, n-1, order(t[n+1]), step_size)*f(y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j]) + end end - yp[n+1] = yp[n+1]/gamma(α(t[n+1]))+ϕ(0) + yp[n+1] = yp[n+1]/gamma(order(t[n+1]))+h(p, 0) y[n+1] = 0 @fastmath @inbounds @simd for j=1:n - y[n+1] = y[n+1]+a(j-1, n-1, α(t[n+1]), h)*f(t[j], y[j], v(ϕ, j, τ, h, y, yp, t)) + y[n+1] = y[n+1]+a(j-1, n-1, order(t[n+1]), step_size)*f(t[j], y[j], v(h, j, τ, step_size, y, yp, t, p)) end - y[n+1] = y[n+1]/gamma(α(t[n+1]))+h^α(t[n+1])*f(t[n+1], yp[n+1], v(ϕ, n+1, τ, h, y, yp, t))/gamma(α(t[n+1])+2)+ϕ(0) + y[n+1] = y[n+1]/gamma(order(t[n+1]))+step_size^order(t[n+1])*f(t[n+1], yp[n+1], v(h, n+1, τ, step_size, y, yp, t, p))/gamma(order(t[n+1])+2)+h(p, 0) end V = copy(t) @fastmath @inbounds @simd for n = 1:maxn-1 - V[n] = v(ϕ, n, τ, h, y, yp, t) + V[n] = v(h, n, τ, step_size, y, yp, t, p) end return V, y end -function solve_fdde_with_multiple_lags_and_variable_order(FDDE::FDDEProblem, h) - @unpack f, ϕ, α, τ, tspan = FDDE - t = collect(0:h:tspan) +function solve_fdde_with_multiple_lags_and_variable_order(FDDE::FDDEProblem, step_size) + @unpack f, h, order, constant_lags, p, tspan = FDDE + τ = constant_lags[1] + t = collect(0:step_size:tspan[2]) maxn = length(t) yp = zeros(maxn) y = copy(t) - y[1] = ϕ(0) + y[1] = h(p, 0) for n in 1:maxn-1 yp[n+1] = 0 for j = 1:n - yp[n+1] = yp[n+1]+generalized_binomials(j-1, n-1, α(t[n+1]), h)*f(t[j], y[j], multiv(ϕ, j, τ, h, y, yp)...) + yp[n+1] = yp[n+1]+generalized_binomials(j-1, n-1, order(t[n+1]), step_size)*f(t[j], y[j], multiv(h, j, τ, step_size, y, yp, p)...) end - yp[n+1] = yp[n+1]/gamma(α(t[n+1]))+ϕ(0) + yp[n+1] = yp[n+1]/gamma(order(t[n+1]))+h(p, 0) y[n+1] = 0 for j=1:n - y[n+1] = y[n+1]+multia(j-1, n-1, α(t[n+1]), h)*f(t[j], y[j], multiv(ϕ, j, τ, h, y, yp)...) + y[n+1] = y[n+1]+multia(j-1, n-1, order(t[n+1]), step_size)*f(t[j], y[j], multiv(h, j, τ, step_size, y, yp, p)...) end - y[n+1] = y[n+1]/gamma(α(t[n+1]))+h^α(t[n+1])*f(t[n+1], yp[n+1], multiv(ϕ, n+1, τ, h, y, yp)...)/gamma(α(t[n+1])+2) + ϕ(0) + y[n+1] = y[n+1]/gamma(order(t[n+1]))+step_size^order(t[n+1])*f(t[n+1], yp[n+1], multiv(h, n+1, τ, step_size, y, yp, p)...)/gamma(order(t[n+1])+2) + h(p, 0) end V = [] for n = 1:maxn - push!(V, multiv(ϕ, n, τ, h, y, yp)) + push!(V, multiv(h, n, τ, step_size, y, yp, p)) end delayed = zeros(Float, length(τ), length(V)) diff --git a/src/delay/product_integral.jl b/src/delay/product_integral.jl index 02ab32ca..7dba4617 100644 --- a/src/delay/product_integral.jl +++ b/src/delay/product_integral.jl @@ -1,7 +1,7 @@ #= # Usage - solve(FDDE::FDDEProblem, h, DelayPI()) + solve(FDDE::FDDEProblem, dt, DelayPI()) Use explicit rectangular product integration algorithm to solve an FDDE problem. @@ -20,17 +20,19 @@ Use explicit rectangular product integration algorithm to solve an FDDE problem. ``` =# -function solve(FDDE::FDDEProblem, h, ::PIEX) - @unpack f, ϕ, α, τ, tspan = FDDE +function solve(FDDE::FDDEProblem, dt, ::PIEX) + @unpack f, order, h, tspan, p, constant_lags = FDDE + τ = constant_lags[1] + iip = SciMLBase.isinplace(FDDE) t0 = tspan[1]; T = tspan[2] - N::Int = ceil(Int, (T-t0)/h) - t = t0 .+ h*collect(0:N) + N::Int = ceil(Int, (T-t0)/dt) + t = t0 .+ dt*collect(0:N) - nn_al = collect(Float64, 0:N).^α - b = [0; nn_al[2:end].-nn_al[1:end-1]]/gamma(α+1) - h_al = h^α + nn_al = collect(Float64, 0:N).^order + b = [0; nn_al[2:end].-nn_al[1:end-1]]/gamma(order+1) + h_al = dt^order - y0 = ϕ(t0) + y0 = h(p, t0) y = zeros(Float64, N+1) g = zeros(Float64, N+1) @@ -39,10 +41,10 @@ function solve(FDDE::FDDEProblem, h, ::PIEX) for n = 1:N tnm1 = t[n] if tnm1 <= τ - y_nm1_tau = ϕ(tnm1-τ) + y_nm1_tau = h(p, tnm1-τ) else - nm1_tau1 = floor(Int, (tnm1-τ)/h) - nm1_tau2 = ceil(Int, (tnm1-τ)/h) + nm1_tau1 = floor(Int, (tnm1-τ)/dt) + nm1_tau2 = ceil(Int, (tnm1-τ)/dt) if nm1_tau1 == nm1_tau2 y_nm1_tau = y[nm1_tau1+1] else @@ -54,7 +56,13 @@ function solve(FDDE::FDDEProblem, h, ::PIEX) end end - g[n] = f(tnm1, y[n], y_nm1_tau) + if iip + tmp = zeros(length(g[1])) + f(tmp, y[n], y_nm1_tau, tnm1, p) + g[n] = tmp + else + g[n] = f(y[n], y_nm1_tau, tnm1, p) + end f_mem = 0 for j = 0:n-1 f_mem = f_mem + g[j+1]*b[n-j+1] diff --git a/src/fodesystem/explicit_pi.jl b/src/fodesystem/explicit_pi.jl index 10a80c49..2544646d 100644 --- a/src/fodesystem/explicit_pi.jl +++ b/src/fodesystem/explicit_pi.jl @@ -5,21 +5,6 @@ Use explicit Product integration method to solve system of FODE. """ - -mutable struct M - an - bn - a0 - halpha1 - halpha2 - mu - mu_tol - r - index_fft - an_fft - bn_fft -end - function solve(prob::FODEProblem, h, ::PIEX) @unpack f, order, u0, tspan, p = prob @@ -230,9 +215,6 @@ function PIEX_system_triangolo(nxi, nxf, t, y, fy, zn, N, METH, problem_size, al return y, fy end - -sysf_vectorfield(t, y, f_fun) = f_fun(t, y) - function PIEX_system_starting_term(t, u0, m_alpha, t0, m_alpha_factorial) ys = zeros(size(u0, 1), 1) for k = 1 : maximum(m_alpha) diff --git a/src/fodesystem/newton_gregory.jl b/src/fodesystem/newton_gregory.jl index effc036d..1de3b0e7 100644 --- a/src/fodesystem/newton_gregory.jl +++ b/src/fodesystem/newton_gregory.jl @@ -285,12 +285,6 @@ function NG_weights(alpha, N) return omega, w, s end -#f_vectorfield(t, y, fdefun) = fdefun(zeros(length(y)), y, 0, t) -function Jf_vectorfield(t, y, Jfdefun) - f = Jfdefun(t, y) - return f -end - function NG_starting_term(t,y0, m_alpha, t0, m_alpha_factorial) ys = zeros(size(y0, 1), 1) for k = 1:Int64(m_alpha) diff --git a/src/fodesystem/trapezoid.jl b/src/fodesystem/trapezoid.jl index a45f6e59..ca2a0642 100644 --- a/src/fodesystem/trapezoid.jl +++ b/src/fodesystem/trapezoid.jl @@ -290,9 +290,6 @@ function TrapWeights(alpha, N) return omega, w, s end -#f_vectorfield(t, y, fdefun) = fdefun(zeros(length(y)), y, 0, t) -Jf_vectorfield(t, y, Jfdefun) = Jfdefun(t, y) - function TrapStartingTerm(t,y0, m_alpha, t0, m_alpha_factorial) ys = zeros(size(y0, 1), 1) for k = 1:m_alpha diff --git a/src/lyapunov.jl b/src/lyapunov.jl index 163fd98b..8843a399 100644 --- a/src/lyapunov.jl +++ b/src/lyapunov.jl @@ -144,14 +144,6 @@ function noncommensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out return FOLE(tspan, LE) end -function jacobian_of_fdefun(f, t, y, p) - ForwardDiff.jacobian(y) do y - du = similar(y) - f(du, y, p, t) - du - end -end - function commensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out, p) ne::Int = length(u0) # System dimension order = order[1] # Since this is the commensurate case, we only need one element in order array @@ -244,19 +236,6 @@ function commensurate_lyapunov(fun, order, t_start, h_norm, t_end, u0, h, out, p return FOLE(tspan, LE) end -mutable struct M - an - bn - a0 - halpha1 - halpha2 - mu - mu_tol - r - index_fft - an_fft - bn_fft -end #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) diff --git a/src/multitermsfode/implicit_pi_trapezoid.jl b/src/multitermsfode/implicit_pi_trapezoid.jl index 54aa2cf7..c433fa1f 100644 --- a/src/multitermsfode/implicit_pi_trapezoid.jl +++ b/src/multitermsfode/implicit_pi_trapezoid.jl @@ -214,8 +214,6 @@ function PITrap_triangolo(nxi, nxf, t, y, fy, zn, N, a0, an, t0, problem_size, u return y, fy end -Jf_vectorfield(t, y, fun)=fun(t, y) - function PITrap_startingterm_multi(t,t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val) ys = zeros(problem_size) diff --git a/src/types/problem_utils.jl b/src/types/problem_utils.jl new file mode 100644 index 00000000..b28839ec --- /dev/null +++ b/src/types/problem_utils.jl @@ -0,0 +1,28 @@ +function Base.summary(io::IO, prob::AbstractFDEProblem) + type_color, no_color = SciMLBase.get_colorizers(io) + print(io, + type_color, nameof(typeof(prob)), + no_color, " with uType ", + type_color, typeof(prob.u0), + no_color, " and tType ", + type_color, + prob.tspan isa Function ? + "Unknown" : (prob.tspan === nothing ? + "Nothing" : typeof(prob.tspan[1])), + no_color, ". In-place: ", + type_color, SciMLBase.isinplace(prob), + no_color) +end + +function Base.show(io::IO, mime::MIME"text\plain", A::FDDEProblem) + summary(io, A) + println(io) + print(io, "timespan: ") + show(io, mime, A.tspan) + println(io) + print(io, "order: ") + show(io, mime, A.order) + println(io) + print(io, "u0: ") + show(io, mime, A.u0) +end \ No newline at end of file diff --git a/src/types/problems.jl b/src/types/problems.jl index cb9a0272..a4459bcb 100644 --- a/src/types/problems.jl +++ b/src/types/problems.jl @@ -19,7 +19,7 @@ MultiTermsFODEProblem(parameters, orders, rightfun, u0, T) = MultiTermsFODEProbl 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 +SciMLBase.isinplace(prob::AbstractFDDEProblem{uType, tType, oType, lType, iip}) where {uType, tType, oType, lType, iip} = iip struct StandardFODEProblem end @@ -136,20 +136,6 @@ function FODEProblem(f, order, u0, tspan, p = SciMLBase.NullParameters(); kwargs end -""" - FDDEProblem(f, ϕ, α, τ, tspan) - -- `f`: The function describing fractional delay differential equations. -- `ϕ`: History function -Construct a fractional delayed differential equation problem. -""" -struct FDDEProblem <: FDEProblem - f::Function - ϕ::Union{Number, Function} - α::Union{Number, Function} - τ::Union{Number, AbstractArray, Function} - tspan::Union{Number, Tuple} -end #= #=FDDEProblem constructor=# function FDDEProblem(f::Function, @@ -165,6 +151,61 @@ end Construct system of fractional delay differential equations problem. """ + +struct StandardFDDEProblem end + +""" + FDDEProblem(f, ϕ, α, τ, tspan) + +- `f`: The function describing fractional delay differential equations. +- `ϕ`: History function +Construct a fractional delayed differential equation problem. +""" +struct FDDEProblem{uType, tType, oType, lType, isinplace, P, F, H, K, PT} <: + AbstractFDDEProblem{uType, tType, oType, lType, isinplace} + f::F + order::oType + u0::uType + h::H + tspan::tType + p::P + constant_lags::lType + kwargs::K + problem_type::PT + + SciMLBase.@add_kwonly function FDDEProblem{iip}(f::SciMLBase.AbstractDDEFunction{iip}, order, u0, h, tspan, + p = SciMLBase.NullParameters(); + constant_lags = (), + problem_type = StandardFDDEProblem(), + kwargs...) where {iip} + _u0 = SciMLBase.prepare_initial_state(u0) + _tspan = SciMLBase.promote_tspan(tspan) + SciMLBase.warn_paramtype(p) + new{typeof(_u0), typeof(_tspan), typeof(order), typeof(constant_lags), + isinplace(f), + typeof(p), typeof(f), typeof(h), typeof(kwargs), typeof(problem_type)}(f, order, _u0, + h, + _tspan, + p, + constant_lags, + kwargs, + problem_type) + end + + function FDDEProblem{iip}(f::SciMLBase.AbstractDDEFunction{iip}, h, tspan::Tuple, + p = SciMLBase.NullParameters(); kwargs...) where {iip} + FDDEProblem{iip}(f, h(p, first(tspan)), h, tspan, p; kwargs...) + end +end + +TruncatedStacktraces.@truncate_stacktrace FDDEProblem 5 1 2 + +FDDEProblem(f, args...; kwargs...) = FDDEProblem(DDEFunction(f), args...; kwargs...) + +function FDDEProblem(f::SciMLBase.AbstractDDEFunction, args...; kwargs...) + FDDEProblem{SciMLBase.isinplace(f)}(f, args...; kwargs...) +end + struct FDDESystem <: FDEProblem f::Function ϕ::AbstractArray diff --git a/test/FDDETests.jl b/test/FDDETests.jl index afe54032..6f082dbb 100644 --- a/test/FDDETests.jl +++ b/test/FDDETests.jl @@ -1,29 +1,28 @@ @testset "Test DelayPECE method with single constant lag" begin - function ϕ(x) - if x == 0 + function h(p, t) + if t == 0 return 19.00001 else return 19.0 end end - function f(t, y, ϕ) - return 3.5*y*(1-ϕ/19) - end + f(y, ϕ, p, t) = 3.5*y*(1-ϕ/19) - h = 0.5 + dt = 0.5 α = 0.97 - τ = 0.8 - T = 1 - fddeprob = FDDEProblem(f, ϕ, α, τ, T) - V, y = solve(fddeprob, h, DelayPECE()) + τ = [0.8] + u0 = [19.00001] + tspan = (0.0, 1.0) + fddeprob = FDDEProblem(f, α, 19.00001, h, constant_lags = [τ], tspan) + V, y = solve(fddeprob, dt, DelayPECE()) @test V≈[19.0, 19.0, 1.0] @test y≈[19.00001, 19.00001, 37.274176448220274] end @testset "Test DelayPECE method with single constant lag with variable order" begin - function ϕ(x) + function h(p, t) if x == 0 return 19.00001 else @@ -31,15 +30,12 @@ end end end - function f(t, y, ϕ) - return 3.5*y*(1-ϕ/19) - end - + f(y, ϕ, p, t) = 3.5*y*(1-ϕ/19) h = 0.5 alpha(t) = 0.99-(0.01/100)*t - τ = 0.8 - T = 1 - fddeprob = FDDEProblem(f, ϕ, alpha, τ, T) + τ = [0.8] + tspan = (0.0, 1.0) + fddeprob = FDDEProblem(f, alpha, h, τ, tspan) V, y = solve(fddeprob, h, DelayPECE()) @test V≈[19.0, 19.0, 1.0] @@ -47,50 +43,61 @@ end end @testset "Test DelayPECE method with single time varying lag" begin - function ϕ(x) - if x == 0 + function h(p, t) + if t == 0 return 19.00001 else return 19.0 end end - function f(t, y, ϕ) + function f(y, ϕ, p, t) return 3.5*y*(1-ϕ/19) end - h = 0.5 + dt = 0.5 alpha = 0.97 + u0 = [19.00001] τ(t) = 0.8 - T = 1 - fddeprob = FDDEProblem(f, ϕ, alpha, τ, T) - V, y = solve(fddeprob, h, DelayPECE()) + tspan = (0.0, 1.0) + fddeprob = FDDEProblem(f, alpha, u0, h, constant_lags = [τ], tspan) + V, y = solve(fddeprob, dt, DelayPECE()) @test V≈[19.0, 19.0, 1.0] @test y≈[19.00001, 19.00001, 37.274176448220274] end @testset "Test Product Integral method" begin - function ϕ(x) - if x == 0 + function ϕ(p, t) + if t == 0 return 19.00001 else return 19.0 end end - function f(t, y, ϕ) + function f(y, ϕ, p, t) return 3.5*y*(1-ϕ/19) end - prob = FDDEProblem(f, ϕ, 0.97, 0.8, (0, 2)) - result = solve(prob, 0.5, PIEX()) + τ = [0.8] + order = 0.97 + u0 = 19.00001 + tspan = (0.0, 2.0) + dt = 0.5 + prob = FDDEProblem(f, order, u0, ϕ, constant_lags = τ, tspan) + result = solve(prob, dt, PIEX()) @test result≈[19.00001, 19.00001, 19.00001, 18.99999190949352, 18.99997456359874] end @testset "Test DelayABM method" begin - h=0.5; T=5; α=0.97; τ=2; q=0.5 - delayabmfun(t, ϕ, y) = 2*ϕ/(1+ϕ^9.65)-y - prob = FDDEProblem(delayabmfun, q, α, τ, T) - x, y=solve(prob, h, DelayABM()) + dt=0.5 + tspan = (0.0, 5.0) + α=0.97 + τ=[2] + h = 0.5 + u0 = 0.5 + delayabmfun(ϕ, y, p, t) = 2*ϕ/(1+ϕ^9.65)-y + prob = FDDEProblem(delayabmfun, α, u0, h, constant_lags=τ, tspan) + x, y=solve(prob, dt, DelayABM()) @test isapprox(x, [1.078559863692747, 1.175963999045738, 1.1661317460354588, 1.128481756921719, 1.0016061526083417, 0.7724564325042358, 0.5974978685646778]; atol=1e-3) @test isapprox(y, [0.8889787467894421, 0.9404487875504524, 0.9667449499617093, 0.9803311436135411, 1.078559863692747, 1.175963999045738, 1.1661317460354588]; atol=1e-3) @@ -173,7 +180,7 @@ end 0.401838 0.401838 0.60134 0.60134]; atol=1e-4) end - +#= @testset "Test DelayPECE method with multiple constant lags" begin α = 0.95; ϕ(x) = 0.5 τ = [2, 2.6] @@ -195,4 +202,5 @@ end @test isapprox(delayed, [0.5 0.5 0.5 0.5 0.5 1.12259 0.62302 1.28025 0.818417 2.53596 -1.65415 0.5 0.5 0.5 0.5 0.5 0.605999 1.2225 0.491575 1.37261 0.47491 3.37398]; atol=1e-3) -end \ No newline at end of file +end +=# \ No newline at end of file From 1a8a696e115ef43be723e9f139291420ed7a5b1f Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Tue, 28 Nov 2023 22:28:00 +0800 Subject: [PATCH 02/14] Remove version bound for compat Signed-off-by: ErikQQY <2283984853@qq.com> --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 6495a195..af2afc6a 100644 --- a/Project.toml +++ b/Project.toml @@ -26,8 +26,8 @@ TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] -ApproxFun = "0.13.25" -ArrayInterface = "7.6" +ApproxFun = "0.13" +ArrayInterface = "7" DiffEqBase = "6" FFTW = "1.0.0" ForwardDiff = "0.10" From e2e6644483e06a3dff1b1266e8880a43ff2b7beb Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Tue, 28 Nov 2023 22:46:59 +0800 Subject: [PATCH 03/14] Fix CI errors Signed-off-by: ErikQQY <2283984853@qq.com> --- src/types/problems.jl | 5 ++--- src/utils.jl | 9 --------- test/auxillary.jl | 20 ++++++++++++-------- 3 files changed, 14 insertions(+), 20 deletions(-) diff --git a/src/types/problems.jl b/src/types/problems.jl index a4459bcb..5c96fb64 100644 --- a/src/types/problems.jl +++ b/src/types/problems.jl @@ -178,12 +178,11 @@ struct FDDEProblem{uType, tType, oType, lType, isinplace, P, F, H, K, PT} <: constant_lags = (), problem_type = StandardFDDEProblem(), kwargs...) where {iip} - _u0 = SciMLBase.prepare_initial_state(u0) _tspan = SciMLBase.promote_tspan(tspan) SciMLBase.warn_paramtype(p) - new{typeof(_u0), typeof(_tspan), typeof(order), typeof(constant_lags), + new{typeof(u0), typeof(_tspan), typeof(order), typeof(constant_lags), isinplace(f), - typeof(p), typeof(f), typeof(h), typeof(kwargs), typeof(problem_type)}(f, order, _u0, + typeof(p), typeof(f), typeof(h), typeof(kwargs), typeof(problem_type)}(f, order, u0, h, _tspan, p, diff --git a/src/utils.jl b/src/utils.jl index 2ca5e1ea..91feb141 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -14,15 +14,6 @@ function Base.show(io::IO, prob::MultiTermsFODEProblem) println("u0: $(prob.u0)") end -function Base.show(io::IO, prob::FDDEProblem) - printstyled(typeof(prob), color=:light_blue) - printstyled(" with order ") - printstyled("$(prob.α)", color=:red) - println() - println("timespan: $(prob.tspan)") - println("history function: $(prob.ϕ)") -end - function Base.show(io::IO, prob::FDDESystem) printstyled(typeof(prob), color=:light_blue) printstyled(" with order ") diff --git a/test/auxillary.jl b/test/auxillary.jl index b6d9092b..758e5319 100644 --- a/test/auxillary.jl +++ b/test/auxillary.jl @@ -80,19 +80,23 @@ end end @testset "Test FDDEProblem show method" begin - # FDDEProblem - function ϕ(x) - if x == 0 + function h(p, t) + if t == 0 return 19.00001 else return 19.0 end end - function delayf(t, y, ϕ) - return 3.5*y*(1-ϕ/19) - end - fddeprob = FDDEProblem(delayf, ϕ, 0.97, 0.8, (0, 2)) - fddesol = solve(fddeprob, 0.5, PIEX()) + + f(y, ϕ, p, t) = 3.5*y*(1-ϕ/19) + + dt = 0.5 + α = 0.97 + τ = [0.8] + u0 = [19.00001] + tspan = (0.0, 1.0) + fddeprob = FDDEProblem(f, α, 19.00001, h, constant_lags = [τ], tspan) + fddesol = solve(fddeprob, dt, DelayPECE()) @test_nowarn show(fddeprob) @test_nowarn show(fddesol) From 7f05492c56b3aafcc393a3098e08aeb3d4277a24 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Tue, 28 Nov 2023 23:00:40 +0800 Subject: [PATCH 04/14] Fix CI errors Signed-off-by: ErikQQY <2283984853@qq.com> --- src/delay/DelayABM.jl | 14 +++++++------- test/FDDETests.jl | 13 +++++++------ 2 files changed, 14 insertions(+), 13 deletions(-) diff --git a/src/delay/DelayABM.jl b/src/delay/DelayABM.jl index c50ed5a4..689744b2 100644 --- a/src/delay/DelayABM.jl +++ b/src/delay/DelayABM.jl @@ -85,8 +85,8 @@ DelayABM method for system of fractional delay differential equations. =# function solve(FDDESys::FDDESystem, dt, ::DelayABM) - @unpack f, h, order, τ, T = FDDESys - len = length(h) + @unpack f, ϕ, α, τ, T = FDDESys + len = length(ϕ) N::Int = round(Int, T/dt) Ndelay = round(Int, τ/dt) t = collect(Float64, 0:dt:(T-τ)) @@ -96,19 +96,19 @@ function solve(FDDESys::FDDESystem, dt, ::DelayABM) # Put the delay term in the array for i=1:len - x[1:Ndelay, i] = h[i]*ones(Ndelay) + x[1:Ndelay, i] = ϕ[i]*ones(Ndelay) end x0 = copy(x[Ndelay, :]) for i=1:len - αi = order[i] + αi = α[i] f(du, x0, x[1, :], 0) x1[Ndelay+1, i] = x0[i] + dt^αi*du[i]/(gamma(αi)*αi) end for i=1:len - αi = order[i] + αi = α[i] f(du, x[Ndelay+1, :], x[1, :], 0) du1 = copy(du) f(du, x0, x[1, :], 0) @@ -118,7 +118,7 @@ function solve(FDDESys::FDDESystem, dt, ::DelayABM) M1=zeros(len); N1=zeros(len) @fastmath @inbounds @simd for n=1:N for i=1:len - αi = order[i] + αi = α[i] f(du, x0, x[1, :], 0) M1[i]=(n^(αi+1)-(n-αi)*(n+1)^αi)*du[i] N1[i]=((n+1)^αi-n^αi)*du[i] @@ -134,7 +134,7 @@ function solve(FDDESys::FDDESystem, dt, ::DelayABM) end for i=1:len - αi = order[i] + αi = α[i] x1[Ndelay+n+1, i] = x0[i]+dt^αi*N1[i]/(gamma(αi)*αi) f(du, x[Ndelay+n+1, :], x[n+1, :], 0) x[Ndelay+n+1, i] = x0[i]+dt^αi*(du[i]+M1[i])/gamma(αi+2) diff --git a/test/FDDETests.jl b/test/FDDETests.jl index 6f082dbb..dc9e5189 100644 --- a/test/FDDETests.jl +++ b/test/FDDETests.jl @@ -12,9 +12,9 @@ dt = 0.5 α = 0.97 τ = [0.8] - u0 = [19.00001] + u0 = 19.00001 tspan = (0.0, 1.0) - fddeprob = FDDEProblem(f, α, 19.00001, h, constant_lags = [τ], tspan) + fddeprob = FDDEProblem(f, α, u0, h, constant_lags = τ, tspan) V, y = solve(fddeprob, dt, DelayPECE()) @test V≈[19.0, 19.0, 1.0] @@ -31,12 +31,13 @@ end end f(y, ϕ, p, t) = 3.5*y*(1-ϕ/19) - h = 0.5 + dt = 0.5 alpha(t) = 0.99-(0.01/100)*t τ = [0.8] tspan = (0.0, 1.0) - fddeprob = FDDEProblem(f, alpha, h, τ, tspan) - V, y = solve(fddeprob, h, DelayPECE()) + u0 = 19.00001 + fddeprob = FDDEProblem(f, alpha, u0, h, constant_lags = τ, tspan) + V, y = solve(fddeprob, dt, DelayPECE()) @test V≈[19.0, 19.0, 1.0] @test y≈[19.00001, 19.00001, 36.698681913021375] @@ -57,7 +58,7 @@ end dt = 0.5 alpha = 0.97 - u0 = [19.00001] + u0 = 19.00001 τ(t) = 0.8 tspan = (0.0, 1.0) fddeprob = FDDEProblem(f, alpha, u0, h, constant_lags = [τ], tspan) From 35f1c43ef7fec4b6fdc77197f7d15b37ac7f089e Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Tue, 28 Nov 2023 23:03:34 +0800 Subject: [PATCH 05/14] Fix FDDEProblem sohw test Signed-off-by: ErikQQY <2283984853@qq.com> --- test/auxillary.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/auxillary.jl b/test/auxillary.jl index 758e5319..b4be5ac2 100644 --- a/test/auxillary.jl +++ b/test/auxillary.jl @@ -93,9 +93,9 @@ end dt = 0.5 α = 0.97 τ = [0.8] - u0 = [19.00001] + u0 = 19.00001 tspan = (0.0, 1.0) - fddeprob = FDDEProblem(f, α, 19.00001, h, constant_lags = [τ], tspan) + fddeprob = FDDEProblem(f, α, u0, h, constant_lags = τ, tspan) fddesol = solve(fddeprob, dt, DelayPECE()) @test_nowarn show(fddeprob) From afcf9d8179981906a1b2cae66f818ec9c216ce63 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Tue, 28 Nov 2023 23:38:30 +0800 Subject: [PATCH 06/14] Fix typos Signed-off-by: ErikQQY <2283984853@qq.com> --- src/delay/DelayABM.jl | 2 +- src/delay/DelayPECE.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/delay/DelayABM.jl b/src/delay/DelayABM.jl index 689744b2..69bab5e6 100644 --- a/src/delay/DelayABM.jl +++ b/src/delay/DelayABM.jl @@ -126,7 +126,7 @@ function solve(FDDESys::FDDESystem, dt, ::DelayABM) @fastmath @inbounds @simd for j=1:n for i=1:len - αi = order[i] + αi = α[i] f(du, x[Ndelay+j, :], x[j, :], 0) M1[i] = M1[i]+((n-j+2)^(αi+1)+(n-j)^(αi+1)-2*(n-j+1)^(αi+1))*du[i] N1[i] = N1[i]+((n-j+1)^αi-(n-j)^αi)*du[i] diff --git a/src/delay/DelayPECE.jl b/src/delay/DelayPECE.jl index 6299337d..382c7221 100644 --- a/src/delay/DelayPECE.jl +++ b/src/delay/DelayPECE.jl @@ -44,7 +44,7 @@ function solve(FDDE::FDDEProblem, step_size, ::DelayPECE) solve_fdde_with_multiple_lags(FDDE, step_size) #TODO: implement this # Varying order fractional delay differential equations elseif FDDE.order isa Function - if length(FDDE.cpnstant_lags[1]) == 1 + if length(FDDE.constant_lags[1]) == 1 # Here is the PECE solver for single lag with variable order solve_fdde_with_single_lag_and_variable_order(FDDE, step_size) else From 30590c6fe55ee418415dca0b55c5576d435754be Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Wed, 29 Nov 2023 01:02:55 +0800 Subject: [PATCH 07/14] Fix variable order test case Signed-off-by: ErikQQY <2283984853@qq.com> --- test/FDDETests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/FDDETests.jl b/test/FDDETests.jl index dc9e5189..594e2d35 100644 --- a/test/FDDETests.jl +++ b/test/FDDETests.jl @@ -36,7 +36,7 @@ end τ = [0.8] tspan = (0.0, 1.0) u0 = 19.00001 - fddeprob = FDDEProblem(f, alpha, u0, h, constant_lags = τ, tspan) + fddeprob = FDDEProblem(f, [alpha], u0, h, constant_lags = τ, tspan) V, y = solve(fddeprob, dt, DelayPECE()) @test V≈[19.0, 19.0, 1.0] From 63ef508d39a4377a76861fd4da54813e4f8017c6 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Wed, 29 Nov 2023 01:21:46 +0800 Subject: [PATCH 08/14] Fix variable order test case Signed-off-by: ErikQQY <2283984853@qq.com> --- test/FDDETests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/FDDETests.jl b/test/FDDETests.jl index 594e2d35..5cbe5989 100644 --- a/test/FDDETests.jl +++ b/test/FDDETests.jl @@ -23,7 +23,7 @@ end @testset "Test DelayPECE method with single constant lag with variable order" begin function h(p, t) - if x == 0 + if t == 0 return 19.00001 else return 19.0 From 709bc9ac72882afd137a396380efb4892e07b922 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Wed, 29 Nov 2023 10:43:18 +0800 Subject: [PATCH 09/14] Fix variable order test case Signed-off-by: ErikQQY <2283984853@qq.com> --- src/delay/DelayPECE.jl | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/src/delay/DelayPECE.jl b/src/delay/DelayPECE.jl index 382c7221..1b34f663 100644 --- a/src/delay/DelayPECE.jl +++ b/src/delay/DelayPECE.jl @@ -43,7 +43,7 @@ function solve(FDDE::FDDEProblem, step_size, ::DelayPECE) # Here is the PECE solver for multiple time varying lags solve_fdde_with_multiple_lags(FDDE, step_size) #TODO: implement this # Varying order fractional delay differential equations - elseif FDDE.order isa Function + elseif FDDE.order[1] isa Function if length(FDDE.constant_lags[1]) == 1 # Here is the PECE solver for single lag with variable order solve_fdde_with_single_lag_and_variable_order(FDDE, step_size) @@ -127,6 +127,7 @@ function a(j::Int, n::Int, order, step_size) return result*step_size^order / (order*(order + 1)) end +#generalized_binomials(j, n, order::Number, step_size) = @. step_size^order/order*((n-j+1)^order - (n-j)^order) generalized_binomials(j, n, order, step_size) = @. step_size^order/order*((n-j+1)^order - (n-j)^order) function v(h, n, τ, step_size, y, yp, t, p) @@ -252,7 +253,7 @@ function solve_fdde_with_single_lag_and_variable_order(FDDE::FDDEProblem, step_s if iip tmp = zeros(length(yp[1])) f(tmp, y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j]) - yp[n+1] = yp[n+1]+generalized_binomials(j-1, n-1, order(t[n+1]), step_size)*tmp + yp[n+1] = yp[n+1] + generalized_binomials(j-1, n-1, order(t[n+1]), step_size)*tmp else yp[n+1] = yp[n+1] + generalized_binomials(j-1, n-1, order(t[n+1]), step_size)*f(y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j]) end @@ -262,10 +263,22 @@ function solve_fdde_with_single_lag_and_variable_order(FDDE::FDDEProblem, step_s y[n+1] = 0 @fastmath @inbounds @simd for j=1:n - y[n+1] = y[n+1]+a(j-1, n-1, order(t[n+1]), step_size)*f(t[j], y[j], v(h, j, τ, step_size, y, yp, t, p)) + if iip + tmp = zeros(length(y[1])) + f(tmp, y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j]) + y[n+1] = y[n+1]+a(j-1, n-1, order(t[n+1]), step_size)*tmp + else + y[n+1] = y[n+1]+a(j-1, n-1, order(t[n+1]), step_size)*f(y[j], v(h, j, τ, step_size, y, yp, t, p), p, t[j]) + end end - y[n+1] = y[n+1]/gamma(order(t[n+1]))+step_size^order(t[n+1])*f(t[n+1], yp[n+1], v(h, n+1, τ, step_size, y, yp, t, p))/gamma(order(t[n+1])+2)+h(p, 0) + if iip + tmp = zeros(length(y[1])) + f(tmp, yp[n+1], v(h, n+1, τ, step_size, y, yp, t, p), p, t[n+1]) + y[n+1] = y[n+1]/gamma(order(t[n+1]))+step_size^order(t[n+1])*tmp/gamma(order(t[n+1])+2)+h(p, 0) + else + y[n+1] = y[n+1]/gamma(order(t[n+1]))+step_size^order(t[n+1])*f(yp[n+1], v(h, n+1, τ, step_size, y, yp, t, p), p, t[n+1])/gamma(order(t[n+1])+2)+h(p, 0) + end end V = copy(t) From d87ec38e73d7d12dac44a2a11537b4e50e6ee7f5 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Wed, 29 Nov 2023 10:58:30 +0800 Subject: [PATCH 10/14] Fix doc failure Signed-off-by: ErikQQY <2283984853@qq.com> --- docs/make.jl | 1 + src/ffode/atangana_seda.jl | 3 +++ src/fodesystem/Euler.jl | 3 +++ src/fodesystem/PIPECE.jl | 3 +++ 4 files changed, 10 insertions(+) diff --git a/docs/make.jl b/docs/make.jl index 10efcb36..65612b4d 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -12,6 +12,7 @@ makedocs(; canonical="https://SciFracX.github.io/FractionalDiffEq.jl", assets=["assets/favicon.ico"], ), + warnonly = :doctest, pages = [ "FractionalDiffEq.jl" => "index.md", "Get Started" => "get_started.md", diff --git a/src/ffode/atangana_seda.jl b/src/ffode/atangana_seda.jl index e6f0c9a6..c625bf20 100644 --- a/src/ffode/atangana_seda.jl +++ b/src/ffode/atangana_seda.jl @@ -1,3 +1,6 @@ +""" +AtanganaSeda methods for +""" struct AtanganaSeda <: FODESystemAlgorithm end function solve(prob::Union{FFMODEProblem, FFMODESystem}, h::Float64, ::AtanganaSeda) diff --git a/src/fodesystem/Euler.jl b/src/fodesystem/Euler.jl index f5484641..c8e92168 100644 --- a/src/fodesystem/Euler.jl +++ b/src/fodesystem/Euler.jl @@ -1,3 +1,6 @@ +""" +The classical Euler method extended for fractional order differential equations. +""" struct Euler <: FODESystemAlgorithm end function solve(prob::FODEProblem, h, ::Euler) diff --git a/src/fodesystem/PIPECE.jl b/src/fodesystem/PIPECE.jl index 794c80c5..1b9497fd 100644 --- a/src/fodesystem/PIPECE.jl +++ b/src/fodesystem/PIPECE.jl @@ -30,6 +30,9 @@ mutable struct M bn_fft end +""" +Predictor-Correct scheme. +""" struct PECE <: FODESystemAlgorithm end #TODO: Rename as PIPECE From 0a6450acd8b27aa844e41235da2c1e31cfc5af14 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Wed, 29 Nov 2023 11:09:12 +0800 Subject: [PATCH 11/14] Fix doc failure Signed-off-by: ErikQQY <2283984853@qq.com> --- docs/src/algorithms.md | 68 ++++++++++++++++++------------------------ 1 file changed, 29 insertions(+), 39 deletions(-) diff --git a/docs/src/algorithms.md b/docs/src/algorithms.md index 2035ba43..df827a35 100644 --- a/docs/src/algorithms.md +++ b/docs/src/algorithms.md @@ -8,56 +8,46 @@ Pages = ["algorithms.md"] ### Single Term FODE -```@docs -FractionalDiffEq.PECE -FractionalDiffEq.Euler -FractionalDiffEq.PIEX -FractionalDiffEq.AtanganaSeda -``` +PECE +Euler +PIEX +AtanganaSeda + ### Multi-Term FODE -```@docs -FractionalDiffEq.FODEMatrixDiscrete -FractionalDiffEq.ClosedForm -FractionalDiffEq.ClosedFormHankelM -FractionalDiffEq.PIPECE -FractionalDiffEq.PIRect -FractionalDiffEq.PITrap -``` +FODEMatrixDiscrete +ClosedForm +ClosedFormHankelM +PIPECE +PIRect +PITrap ### System of FODE -```@docs -FractionalDiffEq.NonLinearAlg -FractionalDiffEq.PECE -FractionalDiffEq.GL -FractionalDiffEq.FLMMBDF -FractionalDiffEq.FLMMNewtonGregory -FractionalDiffEq.FLMMTrap -FractionalDiffEq.PIEX -FractionalDiffEq.NewtonPolynomial -FractionalDiffEq.AtanganaSedaAB -``` +NonLinearAlg +PECE +GL +FLMMBDF +FLMMNewtonGregory +FLMMTrap +PIEX +NewtonPolynomial +AtanganaSedaAB -## Fractional Delay Differential Equatinos +## Fractional Delay Differential Equations + +DelayPECE +DelayPI +MatrixForm +DelayABM -```@docs -FractionalDiffEq.DelayPECE -FractionalDiffEq.DelayPI -FractionalDiffEq.MatrixForm -FractionalDiffEq.DelayABM -``` ## Distributed Order Differential Equations -```@docs -FractionalDiffEq.DOMatrixDiscrete -``` +DOMatrixDiscrete ## Fractional Differences Equations -```@docs -FractionalDiffEq.PECEDifference -FractionalDiffEq.GL -``` \ No newline at end of file +PECEDifference +GL \ No newline at end of file From 7590039efbaee544b75c5122803c84411f6ab446 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Wed, 29 Nov 2023 11:19:48 +0800 Subject: [PATCH 12/14] Fix doc failure Signed-off-by: ErikQQY <2283984853@qq.com> --- docs/make.jl | 1 + docs/src/problems.md | 6 ------ src/types/problems.jl | 13 +++++++++---- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 65612b4d..2c61a266 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -13,6 +13,7 @@ makedocs(; assets=["assets/favicon.ico"], ), warnonly = :doctest, + doctest = false, pages = [ "FractionalDiffEq.jl" => "index.md", "Get Started" => "get_started.md", diff --git a/docs/src/problems.md b/docs/src/problems.md index f54dcbf0..389e7c26 100644 --- a/docs/src/problems.md +++ b/docs/src/problems.md @@ -36,12 +36,6 @@ FractionalDiffEq.FFEODEProblem FractionalDiffEq.FFMODEProblem ``` -## FPDE Problem - -```@docs -FractionalDiffEq.FPDEProblem -``` - ## FDDE Problem ```@docs diff --git a/src/types/problems.jl b/src/types/problems.jl index 5c96fb64..e81305e4 100644 --- a/src/types/problems.jl +++ b/src/types/problems.jl @@ -3,6 +3,10 @@ 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 + +""" +Multiple terms fractional order differential equations. +""" struct MultiTermsFODEProblem <: FDEProblem parameters orders @@ -146,11 +150,7 @@ function FDDEProblem(f::Function, return FDDEProblem(f, ϕ, α, τ, T) end =# -""" - FDDESystem(f, ϕ, α, τ, T) -Construct system of fractional delay differential equations problem. -""" struct StandardFDDEProblem end @@ -205,6 +205,11 @@ function FDDEProblem(f::SciMLBase.AbstractDDEFunction, args...; kwargs...) FDDEProblem{SciMLBase.isinplace(f)}(f, args...; kwargs...) end +""" + FDDESystem(f, ϕ, α, τ, T) + +Construct system of fractional delay differential equations problem. +""" struct FDDESystem <: FDEProblem f::Function ϕ::AbstractArray From 498dcdd4519386bc316072e07618f167c498c963 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Wed, 29 Nov 2023 11:28:39 +0800 Subject: [PATCH 13/14] Fix doc failure Signed-off-by: ErikQQY <2283984853@qq.com> --- docs/src/problems.md | 44 ++++++++++++++------------------------------ 1 file changed, 14 insertions(+), 30 deletions(-) diff --git a/docs/src/problems.md b/docs/src/problems.md index 389e7c26..04985366 100644 --- a/docs/src/problems.md +++ b/docs/src/problems.md @@ -4,15 +4,11 @@ Pages = ["problems.md"] ``` -The general fractional differential equations problem type: - -```@docs -FractionalDiffEq.FDEProblem -``` +The general fractional differential equations problem type. ## FODE Problem -Fractional ordinary problems type, we can classified them as Single-Term and Multi-Term problems: +Fractional ordinary problems type, we can classify them as Single-Term and Multi-Term problems: ```math D^{\alpha}y(t)=f(t, y) @@ -22,43 +18,31 @@ D^{\alpha}y(t)=f(t, y) \sum_{s=0}^pc_sD^{\beta_s}y=f ``` -```@docs -FractionalDiffEq.SingleTermFODEProblem -FractionalDiffEq.MultiTermsFODEProblem -FractionalDiffEq.FODESystem -``` +SingleTermFODEProblem +MultiTermsFODEProblem +FODESystem ### FFODE Problem -```@docs -FractionalDiffEq.FFPODEProblem -FractionalDiffEq.FFEODEProblem -FractionalDiffEq.FFMODEProblem -``` +FFPODEProblem +FFEODEProblem +FFMODEProblem ## FDDE Problem -```@docs -FractionalDiffEq.FDDEProblem -FractionalDiffEq.FDDEMatrixProblem -FractionalDiffEq.FDDESystem -``` +FDDEProblem +FDDEMatrixProblem +FDDESystem ## DODE Problem -```@docs FractionalDiffEq.DODEProblem -``` ## Fractional Difference Problem -```@docs -FractionalDiffEq.FractionalDifferenceProblem -FractionalDiffEq.FractionalDifferenceSystem -``` +FractionalDifferenceProblem +FractionalDifferenceSystem ## FIE Problem -```@docs -FractionalDiffEq.FIEProblem -``` \ No newline at end of file +FIEProblem From 1124d3ba7ee2a4c5804a652f562914d1b9fdefb1 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Wed, 29 Nov 2023 11:44:28 +0800 Subject: [PATCH 14/14] Fix doc failure Signed-off-by: ErikQQY <2283984853@qq.com> --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index 2c61a266..97c9e54a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -12,7 +12,7 @@ makedocs(; canonical="https://SciFracX.github.io/FractionalDiffEq.jl", assets=["assets/favicon.ico"], ), - warnonly = :doctest, + warnonly = [:missing_docs, :cross_references], doctest = false, pages = [ "FractionalDiffEq.jl" => "index.md",