解决Mathematica中登革热DDE模型的NDSolveValue::rdelay延迟错误

我用Wolfram Mathematica构建了包含宿主-媒介交互、季节性效应和疫苗策略的登革热媒介动力学延迟微分方程（DDE）模型，但调用NDSolveValue求解时触发延迟时间相关错误：

NDSolveValue::rdelay: Delayed time {-1. τ} = -1. τ computed at t = 0 did not evaluate to a real number.

已为所有含延迟项的函数添加符合要求的初始历史条件，移除延迟项后系统可正常运行，说明问题仅出在延迟实现环节。现寻求具备Mathematica DDE处理经验的人士分析错误成因，提供延迟项改写或NDSolveValue参数调整的解决方案。

蚊子繁殖数与生命周期参数 #

Rm = 2.69  (*蚊子繁殖数*)
\[Omega] = 0.025 (*幼虫死亡率*)
\[CurlyEpsilon] = 0.1 (*成蚊死亡率*)
\[Alpha] = 1/19 (*幼虫平均发育时间*)
(* 通过设定Rm的值计算b *)
b = Rm / (\[Alpha]/(\[CurlyEpsilon] * (\[Alpha] + \[Omega])))
meq = 0.5 (*人均雌蚊数*)
K0 = meq / (\[Alpha]/\[CurlyEpsilon] (\[Alpha] (b - \
\[CurlyEpsilon])/(\[CurlyEpsilon] \[Omega]) - 
      1)) (*幼虫平均环境容纳量*)
\[CapitalDelta]K = 0.3 (*环境容纳量季节波动幅度，设为0.3以匹配登革热流行周期和振幅*)
POP = 1000000

K[t_] := K0*POP*(1 + \[CapitalDelta]K*Sin[2*Pi*t])

Plot[K[t], {t, 0, 10}]

疫苗参数与函数 #

\[Beta]mh = 1  (*蚊传人的单次叮咬传播概率，需匹配R0*)
\[Kappa]  = 0.6 (*单蚊日均叮咬率*)
TD = 8 (*疫苗诱导感染保护的平均持续时间*)
VEMinus = 0.87 (*血清阴性接种者的最大疫苗感染保护效力*)
VEPlus = 0.4 (*血清阳性接种者的最大疫苗感染保护效力*)
A = 2 (*接种年龄*)
P[0] = 0.45 (*原发感染症状比例*)
P[1] = 0.8 (*继发感染症状比例*)
P[2] = 0.12 (*三、四发感染症状比例*)
P[3] = 0.12 (*三、四发感染症状比例*)
P[4] = 0
P[5] = 0

S[v_, t_, a_, \[CapitalTheta]_] = 100000

(* 人群中血清型i的感染压力\[CapitalLambda]𝑖定义为： *)
\[CapitalLambda][i_, t_] := (\[Kappa] \[Beta]mh/M[t]) Y[t, i]

(* 保护衰减函数h(𝜏)：疫苗接种后保护呈指数衰减，平均持续时间为T *)
h[\[Tau]_] := Piecewise[
  {{Exp[-\[Tau]/TD], \[Tau] < 0.5},
   {Exp[-([Tau] - 0.5)/TD], 0.5 <= \[Tau] < 1},
   {Exp[-([Tau] - 1)/TD], \[Tau] >= 1}}
  ]

Plot[h[\[Tau]], {\[Tau], 0, 10}, PlotLegends -> "Expressions"]

(* 异源疫苗诱导的感染保护随时间衰减，相对风险函数𝑓(𝜏)（𝜏为接种后时间）定义为： *)
f[\[Tau]_, v_] := Piecewise[
  {{1, v == 0},
   {1 - VEMinus*h[\[Tau]], v == 1},
   {1 - VEPlus*h[\[Tau]], v == 2}}
  ]

Plot[{f[\[Tau], 0], f[\[Tau], 1], f[\[Tau], 2]}, {\[Tau], 0, 10}, 
 PlotLegends -> "Expressions"]

(* t时刻，血清型暴露史为theta、年龄a、疫苗状态v的人群中，血清型i的感染发生率 *)
InfIncidence = 
 cI[v_, t_, a_, i_, \[CapitalTheta]_] := \[CapitalLambda][i, t]* 
   f[a - A, v]* S[v, t, a, \[CapitalTheta]]

(* t时刻，血清型暴露史为theta、年龄a、疫苗状态v的人群中，血清型i感染导致的症状性疾病发生率 *)
SympIncidence = 
 cD[v_, t_, a_, i_, \[CapitalTheta]_] := 
  P[Length[\[CapitalTheta]] + (1 - 
       KroneckerDelta[0, v])]*\[CapitalLambda][i, t]*f[a - A, v]* 
   S[v, t, a, \[CapitalTheta]]

蚊子传染性参数与函数 #

\[Theta] = 2 (* 症状性感染的传染性较无症状感染的倍数 *)
\[Beta]hm = 1 (* 人传蚊的单次叮咬传播概率，暂设为1 *)
seroCombs = 
 Subsets[{1, 2, 3, 
   4}] (*生成所有血清型组合*) 
vacStats = {0, 1, 2} (*所有疫苗状态*)
TIP = 6 (* 内在潜伏期 *)
TInf = 4 (* 人群感染期 *)
maxAge = 120 (* 人群最大年龄 *)

innerFunction[i_, t_, \[Tau]_, a_, \[CapitalTheta]_, v_] := 
 If[Not[MemberQ[\[CapitalTheta], i]], 
  cI[v, t - \[Tau], a, i, \[CapitalTheta]] + (\[Theta] - 1)*
    cD[v, t - \[Tau], a, i, \[CapitalTheta]], 0]

(* 蚊子中血清型i的感染压力\[CapitalPsi]𝑖定义为：*)
\[CapitalPsi][i_, t_] := \[Kappa]*\[Beta]hm/POP*
  Integrate[
   Sum[
    innerFunction[i, t, \[Tau], a, \[CapitalTheta], v],
    {\[CapitalTheta], seroCombs},
    {v, vacStats}
    ],
   {a, 0, maxAge},
   {\[Tau], TIP, TIP + TInf}
   ]

微分方程系统 #

\[Omega] = 0.025 (*幼虫死亡率*)
\[Eta] = 1/10 (*外在潜伏期均值*)

eqLarvae = 
 D[L[t], t] == 
  b*M[t] - \[Alpha]*L[t] - \[Omega]*L[t] (1 + L[t]/K[t]) (*幼虫方程*)
eqAdults = 
 D[Am[t], t] == \[Alpha]*L[t] - 
   Sum[\[CapitalPsi][i, t]*Am[t], {i, 4}] - \[CurlyEpsilon]*
    Am[t] (*易感成蚊方程*)
eqInfected = 
 Flatten[Table[
   D[H[t, i, j], 
     t] == (KroneckerDelta[1, j]* \[CapitalPsi][i, t]*
       Am[t]) + (4 \[Eta]*(1 - KroneckerDelta[1, j])*
       H[t, i, j - 1]) - (4 \[Eta] + \[CurlyEpsilon])*H[t, i, j], {i, 
    4}, {j, 4}]] (* 感染但未具有传染性的蚊子方程 *)
eqInfectious = 
 Flatten[Table[
   D[Y[t, i], t] == 4 \[Eta]*H[t, i, 4] - \[CurlyEpsilon]*Y[t, i], {i,
     4}]] (* 传染性蚊子方程 *)

M[t_] := 
 Am[t] + Sum[H[t, i, j], {i, 4}, {j, 4}] + 
  Sum[Y[t, i], {i, 4}] (*成蚊总数，不含幼虫*)

histories = {
  L[t /; t <= 0] == 0,
  Am[t /; t <= 0] == 100000,
  H[t /; t <= 0, i, j] == 0 /. {i -> #, j -> #2} & @@@ 
   Tuples[Range[4], 2],
  Y[t /; t <= 0, i] == 1000 /. i -> # & /@ Range[4]
  }

求解代码 #

solution = NDSolveValue[
  Flatten[{eqLarvae, eqAdults, eqInfected, eqInfectious, histories}],
  Flatten[{L[t], Am[t], H[t, 1, 1], Y[t, 1]}],
  {t, 0, 200},
  Method -> {"EquationSimplification" -> "Residual"}
  ]

内容的提问来源于stack exchange，提问作者Lennert Saerens