在 Julia 中解决 DDE 的问题

如何解决在 Julia 中解决 DDE 的问题

我正在尝试使用 Julia 的 DifferentialEquations.jl 包来解决 DDE 系统。我能够使用下面的代码解决两个对象的问题

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

当我尝试使我的代码 (MultiDelayedBuckle.jl) 可扩展以解决任意数量的对象时，我遇到了一个我无法理解的错误。下面你可以看到我的代码

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,t)  #with secondary delays
  σ,τ=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,j]; idxs = 2*j)^2-h(p,j],Val{1}; idxs = 2*j)*h(p,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,Int} = nothing)
t ≤ 0 || error("history function is only implemented for t ≤ 0")
σ,τ=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,::Type{Val{1}};
                                idxs::Union{nothing,Rbreak,Rb,τ=p

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

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

u2=solve(MKMWD,alg,reltol=tol,abstol=tol)

这是我得到的错误：

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,::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,Vector{Matrix{Float64}},Tuple{},true,Vector{Any},Ddefunction{true,typeof(N_buckle!),Linearalgebra.UniformScaling{Bool},nothing,nothing},typeof(h_N_buckle),Base.Iterators.Pairs{Union{},Union{},NamedTuple{(),Tuple{}}},SciMLBase.StandardDDEProblem},alg::MethodofSteps{Compositealgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,DataType},Rational{Int64},Int64}},NLFunctional{Rational{Int64},Rational{Int64}},false},timeseries_init::Tuple{},ts_init::Tuple{},ks_init::Tuple{}; saveat::Float64,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}},::typeof(SciMLBase.__init),prob::DDEProblem{Matrix{Float64},ks_init::Tuple{}) at solve.jl:67
__init at solve.jl:67 [inlined]
__solve(::DDEProblem{Matrix{Float64},::MethodofSteps{Compositealgorithm{Tuple{Tsit5,false}; kwargs::Base.Iterators.Pairs{Symbol,Real,NTuple{4,Symbol},Float64}}}) at solve.jl:4
(::SciMLBase.var"#__solve##kw")(::NamedTuple{(:default_set,::typeof(SciMLBase.__solve),::DDEProblem{Matrix{Float64},false}) at solve.jl:4
__solve(::DDEProblem{Matrix{Float64},::Tsit5; default_set::Bool,Tuple{Symbol,Symbol,NamedTuple{(:saveat,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,::typeof(solve),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

作为参考，我使用的是 Julia 1.6.0。

解决方法

这是 https://discourse.julialang.org/t/help-with-solving-a-dde/60205/10 的重复帖子，它被缩小为 constant_lag 是一个数组数组，因此解决方案很简单：

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