I have a system of DDEs that model host-parasite dynamics in a seasonal environment and I am having some issues because (I think) the equations are so stiff. Basically, when the number of hosts
S is too high (> about 500) I have to lower the solver tolerances so much that it becomes prohibitively slow (I need
S up to about a 4e6). I believe the low tolerances are required because
u decreases to zero so rapidly when
S is high that the history function (
W_past) just returns zero, which makes
matU = 0.0 and therefore
u stays at 0.0 (i.e., no infective larvae develop) in the model
Re_MWE. It may be the case that I just have to use low tolerances and run the simulation for hours and hours, but if anyone has any suggestions to make this work without having to lower the tolerances too much I would be very grateful. I tried several different stiff solvers (Rosenbrock23, Rodas3,4,4P&5, etc.) but without any great success. Maybe some way of adjusting the step size or addressing the discontinuity is possible (I’ve played around with some solver options but I’m not super clear on exactly what each one does or how best to use them yet)? Any other general performance tips are also appreciated.
Here is a stripped down MWE:
using QuadGK, Roots, DelayDiffEq, Plots #parameters const Temp0K = 15+273.15; const k = 8.62 * 10^(-5); const S = 100.0 #Host density const mu_W_a = 0.0186; const mu_W_E = 0.87; const mu_W_EH = 5.0*mu_W_E; const mu_W_TH = 289.65; const Phi_a = 0.024; const Phi_E = 0.46; const mu_U_a = 0.0001; const mu_U_E = 0.87; const mu_U_EH = 5.0*mu_U_E; const mu_U_TH = 289.65; const tauU_a = 25.1; const tauU_E = 0.87; const tauU_EL = 5.0*tauU_E; const tauU_TL = 282.0218; #Functions for temperature-dependent parameters function g_U(T) tmp = 1 / max(TempKstop, T) 1 / (tauU_a * exp((tauU_E/k) * (tmp - (1 / Temp0K))) * (1 + exp((tauU_EL / k) * (tmp - (1 / tauU_TL))))) end function mu_W(T) tmp = 1 / max(TempKstop, T) mu_W_a * exp((-mu_W_E/k)*((tmp)-(1/Temp0K))) * (1 + exp((mu_W_EH/k)*((-tmp)+(1/mu_W_TH)))); #L1 mort end function mu_U(T) tmp = 1 / max(TempKstop, T) mu_U = mu_U_a * exp((-mu_U_E/k)*((tmp)-(1/Temp0K))) * (1 + exp((mu_U_EH/k)*((-tmp)+(1/mu_U_TH)))); #U mort end function Phi(T) tmp = 1 / max(TempKstop, T) Phi = Phi_a * exp((-Phi_E/k)*((tmp)-(1/Temp0K))); #W penetration into slugs end #Transmission model function Re_MWE(du,u,h,p,t) #State variables W = u; U = u; I = u; DU1 = u; tauU = u; # History calculations ucache = h(p, t - tauU) W_past = ucache U_past = ucache # Compute temperature T = Temp(t) T_past = Temp(t - tauU) # Preliminary calculations gU_ratio = g_U(T) / g_U(T_past) matW = Phi(T) * S * W #Maturation rate of W matU = Phi(T_past) * S * W_past * DU1 * gU_ratio #Model du = - mu_W(T) * W - matW #W) Free-living larvae du = matW - U * mu_U(T) - matU #U) Uninfective larvae du = matU #I) Infective larvae du = DU1 * (mu_U(T_past) * U_past / S * gU_ratio - mu_U(T) * U / S) #Survival probability through delay period du = 1 - gU_ratio #tauU) Change in delay duration nothing end #Setup simualation #Initial conditions and history const InitW = 1.0; const InitU = 0.0; const InitI = 0.0; function DySimulate1(alg, tstart, tfin, tol) gint(u) = quadgk(g_U ∘ Temp, tstart - u, tstart) initTauU = [fzero(u -> gint(u) - p, [0.001, 5000]) for p in [1.0]] #Initial lag for each day (where guadgk = 1.0) mUint = exp(-(quadgk(mu_U ∘ Temp, tstart-initTauU, tstart))) DUinit = mUint DUhist = DUinit u0 = [InitW, InitU, InitI, DUinit,initTauU] h(p,t) = [0.0, 0.0, 0.0, DUhist, initTauU] tspan = (tstart, tfin+tstart) # Declare dependent lags probDynSea = DDEProblem(Re_MWE, u0, h, tspan; dependent_lags = ((u, p, t) -> u,)) solDynSimpSea = solve(probDynSea,alg; maxiters = 1e7, isoutofdomain=(u,p,t) -> any(x -> x < 0.0, u), abstol = tol, reltol = tol, save_idxs = 3) end alg = MethodOfSteps(Rosenbrock23()) # alg = MethodOfSteps(Rosenbrock32()) # alg = MethodOfSteps(Rodas4P()) # alg = MethodOfSteps(Rodas4()) # alg = MethodOfSteps(Rodas5()) # alg = MethodOfSteps(TRBDF2()) # alg = MethodOfSteps(KenCarp4()) #Temperature function const c_K = 263.0; #Mean annual temperature (Kelvins) const d_K = 0.66478 * c_K + -153.14 const t0 = -1.6962 * c_K +551.43 const TempKstop = 273.15-0.0 Temp(t) = c_K + d_K * sin((t-t0)*2*pi/365) MWEanswer = DySimulate1(alg,1.0,365.0*2, 1e-4) p = plot(MWEanswer) # scatter!(p,MWEanswer.t, MWEanswer.u) # p
When both the solver tolerances are
tol = 1e-4 and
S = 100.0 then I get this (as expected):
tol = 1e-4 and
S = 700.0 I get this:
Once I lower the tolerances to
tol = 1e-6 with
S = 700.0 it works again, but is much much slower and requires larger maxiters.
However, with my full model I end up needing to go to at least
tol = 1e-10 to get anything other than the second figure above, and this approach is just too slow. I hope I’ve explained my problem clearly. Thanks.