Delay Differential Equations Networks with different delay per edge

Specific Domains Modelling & Simulations
,
6

solve(prob, MethodOfSteps(BS3()), saveat=0.01, reltol=0, abstol=1e-5). rtol and atol are not options, and you should almost never force_dtmin. I should really just remove that argument. If you’re hitting it, then your definition is probably wrong: the warning even says to think of dtmin as divergence for a reason described in PSA: How to help yourself debug differential equation solving issues .

And so the first thing is to notice is that these are not the same DDE being solved. In fact, the Julia one isn’t even a DDE: it’s not even a well-defined equation!

h(p, t; idxs=nothing) = typeof(idxs) <: Number ? 1.0 : rand(Uniform(0, 2π), n)

The past does not have a value: every single time the past is called the value changes. So what’s happening is that the starting iterations are unable to converge because there is no sense of convergence: there are no values!

But it’s even worse then that: you’re saying that "the u(t) for t<0 is a random value at every t, but u(t)[1] = 1.0 for any past value. So that’s even worse, it’s not even consistent :scream:. And u(t)[1] = 1.0 does not match what is in the Python code which uses a uniform random value for the past.

The fixed code is:

using DelayDiffEq
using Distributions
using Random

function f!(dθ, θ, h, p, t)
    ω, A = p
    n = length(θ)
    lags = reshape(lag, n,n)
    @inbounds for j in 1:n
        coupling = 0.0
        @inbounds for i in 1:n
            coupling += A[i,j]*sin(h(p, t-lags[i,j]; idxs=i) - θ[j])
        end
        dθ[j] = ω[j] + coupling
    end
    nothing
end

n = 10
Random.seed!(1)
ω = rand(n)
A = rand(n,n)
const lag = rand(n*n)
θ₀ = rand(Uniform(0, 2π), n)
p = (ω, A)
const past = rand(Uniform(0, 2π), n)
h(p, t; idxs=nothing) = typeof(idxs) <: Number ? past[idxs] : past

prob = DDEProblem(f!, θ₀, h, (0.0, 1.0), p, constant_lags=lag)

using Profile
@profile solve(prob, MethodOfSteps(BS3()), saveat=0.01, reltol=0, abstol=1e-5)
Juno.profiler()

After fixing that it’s still slow though. The issue seems to be that Julia’s inference is bugging in f!. This could be https://github.com/JuliaLang/julia/issues/35800 or https://github.com/SciML/OrdinaryDiffEq.jl/pull/1473

3 Likes