Different lags for DDEProblem

Specific Domains Optimization (Mathematical)
, ,
1

Hi! I try to use delay differential equations from DifferentialEquations.jl.

I look at docs and make system with different delays, it works but looks not fine:

 hist1 = h(p, t)[1]
 hist2 = h(p, t - 2)[2]
 hist3 = h(p, t - 10)[3]

I tried to make something with constant_lags=lags, but I can’t understand how it works. How i should use DDE with different delays correctly?

Example:

using DifferentialEquations, Plots

h(p, t) = zeros(3)

function pkiv(du, u, h, p, t)
    hist1 = h(p, t)[1]
    hist2 = h(p, t - 2)[2]
    hist3 = h(p, t - 10)[3]
    
    k₁, k₂, k₃ = p
    A, B, C   = u
    du[1] =  - k₁ * A * B  + k₃*hist3
    du[2] =    k₁ * A * B - k₂*hist2
    du[3] =    k₂ * hist2 - k₃*hist3
   
end

u0 = [10.0;1.0;0]
p  = [0.1, 0.15, 0.1]
tspan = (0.0, 60.0)
prob = DDEProblem(pkiv, u0, h, tspan, p; )
alg = MethodOfSteps(Tsit5())
sol = solve(prob,alg)
plot(sol; label = ["C" "P"])
2

Hi, I think your misunderstanding is related to where h comes from. It is generated by the solver.

In your method pkiv you should simply call h(t)

https://diffeq.sciml.ai/stable/types/dde_types/#SciMLBase.DDEProblem

3

You can use, e.g.,

using DifferentialEquations, Plots

h(p, t) = zeros(3)

function pkiv(du, u, h, p, t)
    hist2 = h(p, t - 2)[2]
    hist3 = h(p, t - 10)[3]

    k₁, k₂, k₃ = p
    A, B, C  = u
    du[1] =  - k₁ * A * B  + k₃ * hist3
    du[2] =    k₁ * A * B - k₂ * hist2
    du[3] =    k₂ * hist2 - k₃ * hist3

    nothing
end

u0 = [10.0, 1.0, 0]
p = [0.1, 0.15, 0.1]
tspan = (0.0, 60.0)
prob = DDEProblem(pkiv, u0, h, tspan, p; constant_lags=(2, 10))
alg = MethodOfSteps(Tsit5())
sol = solve(prob, alg)
plot(sol)

I removed hist1 since it was not used. Moreover, h(p, t) will be equal to u as that is the “history” at the current time with the given parameters.

constant_lags should be a collection (i.e., tuple, array, etc.) of the constant delays (if you want to specify them). The documentation also contains a list of other problem-specific options.

By the way, if you only use the DDE solvers, you can reduce loading and compilation times by only loading DelayDiffEq instead of DifferentialEquations.

In this case you can also speed up the evaluation of the history function by extracting the desired indices directly:

using DifferentialEquations, Plots

function h_idxs(p, t; idxs::Union{Int,Nothing}=nothing)
    if idxs === nothing
        zeros(3)
    else
        0.0
    end
end

function pkiv_idxs(du, u, h, p, t)
    hist2 = h(p, t - 2; idxs=2)
    hist3 = h(p, t - 10; idxs=3)

    k₁, k₂, k₃ = p
    A, B, C  = u
    du[1] =  - k₁ * A * B  + k₃ * hist3
    du[2] =    k₁ * A * B - k₂ * hist2
    du[3] =    k₂ * hist2 - k₃ * hist3

    nothing
end

u0 = [10.0, 1.0, 0]
p = [0.1, 0.15, 0.1]
tspan = (0.0, 60.0)
prob_idxs = DDEProblem(pkiv_idxs, u0, h_idxs, tspan, p; constant_lags=(2, 10))
alg = MethodOfSteps(Tsit5())
sol = solve(prob_idxs, alg)
plot(sol)

A quick timing with @time confirms that this reduces allocations (but does not seem to have a significant impact on computation time in this example):

julia> @time solve(prob, alg);
  0.000178 seconds (2.72 k allocations: 232.328 KiB)

julia> @time solve(prob_idxs, alg);
  0.000176 seconds (1.55 k allocations: 141.078 KiB)
2 Likes