Speeding up a differential equation system with nested modulo periodic events

Question

I have the following differential equation system with periodic Events :

tg1 = 3; ts1 = 1; td1 = 1; nlagdil = 1/5; nlagtreat = 1/2;
s[nr_, nm_] := 
  First[ Evaluate[
    NDSolve[{n'[t] == nmu[t]*n[t], n[0] == 1, nmu[0] == 0, 
      timedil[0] == 0, timetreat[0] == tg1,
      WhenEvent[Mod[Mod[t, tg1 + ts1], td1] == 0, {n[t] -> n[t]/10^6, nmu[t] -> 0, timedil[t] -> t}], 
      WhenEvent[t - timedil[t] == nlagdil, {nmu[t] -> nr}], 
      WhenEvent[Mod[Mod[t, tg1 + ts1], tg1] == 0, {n[t] -> n[t], nmu[t] -> nm, timetreat[t] -> t}], 
      WhenEvent[t - timetreat[t] == ts1, {n[t] -> n[t]/100, nmu[t] -> 0}],
      WhenEvent[t - timetreat[t] == (ts1 + nlagtreat), {nmu[t] -> nr}] 
      },
     {n, nmu, timedil, timetreat}, {t, 0, 20}, 
     DiscreteVariables -> {nmu[t] \[Element] {0, nr, nm}, timedil[t], 
       timetreat[t]}, StartingStepSize -> 1/100, 
     Method -> "StiffnessSwitching", AccuracyGoal -> 100, 
     WorkingPrecision -> 100, MaxSteps -> Infinity, 
     MaxStepFraction -> Infinity]]];
soln = s[10, -20];
LogPlot[Evaluate[{n[t]} /. soln], {t, 0, 20}, PlotRange -> All, 
 Frame -> True, FrameLabel -> {"Time", "n"}, 
 PlotTheme -> {"Classic", "LargeLabels"}] 

(Well my actual system is larger, but this can serve as an example)

Problem is that I can only get Mathematica to detect the events with AccuracyGoal -> 100 and WorkingPrecision -> 100 which makes it really slow to solve. Does anybody have any thoughts how I could speed this up? Any other options or settings that would work for this?

EDIT: the answer below points out that the problem has to do with the nested Mod[] conditions. The solution below, unfortunately doesn't work since the proposed real valued function doubleMod[t_,a_,b_] is not fully equivalent to Mod[Mod[t,a],b] e.g. Mod[Mod[1, tg1 + ts1], tg1] == 0 gives False, whereas doubleMod[1, tg1 + ts1, tg1] gives {False,False}. As an alternative way to get rid of the nested Mod[] I tried

s[nr_, nm_] := 
  First@NDSolve[{n'[t] == nmu[t]*n[t], n[0] == 1, nmu[0] == 0, 
     timedil[0] == 0, timetreat[0] == tg1, timestartround[0] == 0,
     WhenEvent[
      Mod[t, tg1 + ts1] == 
       0, {timestartround[t] -> Rationalize[t, 0.01]}],
     WhenEvent[
      Mod[t - timestartround[t], td1] == 0, {n[t] -> n[t]/10^6, 
       nmu[t] -> 0, timedil[t] -> t}],
     WhenEvent[t - timedil[t] == nlagdil, {nmu[t] -> nr}, 
      "Priority" -> 3],
     WhenEvent[
      Mod[t - timestartround[t], tg1] == 0, {n[t] -> n[t], 
       nmu[t] -> nm, timetreat[t] -> t}],
     WhenEvent[
      t - timetreat[t] == ts1, {n[t] -> n[t]/100, nmu[t] -> 0}],
     WhenEvent[t - timetreat[t] == (ts1 + nlagtreat), {nmu[t] -> nr}]},
    {n, nmu, timestartround, timedil, timetreat}, {t, 0, 20},
    DiscreteVariables -> {nmu[t], timestartround[t], timedil[t], 
      timetreat[t]}, Method -> "StiffnessSwitching", 
    StartingStepSize -> 1/100, AccuracyGoal -> 70, 
    WorkingPrecision -> 50, MaxSteps -> Infinity, 
    MaxStepFraction -> Infinity];
AbsoluteTiming[soln = s[10, -20];]
LogPlot[Evaluate[{n[t]} /. soln], {t, 0, 20}, PlotRange -> All, 
 Frame -> True, FrameLabel -> {"Time", "n"}, 
 PlotTheme -> {"Classic", "LargeLabels"}] 

But this still requires the use of method StiffnessSwitching, and requires high accuracy and working precision, so is also really slow. Does anybody have any other thoughts about how to fix this? Basically I have several growth cycles (nr=tg1/td1) followed by an antibiotic treatment cycle of length ts1, which are then further iterated. Other possibility perhaps could be to recursively iterate this system, but I was wondering how to do this (ie how to patch together the different sections of the differential equation) and still obtain the same plot as shown above? (Obviously my real system is much more complicated, but I am mainly testing the events here)

Answer 1

The functions that are giving you trouble are the nested Mods, like this one:

Mod[Mod[t, tg1 + ts1], td1] == 0

The documentation for WhenEvent says:

• WhenEvent typically finds a function f of the solution variables such that the roots of f occur where the event becomes True...
• In the case that a root function cannot be determined, the bracketing interval is found using bisection.

However, the Mod functions never actually pass through zero, so the numerical algorithm has trouble finding their locations: since it doesn't see the sign of the function change, it thinks that it couldn't have passed through zero and never attempts bisection. In fact, the first WhenEvent should be occurring every 1 time unit (its condition is Mod[Mod[t,4],1], equivalent to Mod[t,1]) but clearly is being skipped in your plot.

Fortunately, there is a solution. The documentation list a special type of event condition:

Mod[t, Δt] == 0: sample at regular intervals Δt in the time variable t.

In this case, we can split the nested Mods apart using a helper function:

doubleMod[t_, a_, b_] /; Mod[a, b] == 0 := Mod[t, b] == 0
doubleMod[t_, a_, b_] := 
 Table[Mod[t - (i - 1) b, a] == 0, {i, Ceiling[a/b]}]

Example usage:

Mod[Mod[t, tg1 + ts1], td1] == 0
doubleMod[t, tg1 + ts1, td1]
(* Mod[Mod[t, 4], 1] == 0 *)
(* Mod[t, 1] == 0 *)

Mod[Mod[t, tg1 + ts1], tg1] == 0
doubleMod[t, tg1 + ts1, tg1]
(* Mod[Mod[t, 4], 3] == 0 *)
(* {Mod[t, 4] == 0, Mod[-3 + t, 4] == 0} *)

Your original code takes around 1.2 seconds on my computer, but when revised using the doubleMod helper:

s[nr_, nm_] := 
  First@NDSolve[{n'[t] == nmu[t]*n[t], n[0] == 1, nmu[0] == 0, 
     timedil[0] == 0, timetreat[0] == tg1, 
     WhenEvent[
      Evaluate@doubleMod[t, tg1 + ts1, td1], {n[t] -> n[t]/10^6, 
       nmu[t] -> 0, timedil[t] -> t}], 
     WhenEvent[t - timedil[t] == nlagdil, {nmu[t] -> nr}], 
     WhenEvent[
      Evaluate@doubleMod[t, tg1 + ts1, tg1], {n[t] -> n[t], 
       nmu[t] -> nm, timetreat[t] -> t}], 
     WhenEvent[
      t - timetreat[t] == ts1, {n[t] -> n[t]/100, nmu[t] -> 0}], 
     WhenEvent[
      t - timetreat[t] == (ts1 + nlagtreat), {nmu[t] -> nr}]}, {n, 
     nmu, timedil, timetreat}, {t, 0, 20}, 
    DiscreteVariables -> {nmu[t], timedil[t], timetreat[t]}];
AbsoluteTiming[soln = s[10, -20];]

It takes only 0.016 seconds (and hits all the events). No extra options are needed. (I usually find that if I need to pass lots of extra options to NDSolve then I'm doing something wrong.)

Note that in this case, the events do not depend on n[t] and n obeys a trivial exponential model. It may be even faster to remove n from NDSolve, compile a list of events that set n and nmu, and compute n from the breakpoints afterward. Removing n, however, leaves only DiscreteVariables remaining so you don't even really need NDSolve for this problem. If you need speed, it will be much faster to roll your own algorithm: at each step, compute when each event will occur next, step directly to the nearest event, increment your variables appropriately, then repeat.

Update

Based on our discussion, I've come up with a method of programmatically generating the WhenEvents to avoid the need for the nested Mods in the first place (which weren't behaving like you expected anyway). First we set up the definitions of your tests:

The tests repeat in the pattern defined by testSchedule, which is a list of indexes in the following arrays.

testSchedule = {1, 1, 1, 2, 3};

Then we define the tests; for each test, dilution by dilutionFactor occurs immediately, then after growthDelay time units, the growth factor is set to growthFactor for growthTime time units before beginning the next test.

dilutionFactor = {1/10^6, 1, 1/100};
growthDelay = {0.2, 0.5, 0.5};
growthFactor = {10, -20, -20};
growthTime = {1, 1, 1};
(*sanity checks*)
Length[dilutionFactor] == Length[growthDelay] ==
  Length[growthFactor] == Length[growthTime]
And @@ Map[IntegerQ[#] && 1 <= # <= Length[dilutionFactor] &, testSchedule]
And @@ Thread[growthDelay > 0]
And @@ Thread[growthTime > 0]

Now we can compute the event times:

With[{df = dilutionFactor[[testSchedule]],
  gd = growthDelay[[testSchedule]],
  gf = growthFactor[[testSchedule]],
  gt = growthTime[[testSchedule]]},
 tDilution = Most[Accumulate[gd] + Accumulate[gt]]~Prepend~0
    tGrowth = Accumulate[gd] + Most[Accumulate[gt]]~Prepend~0
      tRepeat = Total[gd] + Total[gt]
       events = Join @@ Table[{
        Evaluate /@ 
         WhenEvent[
          Mod[t, tRepeat] == tDilution[[i]], {n[t] -> n[t]*df[[i]], 
           nmu[t] -> 0}],
        Evaluate /@ 
         WhenEvent[
          Mod[t, tRepeat] == tGrowth[[i]], {nmu[t] -> gf[[i]]}]
        }, {i, Length[testSchedule]}]
 ]

Finally we solve the differential equation.

soln = NDSolve[{n'[t] == n[t] nmu[t], n[0] == 1, nmu[0] == 0}~Join~
   events, {n, nmu}, {t, 0, 20}, DiscreteVariables -> {nmu[t]}, 
  MaxStepSize -> (Min[growthDelay, growthTime]/10)]
LogPlot[Evaluate[{n[t]} /. soln], {t, 0, 2 tRepeat}, 
 PlotPoints -> 1000]
Plot[Evaluate[{nmu[t]} /. soln], {t, 0, 2 tRepeat}, 
 PlotPoints -> 1000]