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.