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)