Issue with solving a DDE in Julia

Question

I am trying to use Julia’s DifferentialEquations.jl package to solve a DDE system. I was able to solve my problem for two objects using the code below

clearconsole()
@time using DifferentialEquations
const two_bubble =
  let σ=0.0725,  ρ=998,  μ=1e-3,  c=1481,  pvtinf=0,  pinf0=1.01e5,  k=7/5,
      dist=10e-6, τ=dist/c,  L=[Inf dist;dist Inf], x=[0 dist],
      ps=100e3,    f=1e6,    R=[1e-6, 1e-6],
      SP=0.01, cycle=5;

    global function two_bubble!(du, u, h, p, t)



    A1=(1+(1-3*k)*u[2]/c)*((pinf0-pvtinf)/ρ+2*σ/(ρ*R[1]))*(R[1]/u[1])^(3*k)-2*σ/(ρ*u[1])-4*μ*u[2]/(ρ*u[1])-(1+u[2]/c)*(pinf0-pvtinf+ps*sin(2*pi*f*t-(2*pi*f/c)*x[1]))/ρ-ps*u[1]*cos(2*pi*f*t-(2*pi*f/c)*x[1])*2*pi*f/(ρ*c);
    A2=(1+(1-3*k)*u[4]/c)*((pinf0-pvtinf)/ρ+2*σ/(ρ*R[2]))*(R[2]/u[3])^(3*k)-2*σ/(ρ*u[3])-4*μ*u[4]/(ρ*u[3])-(1+u[4]/c)*(pinf0-pvtinf+ps*sin(2*pi*f*t-(2*pi*f/c)*x[2]))/ρ-ps*u[3]*cos(2*pi*f*t-(2*pi*f/c)*x[2])*2*pi*f/(ρ*c);

    du[1]= u[2]
    du[2]= (A1-(1.5*(1-u[2]/(3*c)))*u[2]^2-(2*h(p, t - τ; idxs = 3)*h(p, t - τ; idxs = 4)^2-h(p, t - τ, Val{1}; idxs = 4)*h(p, t - τ; idxs = 3)^2)/L[1, 2])/((1-u[2]/c)*u[1]+4*μ/(ρ*c))
    du[3]= u[4]
    du[4]= (A2-(1.5*(1-u[4]/(3*c)))*u[4]^2-(2*h(p, t - τ; idxs = 1)*h(p, t - τ; idxs = 2)^2-h(p, t - τ, Val{1}; idxs = 2)*h(p, t - τ; idxs = 1)^2)/L[1, 2])/((1-u[4]/c)*u[3]+4*μ/(ρ*c));

      nothing
    end

    global function h_two_bubble(p, t; idxs::Union{Nothing,Int} = nothing)
      t ≤ 0 || error("history function is only implemented for t ≤ 0")
R=[1e-6, 1e-6]
      if idxs === nothing
        [R[1],0,R[2],0]
      elseif idxs == 1
      R[1]
      elseif idxs == 2
      0
      elseif idxs == 3
      R[2]
      elseif idxs == 4
      0
      else
        error("delay differential equation consists of two components")
      end
    end

    global function h_two_bubble(p, t, ::Type{Val{1}};
                                    idxs::Union{Nothing,Int} = nothing)
      t ≤ 0 || error("history function is only implemented for t ≤ 0")

      if idxs === nothing
        [0,0,0,0]
      elseif idxs == 1 || idxs == 3
        0
      elseif idxs == 2 || idxs==4
        0
      else
        error("delay differential equation consists of two components")
      end
    end


    DDEProblem(two_bubble!, h_two_bubble, (0.0, cycle/f);
               constant_lags = [τ], neutral = true)

end



f=1e6;
SP=0.01

@time u=  solve(two_bubble,Tsit5(),saveat=SP/f,span=SP/f)
r=convert(Array,u)
t=u.t

When I am trying to make my code (MultiDelayedBuckle.jl) scalable to solve for arbitrary number of objects I am encountering an error I am not able to understand. Below you can see my code

clearconsole()


@time using DifferentialEquations
Sp=11; Sf=2.5e-6
SR0=0.01; kai=Sp/2; ks=Sf/(12*pi); RBB=1.02
      σ=0.0725;  ρ=998;  μ=1e-3;  c=1481;  pvtinf=0;  pinf0=1.01e5;  k=7/5;
      ps=1e3; f=1e6; R0=1.7e-6; Nb=2; d=100e-6; Mn=20e-6;
      SP=0.01; cycle=80;



R=zeros(1,Nb)
for i=1:Nb
  R[i]=R0
end
rbreak=R.*RBB; rb=R./sqrt(SR0/kai +1);


X=zeros(1,Nb)
Y=zeros(1,Nb)
Z=zeros(1,Nb)
L=zeros(Nb,Nb)
X[1]=d*rand(1)[1]
Y[1]=d*rand(1)[1]
Z[1]=d*rand(1)[1]
L[1,1]=Inf

k=1
while minimum(L)<Mn
global  i=Int64(2)
println("Attempt $k
")
while i<=Nb
X[i]=d*rand(1)[1]
Y[i]=d*rand(1)[1]
Z[i]=d*rand(1)[1]
for j=1:i
  if i==j
  L[i,j]=Inf
  else
  L[i,j]=sqrt(((X[i]-X[j])^2)+((Y[i]-Y[j])^2)+((Z[i]-Z[j])^2));
  L[j,i]=L[i,j];
  end
end
      global i=i+1
end
global k=k+1
end




τ=L/c;



p=[σ,ρ,μ,c,pvtinf,pinf0,k,ps,f,R,Nb,L,τ,SR0,kai,ks,rbreak,rb,X,τ]
################################################################################################################################################################################################################################################
function N_buckle!(du, u, h, p, t)  #with secondary delays
  σ,ρ,μ,c,pvtinf,pinf0,k,ps,f,R,Nb,L,τ,SR0,kai,ks,rbreak,rb,X,τ=p

  TT=zeros(1,Nb)
  P=zeros(1,Nb)
  SR=zeros(1,Nb)
  Psc=zeros(1,Nb)
  for i=1:Nb
   if 2*pi*f*t-(2*pi*f/c)*X[i]<0
        TT[i]=0
    else
        TT[i]=1
    end
    if u[2*i-1]>=rb[i] && u[2*i-1]<rbreak[i]
        SR[i]=kai*((u[2*i-1]/rb[i])^2-1)
    elseif u[2*i-1]>=rbreak[i]
        SR[i]=σ
    end

    P[i]=((pinf0+2*SR0)*(R[i]/u[2*i-1])^(3*k))*(1-3*k*u[2*i]/c)-(4*μ*u[2*i]/u[1])-(2*SR[i]/u[2*i-1])-(4*ks*u[2*i]/(u[2*i-1]^2))-pinf0-ps*TT[i]*sin(2*pi*f*t-(2*pi*f/c)*X[i]);


    Psc[i]=0
    for j=1:Nb
         if j~=i
         Pcs[i]=Pcs[i]+(2*h(p, t - τ[i,j]; idxs = 2*j-1)*h(p, t - τ[i,j]; idxs = 2*j)^2-h(p, t - τ[i,j], Val{1}; idxs = 2*j)*h(p, t - τ[i,j]; idxs = 2*j-1)^2)/(L[i,j])
         end
    end

du[2*i-1]=u[2*i]
du[2*i]=(P[2*i-1]/(ρ*u[2*i-1]))-(1.5*(u[2*i]^2)/u[2*i-1])-(Psc[i])/u[2*i-1]
end

  nothing
end
##################################################################################################
##################################################################################################
##################################################################################################
function h_N_buckle(p, t; idxs::Union{Nothing,Int} = nothing)
t ≤ 0 || error("history function is only implemented for t ≤ 0")
σ,ρ,μ,c,pvtinf,pinf0,k,ps,f,R,Nb,L,τ,SR0,kai,ks,rbreak,rb,X,τ=p

  u0=zeros(1,2*Nb)
  for i=1:2
  u0[2*i-1]=R[i]
  u0[2*i]=0
  end


if idxs==nothing
  u0
else
for i=1:2*Nb
 if idxs==2*i-1
  R[i]
elseif idxs==2*i
  0
 end
 end
end
end
##################################################################################################
##################################################################################################
##################################################################################################
function h_N_buckle(p, t, ::Type{Val{1}};
                                idxs::Union{Nothing,Int} = nothing)
t ≤ 0 || error("history function is only implemented for t ≤ 0")
σ,ρ,μ,c,pvtinf,pinf0,k,ps,f,R,Nb,L,τ,SR0,kai,ks,Rbreak,Rb,X,τ=p

  if idxs === nothing
    zeros(1,2*Nb)
     else
     for i=1:Nb
     0
     end
  end
end
##################################################################################################
##################################################################################################
##################################################################################################
MKMWD=DDEProblem(N_buckle!, h_N_buckle, (0.0, cycle/f),p;
               constant_lags = [τ], neutral = true)

SP=0.01;
tol=1e-4
alg=Tsit5()

u2=solve(MKMWD,alg,saveat=SP/f,reltol=tol,abstol=tol)

and here is the error I am getting:

LoadError: MethodError: no method matching isless(::Matrix{Float64}, ::Float64)
Closest candidates are:
  isless(!Matched::Union{StatsBase.PValue, StatsBase.TestStat}, ::AbstractFloat) at C:\Users\hhagh\.julia\packages\StatsBase\Lc3YW\src\statmodels.jl:504
  isless(!Matched::Union{StatsBase.PValue, StatsBase.TestStat}, ::Real) at C:\Users\hhagh\.julia\packages\StatsBase\Lc3YW\src\statmodels.jl:495
  isless(!Matched::Discontinuity, ::Number) at C:\Users\hhagh\.julia\packages\DelayDiffEq\otmBn\src\discontinuity_type.jl:21
  ...
in expression starting at D:\Julia\Working\MultiDelayedBuckle.jl:147
<(x::Matrix{Float64}, y::Float64) at operators.jl:279
initialize_tstops_d_discontinuities_propagated(#unused#::Type{Float64}, tstops::Tuple{}, d_discontinuities::Tuple{}, tspan::Tuple{Float64, Float64}, order_discontinuity_t0::Int64, alg_maximum_order::Int64, constant_lags::Vector{Matrix{Float64}}, neutral::Bool) at solve.jl:438
__init(prob::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, alg::MethodOfSteps{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, NLFunctional{Rational{Int64}, Rational{Int64}}, false}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}; saveat::Float64, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Float64, reltol::Float64, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, qoldinit::Rational{Int64}, fullnormalize::Bool, failfactor::Int64, beta1::Nothing, beta2::Nothing, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, discontinuity_interp_points::Int64, discontinuity_abstol::Float64, discontinuity_reltol::Int64, kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:default_set,), Tuple{Bool}}}) at solve.jl:187
(::SciMLBase.var"#__init##kw")(::NamedTuple{(:default_set, :saveat, :reltol, :abstol), Tuple{Bool, Float64, Float64, Float64}}, ::typeof(SciMLBase.__init), prob::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, alg::MethodOfSteps{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, NLFunctional{Rational{Int64}, Rational{Int64}}, false}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}) at solve.jl:67
__init at solve.jl:67 [inlined]
__solve(::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, ::MethodOfSteps{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, NLFunctional{Rational{Int64}, Rational{Int64}}, false}; kwargs::Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:default_set, :saveat, :reltol, :abstol), Tuple{Bool, Float64, Float64, Float64}}}) at solve.jl:4
(::SciMLBase.var"#__solve##kw")(::NamedTuple{(:default_set, :saveat, :reltol, :abstol), Tuple{Bool, Float64, Float64, Float64}}, ::typeof(SciMLBase.__solve), ::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, ::MethodOfSteps{CompositeAlgorithm{Tuple{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}}, AutoSwitch{Tsit5, Rosenbrock23{0, false, DefaultLinSolve, DataType}, Rational{Int64}, Int64}}, NLFunctional{Rational{Int64}, Rational{Int64}}, false}) at solve.jl:4
__solve(::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, ::Tsit5; default_set::Bool, kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:saveat, :reltol, :abstol), Tuple{Float64, Float64, Float64}}}) at default_solve.jl:7
__solve at default_solve.jl:3 [inlined]
#solve_call#56 at solve.jl:61 [inlined]
solve_call at solve.jl:48 [inlined]
#solve_up#58 at solve.jl:82 [inlined]
solve_up at solve.jl:75 [inlined]
#solve#57 at solve.jl:70 [inlined]
(::CommonSolve.var"#solve##kw")(::NamedTuple{(:saveat, :reltol, :abstol), Tuple{Float64, Float64, Float64}}, ::typeof(solve), prob::DDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, Vector{Matrix{Float64}}, Tuple{}, true, Vector{Any}, DDEFunction{true, typeof(N_buckle!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(h_N_buckle), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, args::Tsit5) at solve.jl:68
top-level scope at MultiDelayedBuckle.jl:147
eval at boot.jl:360 [inlined]
include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String) at loading.jl:1094

For reference I am using Julia 1.6.0.

Answer 1

This was a duplicate post of https://discourse.julialang.org/t/help-with-solving-a-dde/60205/10 which was narrowed down to constant_lag being an array of arrays, and thus the solution is simply:

MKMWD=DDEProblem(N_buckle!, h_N_buckle, (0.0, cycle/f),p;
               constant_lags = τ, neutral = true)