# Help with solving a DDE

Hi!

I am trying to use Julia’s DifferentialEquations.jl package to solve DDE. 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

``````

I understand that the written code is not optimized at all but at this point I am only trying to get it to work and work on optimization after.

``````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)
``````

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

``````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:155
<(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:155
eval at boot.jl:360 [inlined]
``````

Any help would be greatly appreciated!

The error is pointing to line 155 and you shared a 147 line code, so this cannot be the same code you’re actually running.

My guess is what’s going on is that you needed `.<`.

P.S. Just do `zeros(Nb)` and simplify your code.

1 Like

You are right, here is the correct error message corresponding to correct line numbers

``````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]
``````

I tried this and it doesn’t work and now it returns the `history function is only implemented for t ≤ 0` that is defined in the code. I also tried `.≤` and it didn’t help either and resulted in the original error.

@ChrisRackauckas Do you have any other ideas that might be of help?

What’s the simplified code to copy paste and test?

This

``````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)
``````

Can you simplify it a bit first? `zeros(1,n)` should be `zeros(n)`. Etc. The tabbing is all off. That’s just really disrespectful of my time and I don’t have the time to do all of the cleanup before looking at the real issue. Does there really need to be all of those parameters in order to get the issue? You can’t delete a single parameter to recreate the issue? I don’t believe that.

2 Likes

oh I see, Yes that’s fair.

I managed to recreate the issue for the two object case and learned why it happens. just not sure how to address it
same error happens when I am addressing my delays in form of `τ[1,2], τ[2,1]` (see lines 36 and 38 where I am use them and line 5 where I am define them)

``````clearconsole()
@time using DifferentialEquations

σ=0.0725;  ρ=998;  μ=1e-3;  c=1481;  pvtinf=0;  pinf0=1.01e5;  k=7/5;
dist=10e-6;  L=[Inf dist;dist Inf]; τ=L/c; x=[0 dist];
ps=1e3; f=5e6; R=[2e-6, 2e-6]; Nb=2;
SR0=0.02; kai=0.9; ks=9e-9; RBB=1.01; rbreak=R.*RBB; rb=R./sqrt(SR0/kai +1);
cycle=80;
SP=0.01; tol=1e-5; alg=Tsit5()

################################################################################################################################################################################################################################################
function two_buckle!(du, u, h, p, t)
f=p[1];
TT=zeros(Nb)
P=zeros(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

SR=zeros(1,2)
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]);

end

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

nothing
end
########################################################################################################################################################################################################################################################################################
u0=zeros(2*Nb)
for i=1:2
u0[2*i-1]=R[i]
u0[2*i]=0
end
function h_two_buckle(p, t; idxs::Union{Nothing,Int} = nothing)
t ≤ 0 || error("history function is only implemented for t ≤ 0")

if idxs === nothing
u0
elseif mod(idxs,2)==1
R[Int64((idxs+1)/2)]
elseif mod(idxs,2)==0
0
end
end
########################################################################################################################################################################################################################################################################################
function h_two_buckle(p, t, ::Type{Val{1}};
idxs::Union{Nothing,Int} = nothing)
t ≤ 0 || error("history function is only implemented for t ≤ 0")

if idxs === nothing
zeros(Nb)
elseif mod(idxs,2)==1
0
elseif mod(idxs,2)==0
0
end
end
########################################################################################################################################################################################################################################################################################
p=[f]
MKMWD=DDEProblem(two_buckle!, h_two_buckle, (0.0, cycle/f),p;
constant_lags = [τ], neutral = true)

u=  DifferentialEquations.solve(MKMWD,alg,saveat=SP/f,reltol=tol,abstol=tol)

r=convert(Array,u)
t=u.t

``````

but when I use the delays like below (see lines 5, 36 and 38) it runs with no issues

``````clearconsole()
@time using DifferentialEquations

σ=0.0725;  ρ=998;  μ=1e-3;  c=1481;  pvtinf=0;  pinf0=1.01e5;  k=7/5;
dist=10e-6;  L=[Inf dist;dist Inf]; τ=dist/c; x=[0 dist];
ps=1e3; f=5e6; R=[2e-6, 2e-6]; Nb=2;
SR0=0.02; kai=0.9; ks=9e-9; RBB=1.01; rbreak=R.*RBB; rb=R./sqrt(SR0/kai +1);
cycle=80;
SP=0.01; tol=1e-5; alg=Tsit5()

################################################################################################################################################################################################################################################
function two_buckle!(du, u, h, p, t)
f=p[1];
TT=zeros(Nb)
P=zeros(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

SR=zeros(1,2)
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]);

end

du[1]= u[2]
du[2]= (P[1]/(ρ*u[1]))-(1.5*(u[2]^2)/u[1])-(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)/(u[1]*L[1, 2])
du[3]= u[4]
du[4]= (P[2]/(ρ*u[3]))-(1.5*(u[4]^2)/u[3])-(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)/(u[3]*L[2, 1])

nothing
end
########################################################################################################################################################################################################################################################################################
u0=zeros(2*Nb)
for i=1:2
u0[2*i-1]=R[i]
u0[2*i]=0
end
function h_two_buckle(p, t; idxs::Union{Nothing,Int} = nothing)
t ≤ 0 || error("history function is only implemented for t ≤ 0")

if idxs === nothing
u0
elseif mod(idxs,2)==1
R[Int64((idxs+1)/2)]
elseif mod(idxs,2)==0
0
end
end
########################################################################################################################################################################################################################################################################################
function h_two_buckle(p, t, ::Type{Val{1}};
idxs::Union{Nothing,Int} = nothing)
t ≤ 0 || error("history function is only implemented for t ≤ 0")

if idxs === nothing
zeros(Nb)
elseif mod(idxs,2)==1
0
elseif mod(idxs,2)==0
0
end
end
########################################################################################################################################################################################################################################################################################
p=[f]
MKMWD=DDEProblem(two_buckle!, h_two_buckle, (0.0, cycle/f),p;
constant_lags = [τ], neutral = true)

u=  DifferentialEquations.solve(MKMWD,alg,saveat=SP/f,reltol=tol,abstol=tol)

r=convert(Array,u)
t=u.t

``````

This is possible for the case of two objects in my equations since `τ[1,2]=τ[2,1]` but I cannot generalize it for the case of N objects. I also need to add the `τ[i,j]`s are constant and do not change throughout my solution.

Sorry it took a bit to get back to this. The answer is rather simple:

``````MKMWD=DDEProblem(two_buckle!, h_two_buckle, (0.0, cycle/f),p;
constant_lags = τ, neutral = true)
``````

You had `constant_lags = [τ]`. `constant_lags` is supposed to be an array of constant lag times. Your `τ` is already an array:

``````c=1481
L=[Inf dist;dist Inf]; τ=L/c
``````

so putting that array inside of another array confused it, and it was trying to act on `τ` as though it was a scalar. We could probably throw a better error for that, but hopefully that gets you working.

Thank you this solved my issue.