# Speeding up a differential equation system with nested modulo periodic events

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)

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

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

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


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]


• I forgot to mention, removing the restriction on the values of nmu sped the code up by 2x (from 40x faster to 75x faster). – 2012rcampion Mar 29 '15 at 4:03
• Thanks millions for this! I notice there must still be a slight mistake though, as the graph, with identical parameters, is not exactly the same - in my example I have short delay nlagdil followed by exponential increase at rate nr and 10^6 fold dilution, repeated 3 times, followed by treatment and exponential decay at rate nm and 100 fold dilution, followed by long delay nlagtreat and exponential growth at rate nlagdil, and then repeated. In your code, the mortality doesn't seem to go on for long enough and the long delay nlagtreat is missing. (Btw my actual eqn are a system of diff eqns) – Tom Wenseleers Mar 29 '15 at 9:11
• So basically I have three main types of events here - I have regrowth after dilution (here at t=0, t=1, t=2, t=5, t=6 etc), treatments (here at t=3, t=7, etc) and regrowth after treatment and dilution (here at t=4 and t=8) and the lags are different after dilution vs after treatment. Would it help by any chance to put the time values of these events in a list first, and then checking if t is a member of any of these lists to check event locations? – Tom Wenseleers Mar 29 '15 at 9:31
• The issue apparently is something with the doubleMod function - e.g. Mod[Mod[1, tg1 + ts1], tg1] == 0 gives False, whereas doubleMod[1, tg1 + ts1, tg1] gives {False,False} - I tried changing the doubleMod function to doubleMod[t_, a_, b_] /; Mod[a, b] == 0 := Mod[t, b] == 0 doubleMod[t_, a_, b_] := Max[Boole[Table[Mod[t - (i - 1) b, a] == 0, {i, Ceiling[a/b]}]]] == 1 but no joy yet... Any thoughts? – Tom Wenseleers Mar 29 '15 at 11:17
• I was thinking doubleMod[t_, a_, b_] := MemberQ[Table[Mod[t - (i - 1) b, a] == 0, {i, Ceiling[a/b]}], True] should be able to fix things, but it doesn't - any thoughts why this event then doesn't get detected, because if you plot it the value of the condition over time it shows the right values? As in ListPlot[Table[{t, Boole[doubleMod[t, tg1 + ts1, tg1]]}, {t, 0, 20, 1/10}]] which is then identical to ListPlot[Table[{t, Boole[Mod[Mod[t, tg1 + ts1], tg1] == 0]}, {t, 0, 20, 1/10}]] – Tom Wenseleers Mar 29 '15 at 11:51