Skip to content

Commit

Permalink
Several fix for FDDE solvers
Browse files Browse the repository at this point in the history
  • Loading branch information
@ErikQQY
ErikQQY committed Jan 3, 2025
1 parent 99abd37 commit f8bdf65
Show file tree
Hide file tree
Showing 7 changed files with 104 additions and 87 deletions.
138 changes: 70 additions & 68 deletions src/delay/pece.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
p
y
yp
L

dt
kwargs
Expand Down Expand Up @@ -52,28 +53,18 @@ function SciMLBase.__init(prob::FDDEProblem, alg::DelayPECE; dt = 0.0, kwargs...
t0 = tspan[1]
tfinal = tspan[2]
T = eltype(u0)
mesh = collect(Float64, t0:dt:tfinal)
size(mesh, 1)
yp = collect(Float64, t0:dt:(tfinal + dt))
N = length(t0:dt:(tfinal + dt))
yp = _generate_similar_array(u0, N, u0)
y = _generate_similar_array(u0, N - 1, u0)
y[1] = (length(u0) == 1) ? _generate_similar_array(u0, 1, h(p, 0)) : u0#_generate_similar_array(u0, 1, h(p, 0))

return DelayPECECache{iip, T}(prob, alg, mesh, u0, order[1], τ, p, y, yp, dt, kwargs)
N = Int(cld(tfinal - t0, dt))
L = length(collect(-τ:dt:t0))
mesh = collect(range(t0; stop = tfinal, length = N + 1))
yp = _generate_similar_array(u0, N+1, u0)
y = _generate_similar_array(u0, L+N+1, u0)
y[L:-1:1] .= ifelse((length(u0) == 1), map(x->[h(p, x)], collect(-τ:dt:t0)), map(x->h(p, x), collect(-τ:dt:(t0))))
y[L+1] = ifelse((length(u0) == 1), [u0], u0)#_generate_similar_array(u0, 1, h(p, 0))

return DelayPECECache{iip, T}(prob, alg, mesh, u0, order, τ, p, y, yp, L, dt, kwargs)
end

@inline function _generate_similar_array(u0, N, src)
if N ≠ 1
if length(u0) == 1
return [similar([src]) for i in 1:N]
else
return [similar(src) for i in 1:N]
end
else
return [src]
end
end
@inline _generate_similar_array(u0, N, src) = ifelse(length(u0) == 1, [[zero(u0)] for _ in 1:N], [zero(src) for _ in 1:N])

@inline function OrderWrapper(order, t)
if order isa Function
Expand All @@ -84,65 +75,70 @@ end
end

function SciMLBase.solve!(cache::DelayPECECache{iip, T}) where {iip, T}
(; prob, alg, mesh, u0, order, p, dt) = cache
maxn = length(mesh)
if length(cache.constant_lags) == 1
return solve_fdde_with_single_lag(cache)
else
return solve_fdde_with_multiple_lags(cache)
end
end

function solve_fdde_with_single_lag(cache::DelayPECECache{iip, T}) where {iip, T}
(; prob, alg, mesh, u0, order, L, p, dt) = cache
N = length(mesh)
l = length(u0)
initial = (length(u0) == 1) ? _generate_similar_array(u0, 1, prob.h(p, 0)) : u0
for n in 1:(maxn - 1)
order = OrderWrapper(order, mesh[n + 1])
cache.yp[n + 1] = zeros(T, l)
for j in 1:n
initial = ifelse(l == 1, [u0], u0)
tmp = zeros(T, l)
for i in 0:N-2
order = OrderWrapper(order, mesh[i + 1])
# compute the yp part
for j in 0:i
if iip
tmp = zeros(length(cache.yp[1]))
prob.f(tmp, cache.y[j], v(cache, j), p, mesh[j])
cache.yp[n + 1] = cache.yp[n + 1] +
generalized_binomials(j - 1, n - 1, order, dt) * tmp
prob.f(tmp, cache.y[L+j+1], v(cache, j), p, mesh[j + 1])
cache.yp[i + 1] = cache.yp[i + 1] +
generalized_binomials(j, i, order, dt) * tmp
else
length(u0) == 1 ?
(tmp = prob.f(cache.y[j][1], v(cache, j)[1], p, mesh[j])) :
(tmp = prob.f(cache.y[j], v(cache, j), p, mesh[j]))
@. cache.yp[n + 1] = cache.yp[n + 1] +
generalized_binomials(j - 1, n - 1, order, dt) * tmp
(tmp = prob.f(cache.y[L+j+1][1], v(cache, j)[1], p, mesh[j + 1])) :
(tmp = prob.f(cache.y[L+j+1], v(cache, j), p, mesh[j+1]))
@. cache.yp[i + 1] = cache.yp[i + 1] +
generalized_binomials(j, i, order, dt) * tmp
end
end
cache.yp[n + 1] = cache.yp[n + 1] / gamma(order) + initial
cache.yp[i + 1] = cache.yp[i + 1] / gamma(order) + initial

cache.y[n + 1] = zeros(T, l)

@fastmath @inbounds @simd for j in 1:n
# compute the a part
for j in 0:i
if iip
tmp = zeros(T, length(cache.y[1]))
prob.f(tmp, cache.y[j], v(cache, j), p, mesh[j])
cache.y[n + 1] = cache.y[n + 1] + a(j - 1, n - 1, order, dt) * tmp
prob.f(tmp, cache.y[L + j + 1], v(cache, j), p, mesh[j+1])
cache.y[L + i + 2] += a(j, i, order, dt) * tmp
else
length(u0) == 1 ?
(tmp = prob.f(cache.y[j][1], v(cache, j)[1], p, mesh[j])) :
(tmp = prob.f(cache.y[j], v(cache, j), p, mesh[j]))
@. cache.y[n + 1] = cache.y[n + 1] + a(j - 1, n - 1, order, dt) * tmp
(tmp = prob.f(cache.y[L + j + 1][1], v(cache, j)[1], p, mesh[j+1])) :
(tmp = prob.f(cache.y[L + j + 1], v(cache, j), p, mesh[j]))
@. cache.y[L + i + 2] += a(j, i, order, dt) * tmp
end
end

if iip
tmp = zeros(T, l)
prob.f(tmp, cache.yp[n + 1], v(cache, n + 1), p, mesh[n + 1])
cache.y[n + 1] = cache.y[n + 1] / gamma(order) +
dt^order * tmp / gamma(order + 2) +
initial
prob.f(tmp, cache.yp[i + 1], v(cache, i + 1), p, mesh[i + 2])
cache.y[L + i + 2] = cache.y[L + i + 2] / gamma(order) +
dt^order * tmp / gamma(order + 2) + initial
else
length(u0) == 1 ?
(tmp = prob.f(cache.yp[n + 1][1], v(cache, n + 1)[1], p, mesh[n + 1])) :
(tmp = prob.f(cache.yp[n + 1], v(cache, n + 1), p, mesh[n + 1]))
@. cache.y[n + 1] = cache.y[n + 1] / gamma(order) +
dt^order * tmp / gamma(order + 2) +
initial
(tmp = prob.f(cache.yp[i + 1][1], v(cache, i + 1)[1], p, mesh[i + 2])) :
(tmp = prob.f(cache.yp[i + 1], v(cache, i + 1), p, mesh[i + 2]))
@. cache.y[L + i + 2] = cache.y[L + i + 2] / gamma(order) +
dt^order * tmp / gamma(order + 2) + initial
end
end

return DiffEqBase.build_solution(prob, alg, mesh, cache.y)
u = cache.y[L+1:end]

return DiffEqBase.build_solution(prob, alg, mesh, u)
end

function a(j::I, n::I, order, dt) where {I <: Integer}
if j == n + 1
if j == n
result = 1
elseif j == 0
result = n^(order + 1) - (n - order) * (n + 1)^order
Expand All @@ -153,11 +149,11 @@ function a(j::I, n::I, order, dt) where {I <: Integer}
end

function generalized_binomials(j, n, order, dt)
@. dt^order / order * ((n - j + 1)^order - (n - j)^order)
@. dt^order / order * ((n - j + 1)^order - (n - j )^order)
end

function v(cache::DelayPECECache{iip, T}, n) where {iip, T}
(; prob, mesh, dt, constant_lags, p) = cache
(; prob, mesh, dt, constant_lags, L, p) = cache
τ = constant_lags
if typeof(τ) <: Function
m = floor.(Int, τ.(mesh) / dt)
Expand All @@ -172,26 +168,32 @@ function v(cache::DelayPECECache{iip, T}, n) where {iip, T}
end
end
else
if τ >= (n - 1) * dt
if length(prob.u0) == 1
return [prob.h(p, (n - 1) * dt - τ)]
if isinteger(τ/dt)#y-τ fall in grid point
# case when δ == 0
if τ/dt < n
return cache.y[L + n - Int(τ/dt) + 1]
else
return prob.h(p, (n - 1) * dt - τ)
if length(prob.u0) == 1
return [prob.h(p, n * dt - τ)]
else
return prob.h(p, n * dt - τ)
end
end
else
m = floor(Int, τ / dt)
δ = m - τ / dt
# interpolate by the two nearest points
if m > 1
return δ * cache.y[n - m + 2] + (1 - δ) * cache.y[n - m + 1]
elseif m == 1
return δ * cache.yp[n + 1] + (1 - δ) * cache.y[n]
return δ * cache.y[n - m + L + 2] + (1 - δ) * cache.y[n - m + L + 1]
elseif m == 1 # h is larger than τ
return δ * cache.yp[n + 1] + (1 - δ) * cache.y[L + n + 1]
end
end
end
end

function solve_fdde_with_multiple_lags(FDDE::FDDEProblem, dt)
(; f, h, order, constant_lags, p, tspan) = FDDE
function solve_fdde_with_multiple_lags(cache::DelayPECECache)
(; f, h, order, constant_lags, p, tspan) = cache
τ = constant_lags[1]
mesh = collect(0:dt:tspan[2])
maxn = length(mesh)
Expand Down
4 changes: 2 additions & 2 deletions src/delay/product_integral.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,9 +75,9 @@ function SciMLBase.solve!(cache::DelayPIEXCache{iip, T}) where {iip, T}
prob.f(tmp, y[n], y_nm1_tau, p, tnm1)
g[n] = tmp
else
length(u0) == 1 ? (tmp = prob.f(y[n], y_nm1_tau[1], p, tnm1)) :
length(u0) == 1 ? (tmp = prob.f(y[n][1], y_nm1_tau[1], p, tnm1)) :
(tmp = prob.f(y[n], y_nm1_tau, p, tnm1))
g[n] = tmp
g[n] = (tmp isa AbstractVector) ? tmp : [tmp]
end
f_mem = zeros(T, l)
for j in 0:(n - 1)
Expand Down
14 changes: 7 additions & 7 deletions src/multitermsfode/explicit_pi.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,17 +67,17 @@ function SciMLBase.__init(

y = Vector{T}(undef, N + 1)
fy = similar(y)
zn = zeros(Float64, NNr + 1, orders_length)
zn = zeros(T, NNr + 1, orders_length)

nvett = collect(0:(NNr + 1))
bn = zeros(orders_length, NNr + 1)
bn = [Vector{T}(undef, NNr+1) for _ in 1:orders_length]
for i in 1:orders_length
nbeta = nvett .^ bet[i]
bn[i, :] = (nbeta[2:end] - nbeta[1:(end - 1)]) * dt^bet[i] / gamma(bet[i] + 1)
bn[i] = (nbeta[2:end] - nbeta[1:(end - 1)]) * dt^bet[i] / gamma(bet[i] + 1)
end
C = 0
for i in 1:(orders_length - 1)
C = C + lam_rat_i[i] * bn[i, 1]
C = C + lam_rat_i[i] * bn[i][1]
end

mesh = t0 .+ collect(0:N) * dt
Expand Down Expand Up @@ -175,7 +175,7 @@ function MTPX_multiterms_quadrato(cache::MTPIEXCache{T}, nxi, nxf, nyi, nyf) whe
funz_end = nyf + 1

for i in 1:orders_length
vett_coef = bn[i, coef_beg:coef_end]
vett_coef = bn[i][coef_beg:coef_end]
if i < orders_length
vett_funz = [permutedims(cache.y[funz_beg:funz_end]) zeros(1,
funz_end - funz_beg + 1)]
Expand Down Expand Up @@ -206,13 +206,13 @@ function MTPX_multiterms_triangolo(cache::MTPIEXCache{T}, nxi, nxf, t0) where {T
for i in 1:(orders_length - 1)
temp = zn[n + 1, i]
for j in j_beg:(n - 1)
temp = temp + bn[i, n - j] * cache.y[j + 1]
temp = temp + bn[i][n - j] * cache.y[j + 1]
end
Phi_n = Phi_n - lam_rat_i[i] * temp
end
temp = zn[n + 1, orders_length]
for j in j_beg:(n - 1)
temp = temp + bn[orders_length, n - j] * cache.fy[j + 1]
temp = temp + bn[orders_length][n - j] * cache.fy[j + 1]
end
Phi_n = Phi_n + temp / highest_order_parameter

Expand Down
4 changes: 2 additions & 2 deletions src/multitermsfode/implicit_pi_pece.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,9 +75,9 @@ function SciMLBase.__init(

y = Vector{T}(undef, N+1)
fy = similar(y)
zn_pred = zeros(NNr + 1, orders_length)
zn_pred = zeros(T, NNr + 1, orders_length)
if mu > 0
zn_corr = zeros(NNr + 1, orders_length)
zn_corr = zeros(T, NNr + 1, orders_length)
else
zn_corr = 0
end
Expand Down
4 changes: 1 addition & 3 deletions src/multitermsfode/implicit_pi_rectangle.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,6 @@ function SciMLBase.__init(prob::MultiTermsFODEProblem, alg::MTPIRect;
zn = zeros(NNr + 1, orders_length)

nvett = collect(0:(NNr + 1))
#bn = zeros(orders_length, NNr + 1)
bn = [Vector{T}(undef, NNr+1) for _ in 1:orders_length]
for i in 1:orders_length
nbeta = nvett .^ bet[i]
Expand All @@ -91,7 +90,7 @@ function SciMLBase.__init(prob::MultiTermsFODEProblem, alg::MTPIRect;
C = C + lam_rat_i[i] * bn[i][1]
end

mesh = collect(0:N) * dt
mesh = t0 .+ collect(0:N) * dt
y[1] = u0[1]
fy[1] = f(u0[1], p, t0)

Expand All @@ -103,7 +102,6 @@ end
function SciMLBase.solve!(cache::MTPIRectCache{T}) where {T}
(; prob, alg, mesh, u0, y, r, N, Qr) = cache
t0 = mesh[1]
tfinal = mesh[end]
MTPIRect_triangolo(cache, 1, r - 1, t0)

ff = zeros(1, 2^(Qr + 2))
Expand Down
10 changes: 5 additions & 5 deletions src/multitermsfode/implicit_pi_trapezoid.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ function SciMLBase.__init(prob::MultiTermsFODEProblem, alg::MTPITrap;

y = Vector{T}(undef, N + 1)
fy = similar(y)
zn = zeros(NNr + 1, orders_length)
zn = zeros(T, NNr + 1, orders_length)

nvett = collect(0:(NNr + 1))
an = [Vector{T}(undef, NNr + 1) for _ in 1:orders_length]#zeros(orders_length, NNr + 1)
Expand Down Expand Up @@ -140,10 +140,10 @@ function MTPITrap_disegna_blocchi(cache::MTPITrapCache{T}, L, ff, nx0, nu0, t0)
nyi::Int = nu0
nyf::Int = nu0 + L * r - 1
is = 1
s_nxi = zeros(N)
s_nxf = zeros(N)
s_nyi = zeros(N)
s_nyf = zeros(N)
s_nxi = Vector{T}(undef, N)
s_nxf = similar(s_nxi)
s_nyi = similar(s_nxi)
s_nyf = similar(s_nxi)
s_nxi[1] = nxi
s_nxf[1] = nxf
s_nyi[1] = nyi
Expand Down
17 changes: 17 additions & 0 deletions test/fdde.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,23 @@
atol = 1e-7)
end

@testset "Test problem in wang" begin
using SpecialFunctions

f(u, h, p, t) = 2/gamma(3-alpha)*t^(2-alpha) - 1/gamma(2-alpha)*t^(1-alpha) + 2*tau*t - tau^2 - tau - u + h
function h(p, t)
return 0.0
end
alpha = 0.5
tau = 0.5
u0 = 0.0
tspan = (0.0, 5.0)
analytical(u, h, p, t) = [t^2-t]
#fun = DDEFunction(f, analytic = analytical)
prob1 = FDDEProblem(f, alpha, u0, h, constant_lags = [tau], tspan)
sol = solve(prob1, DelayPECE(), dt=0.01)
end

@testset "Test DelayPECE method with single constant lag with variable order" begin
function h(p, t)
if t == 0
Expand Down

1 comment on commit f8bdf65

@github-actions
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Benchmark Results

Benchmark suite Current: f8bdf65 Previous: 99abd37 Ratio
FLMM/Trapezoid 15427749 ns 16208311 ns 0.95
FLMM/NewtonGregory 15470838 ns 16027141.5 ns 0.97
FLMM/BDF 13882942 ns 14540828 ns 0.95

This comment was automatically generated by workflow using github-action-benchmark.

Please sign in to comment.