Skip to content

Commit f8bdf65

Browse files
ErikQQYErikQQY
committed
Several fix for FDDE solvers
1 parent 99abd37 commit f8bdf65

File tree

7 files changed

+104
-87
lines changed

7 files changed

+104
-87
lines changed

‎src/delay/pece.jl

Lines changed: 70 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
p
99
y
1010
yp
11+
L
1112

1213
dt
1314
kwargs
@@ -52,28 +53,18 @@ function SciMLBase.__init(prob::FDDEProblem, alg::DelayPECE; dt = 0.0, kwargs...
5253
t0 = tspan[1]
5354
tfinal = tspan[2]
5455
T = eltype(u0)
55-
mesh = collect(Float64, t0:dt:tfinal)
56-
size(mesh, 1)
57-
yp = collect(Float64, t0:dt:(tfinal + dt))
58-
N = length(t0:dt:(tfinal + dt))
59-
yp = _generate_similar_array(u0, N, u0)
60-
y = _generate_similar_array(u0, N - 1, u0)
61-
y[1] = (length(u0) == 1) ? _generate_similar_array(u0, 1, h(p, 0)) : u0#_generate_similar_array(u0, 1, h(p, 0))
62-
63-
return DelayPECECache{iip, T}(prob, alg, mesh, u0, order[1], τ, p, y, yp, dt, kwargs)
56+
N = Int(cld(tfinal - t0, dt))
57+
L = length(collect(-τ:dt:t0))
58+
mesh = collect(range(t0; stop = tfinal, length = N + 1))
59+
yp = _generate_similar_array(u0, N+1, u0)
60+
y = _generate_similar_array(u0, L+N+1, u0)
61+
y[L:-1:1] .= ifelse((length(u0) == 1), map(x->[h(p, x)], collect(-τ:dt:t0)), map(x->h(p, x), collect(-τ:dt:(t0))))
62+
y[L+1] = ifelse((length(u0) == 1), [u0], u0)#_generate_similar_array(u0, 1, h(p, 0))
63+
64+
return DelayPECECache{iip, T}(prob, alg, mesh, u0, order, τ, p, y, yp, L, dt, kwargs)
6465
end
6566

66-
@inline function _generate_similar_array(u0, N, src)
67-
if N ≠ 1
68-
if length(u0) == 1
69-
return [similar([src]) for i in 1:N]
70-
else
71-
return [similar(src) for i in 1:N]
72-
end
73-
else
74-
return [src]
75-
end
76-
end
67+
@inline _generate_similar_array(u0, N, src) = ifelse(length(u0) == 1, [[zero(u0)] for _ in 1:N], [zero(src) for _ in 1:N])
7768

7869
@inline function OrderWrapper(order, t)
7970
if order isa Function
@@ -84,65 +75,70 @@ end
8475
end
8576

8677
function SciMLBase.solve!(cache::DelayPECECache{iip, T}) where {iip, T}
87-
(; prob, alg, mesh, u0, order, p, dt) = cache
88-
maxn = length(mesh)
78+
if length(cache.constant_lags) == 1
79+
return solve_fdde_with_single_lag(cache)
80+
else
81+
return solve_fdde_with_multiple_lags(cache)
82+
end
83+
end
84+
85+
function solve_fdde_with_single_lag(cache::DelayPECECache{iip, T}) where {iip, T}
86+
(; prob, alg, mesh, u0, order, L, p, dt) = cache
87+
N = length(mesh)
8988
l = length(u0)
90-
initial = (length(u0) == 1) ? _generate_similar_array(u0, 1, prob.h(p, 0)) : u0
91-
for n in 1:(maxn - 1)
92-
order = OrderWrapper(order, mesh[n + 1])
93-
cache.yp[n + 1] = zeros(T, l)
94-
for j in 1:n
89+
initial = ifelse(l == 1, [u0], u0)
90+
tmp = zeros(T, l)
91+
for i in 0:N-2
92+
order = OrderWrapper(order, mesh[i + 1])
93+
# compute the yp part
94+
for j in 0:i
9595
if iip
96-
tmp = zeros(length(cache.yp[1]))
97-
prob.f(tmp, cache.y[j], v(cache, j), p, mesh[j])
98-
cache.yp[n + 1] = cache.yp[n + 1] +
99-
generalized_binomials(j - 1, n - 1, order, dt) * tmp
96+
prob.f(tmp, cache.y[L+j+1], v(cache, j), p, mesh[j + 1])
97+
cache.yp[i + 1] = cache.yp[i + 1] +
98+
generalized_binomials(j, i, order, dt) * tmp
10099
else
101100
length(u0) == 1 ?
102-
(tmp = prob.f(cache.y[j][1], v(cache, j)[1], p, mesh[j])) :
103-
(tmp = prob.f(cache.y[j], v(cache, j), p, mesh[j]))
104-
@. cache.yp[n + 1] = cache.yp[n + 1] +
105-
generalized_binomials(j - 1, n - 1, order, dt) * tmp
101+
(tmp = prob.f(cache.y[L+j+1][1], v(cache, j)[1], p, mesh[j + 1])) :
102+
(tmp = prob.f(cache.y[L+j+1], v(cache, j), p, mesh[j+1]))
103+
@. cache.yp[i + 1] = cache.yp[i + 1] +
104+
generalized_binomials(j, i, order, dt) * tmp
106105
end
107106
end
108-
cache.yp[n + 1] = cache.yp[n + 1] / gamma(order) + initial
107+
cache.yp[i + 1] = cache.yp[i + 1] / gamma(order) + initial
109108

110-
cache.y[n + 1] = zeros(T, l)
111-
112-
@fastmath @inbounds @simd for j in 1:n
109+
# compute the a part
110+
for j in 0:i
113111
if iip
114-
tmp = zeros(T, length(cache.y[1]))
115-
prob.f(tmp, cache.y[j], v(cache, j), p, mesh[j])
116-
cache.y[n + 1] = cache.y[n + 1] + a(j - 1, n - 1, order, dt) * tmp
112+
prob.f(tmp, cache.y[L + j + 1], v(cache, j), p, mesh[j+1])
113+
cache.y[L + i + 2] += a(j, i, order, dt) * tmp
117114
else
118115
length(u0) == 1 ?
119-
(tmp = prob.f(cache.y[j][1], v(cache, j)[1], p, mesh[j])) :
120-
(tmp = prob.f(cache.y[j], v(cache, j), p, mesh[j]))
121-
@. cache.y[n + 1] = cache.y[n + 1] + a(j - 1, n - 1, order, dt) * tmp
116+
(tmp = prob.f(cache.y[L + j + 1][1], v(cache, j)[1], p, mesh[j+1])) :
117+
(tmp = prob.f(cache.y[L + j + 1], v(cache, j), p, mesh[j]))
118+
@. cache.y[L + i + 2] += a(j, i, order, dt) * tmp
122119
end
123120
end
124121

125122
if iip
126-
tmp = zeros(T, l)
127-
prob.f(tmp, cache.yp[n + 1], v(cache, n + 1), p, mesh[n + 1])
128-
cache.y[n + 1] = cache.y[n + 1] / gamma(order) +
129-
dt^order * tmp / gamma(order + 2) +
130-
initial
123+
prob.f(tmp, cache.yp[i + 1], v(cache, i + 1), p, mesh[i + 2])
124+
cache.y[L + i + 2] = cache.y[L + i + 2] / gamma(order) +
125+
dt^order * tmp / gamma(order + 2) + initial
131126
else
132127
length(u0) == 1 ?
133-
(tmp = prob.f(cache.yp[n + 1][1], v(cache, n + 1)[1], p, mesh[n + 1])) :
134-
(tmp = prob.f(cache.yp[n + 1], v(cache, n + 1), p, mesh[n + 1]))
135-
@. cache.y[n + 1] = cache.y[n + 1] / gamma(order) +
136-
dt^order * tmp / gamma(order + 2) +
137-
initial
128+
(tmp = prob.f(cache.yp[i + 1][1], v(cache, i + 1)[1], p, mesh[i + 2])) :
129+
(tmp = prob.f(cache.yp[i + 1], v(cache, i + 1), p, mesh[i + 2]))
130+
@. cache.y[L + i + 2] = cache.y[L + i + 2] / gamma(order) +
131+
dt^order * tmp / gamma(order + 2) + initial
138132
end
139133
end
140134

141-
return DiffEqBase.build_solution(prob, alg, mesh, cache.y)
135+
u = cache.y[L+1:end]
136+
137+
return DiffEqBase.build_solution(prob, alg, mesh, u)
142138
end
143139

144140
function a(j::I, n::I, order, dt) where {I <: Integer}
145-
if j == n + 1
141+
if j == n
146142
result = 1
147143
elseif j == 0
148144
result = n^(order + 1) - (n - order) * (n + 1)^order
@@ -153,11 +149,11 @@ function a(j::I, n::I, order, dt) where {I <: Integer}
153149
end
154150

155151
function generalized_binomials(j, n, order, dt)
156-
@. dt^order / order * ((n - j + 1)^order - (n - j)^order)
152+
@. dt^order / order * ((n - j + 1)^order - (n - j )^order)
157153
end
158154

159155
function v(cache::DelayPECECache{iip, T}, n) where {iip, T}
160-
(; prob, mesh, dt, constant_lags, p) = cache
156+
(; prob, mesh, dt, constant_lags, L, p) = cache
161157
τ = constant_lags
162158
if typeof(τ) <: Function
163159
m = floor.(Int, τ.(mesh) / dt)
@@ -172,26 +168,32 @@ function v(cache::DelayPECECache{iip, T}, n) where {iip, T}
172168
end
173169
end
174170
else
175-
if τ >= (n - 1) * dt
176-
if length(prob.u0) == 1
177-
return [prob.h(p, (n - 1) * dt - τ)]
171+
if isinteger(τ/dt)#y-τ fall in grid point
172+
# case when δ == 0
173+
if τ/dt < n
174+
return cache.y[L + n - Int(τ/dt) + 1]
178175
else
179-
return prob.h(p, (n - 1) * dt - τ)
176+
if length(prob.u0) == 1
177+
return [prob.h(p, n * dt - τ)]
178+
else
179+
return prob.h(p, n * dt - τ)
180+
end
180181
end
181182
else
182183
m = floor(Int, τ / dt)
183184
δ = m - τ / dt
185+
# interpolate by the two nearest points
184186
if m > 1
185-
return δ * cache.y[n - m + 2] + (1 - δ) * cache.y[n - m + 1]
186-
elseif m == 1
187-
return δ * cache.yp[n + 1] + (1 - δ) * cache.y[n]
187+
return δ * cache.y[n - m + L + 2] + (1 - δ) * cache.y[n - m + L + 1]
188+
elseif m == 1 # h is larger than τ
189+
return δ * cache.yp[n + 1] + (1 - δ) * cache.y[L + n + 1]
188190
end
189191
end
190192
end
191193
end
192194

193-
function solve_fdde_with_multiple_lags(FDDE::FDDEProblem, dt)
194-
(; f, h, order, constant_lags, p, tspan) = FDDE
195+
function solve_fdde_with_multiple_lags(cache::DelayPECECache)
196+
(; f, h, order, constant_lags, p, tspan) = cache
195197
τ = constant_lags[1]
196198
mesh = collect(0:dt:tspan[2])
197199
maxn = length(mesh)

‎src/delay/product_integral.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -75,9 +75,9 @@ function SciMLBase.solve!(cache::DelayPIEXCache{iip, T}) where {iip, T}
7575
prob.f(tmp, y[n], y_nm1_tau, p, tnm1)
7676
g[n] = tmp
7777
else
78-
length(u0) == 1 ? (tmp = prob.f(y[n], y_nm1_tau[1], p, tnm1)) :
78+
length(u0) == 1 ? (tmp = prob.f(y[n][1], y_nm1_tau[1], p, tnm1)) :
7979
(tmp = prob.f(y[n], y_nm1_tau, p, tnm1))
80-
g[n] = tmp
80+
g[n] = (tmp isa AbstractVector) ? tmp : [tmp]
8181
end
8282
f_mem = zeros(T, l)
8383
for j in 0:(n - 1)

‎src/multitermsfode/explicit_pi.jl

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -67,17 +67,17 @@ function SciMLBase.__init(
6767

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

7272
nvett = collect(0:(NNr + 1))
73-
bn = zeros(orders_length, NNr + 1)
73+
bn = [Vector{T}(undef, NNr+1) for _ in 1:orders_length]
7474
for i in 1:orders_length
7575
nbeta = nvett .^ bet[i]
76-
bn[i, :] = (nbeta[2:end] - nbeta[1:(end - 1)]) * dt^bet[i] / gamma(bet[i] + 1)
76+
bn[i] = (nbeta[2:end] - nbeta[1:(end - 1)]) * dt^bet[i] / gamma(bet[i] + 1)
7777
end
7878
C = 0
7979
for i in 1:(orders_length - 1)
80-
C = C + lam_rat_i[i] * bn[i, 1]
80+
C = C + lam_rat_i[i] * bn[i][1]
8181
end
8282

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

177177
for i in 1:orders_length
178-
vett_coef = bn[i, coef_beg:coef_end]
178+
vett_coef = bn[i][coef_beg:coef_end]
179179
if i < orders_length
180180
vett_funz = [permutedims(cache.y[funz_beg:funz_end]) zeros(1,
181181
funz_end - funz_beg + 1)]
@@ -206,13 +206,13 @@ function MTPX_multiterms_triangolo(cache::MTPIEXCache{T}, nxi, nxf, t0) where {T
206206
for i in 1:(orders_length - 1)
207207
temp = zn[n + 1, i]
208208
for j in j_beg:(n - 1)
209-
temp = temp + bn[i, n - j] * cache.y[j + 1]
209+
temp = temp + bn[i][n - j] * cache.y[j + 1]
210210
end
211211
Phi_n = Phi_n - lam_rat_i[i] * temp
212212
end
213213
temp = zn[n + 1, orders_length]
214214
for j in j_beg:(n - 1)
215-
temp = temp + bn[orders_length, n - j] * cache.fy[j + 1]
215+
temp = temp + bn[orders_length][n - j] * cache.fy[j + 1]
216216
end
217217
Phi_n = Phi_n + temp / highest_order_parameter
218218

‎src/multitermsfode/implicit_pi_pece.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -75,9 +75,9 @@ function SciMLBase.__init(
7575

7676
y = Vector{T}(undef, N+1)
7777
fy = similar(y)
78-
zn_pred = zeros(NNr + 1, orders_length)
78+
zn_pred = zeros(T, NNr + 1, orders_length)
7979
if mu > 0
80-
zn_corr = zeros(NNr + 1, orders_length)
80+
zn_corr = zeros(T, NNr + 1, orders_length)
8181
else
8282
zn_corr = 0
8383
end

‎src/multitermsfode/implicit_pi_rectangle.jl

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -80,7 +80,6 @@ function SciMLBase.__init(prob::MultiTermsFODEProblem, alg::MTPIRect;
8080
zn = zeros(NNr + 1, orders_length)
8181

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

94-
mesh = collect(0:N) * dt
93+
mesh = t0 .+ collect(0:N) * dt
9594
y[1] = u0[1]
9695
fy[1] = f(u0[1], p, t0)
9796

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

109107
ff = zeros(1, 2^(Qr + 2))

‎src/multitermsfode/implicit_pi_trapezoid.jl

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -78,7 +78,7 @@ function SciMLBase.__init(prob::MultiTermsFODEProblem, alg::MTPITrap;
7878

7979
y = Vector{T}(undef, N + 1)
8080
fy = similar(y)
81-
zn = zeros(NNr + 1, orders_length)
81+
zn = zeros(T, NNr + 1, orders_length)
8282

8383
nvett = collect(0:(NNr + 1))
8484
an = [Vector{T}(undef, NNr + 1) for _ in 1:orders_length]#zeros(orders_length, NNr + 1)
@@ -140,10 +140,10 @@ function MTPITrap_disegna_blocchi(cache::MTPITrapCache{T}, L, ff, nx0, nu0, t0)
140140
nyi::Int = nu0
141141
nyf::Int = nu0 + L * r - 1
142142
is = 1
143-
s_nxi = zeros(N)
144-
s_nxf = zeros(N)
145-
s_nyi = zeros(N)
146-
s_nyf = zeros(N)
143+
s_nxi = Vector{T}(undef, N)
144+
s_nxf = similar(s_nxi)
145+
s_nyi = similar(s_nxi)
146+
s_nyf = similar(s_nxi)
147147
s_nxi[1] = nxi
148148
s_nxf[1] = nxf
149149
s_nyi[1] = nyi

‎test/fdde.jl

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,23 @@
2727
atol = 1e-7)
2828
end
2929

30+
@testset "Test problem in wang" begin
31+
using SpecialFunctions
32+
33+
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
34+
function h(p, t)
35+
return 0.0
36+
end
37+
alpha = 0.5
38+
tau = 0.5
39+
u0 = 0.0
40+
tspan = (0.0, 5.0)
41+
analytical(u, h, p, t) = [t^2-t]
42+
#fun = DDEFunction(f, analytic = analytical)
43+
prob1 = FDDEProblem(f, alpha, u0, h, constant_lags = [tau], tspan)
44+
sol = solve(prob1, DelayPECE(), dt=0.01)
45+
end
46+
3047
@testset "Test DelayPECE method with single constant lag with variable order" begin
3148
function h(p, t)
3249
if t == 0

0 commit comments

Comments
 (0)