Skip to main content

differential equations - Having trouble using WhenEvent with NDSolve


First of all, this is my first time on this forum. Please be patient with me because if I make some mistakes in this post. I'm not used to doing this!


I'm trying to modify a parameter during NDSolve iterations. I tried to follow advice given in the following exchanges:




But none of it seems to fit my problem. Perhaps I'm doing something wrong.


Here is my problem:


I'm currently trying to create a model with the following differential equations:


dEc = p1*Ec[t]*(1 - (Ec[t]/Emax)) - d1*Ec[t] - i1*Ic[t]*Ec[t];
dIc = i1*Ic[t]*Ec[t] - u1*Ic[t];
dVp = ist*b1*u1*Ic[t] - c1*Vp[t];

The parameters are defined as following:


parameters = 
{Emax -> 10000, p1 -> 0.6, d1 -> 0.003, u1 -> 0.33, c1 -> 10,

b1 -> 6000, i1 -> 0.0000002, ist -> 0.0001};

Then, I'm doing the following NDSolve:


dynamicsmodel = 
NDSolve[
Evaluate[
{Ec'[t] == dEc, Ic'[t] == dIc, Vp'[t] == dVp,
Ec[0] == 10000, Ic[0] == 0, Vp[0] == 10}] /. parameters,
{Ec[t], Ic[t], Vp[t]},
{t, 0, 300}]


It seems to work perfectly. But I would like to add an event listener that would change the ist parameter when t >= 200. To do so, I tried with the following code:


dynamicsmodel = 
NDSolve[
Evaluate[
{Ec'[t] == dEc, Ic'[t] == dIc, Vp'[t] == dVp,
Ec[0] == 10000, Ic[0] == 0, Vp[0] == 10}] /. parameters,
WhenEvent[t >= 200, parameters[[8]] = ist -> 0.01],
{Ec[t], Ic[t], Vp[t]},
{t, 0, 300}]


But it doesn't seem to work! I get the following error:



To avoid possible ambiguity, the arguments of the dependent variable in WhenEvent[t >= 200, parameters[[8]] = ist -> 0.01] should literally match the independent variables."



I may have misplaced the WhenEvent. I tried to put it inside of the Evaluate expression, but that produces errors.


Please note that I simplified the code to make it more readable. I hope I didn't introduce errors.



Answer



I believe you need to take ist out of parameters and make it a DiscreteVariable in the NDSolve. This seems to work:


dEc := p1*Ec[t]*(1 - (Ec[t]/Emax)) - d1*Ec[t] - i1*Ic[t]*Ec[t];

dIc := i1*Ic[t]*Ec[t] - u1*Ic[t];
dVp := ist[t]*b1*u1*Ic[t] - c1*Vp[t];

parameters = {Emax -> 10000, p1 -> 0.6, d1 -> 0.003, u1 -> 0.33, c1 -> 10, b1 -> 6000, i1 -> 0.0000002};

dynamicsmodel = NDSolve[{
Ec'[t] == dEc, Ic'[t] == dIc, Vp'[t] == dVp,
Ec[0] == 10000, Ic[0] == 0, Vp[0] == 10, ist[0] == 0.0001,
WhenEvent[t == 200, ist[t] -> 0.01]} /. parameters,
{Ec, Ic, Vp, ist}, {t, 0, 300}, DiscreteVariables -> {ist}][[1]];


GraphicsGrid[{{
Plot[Evaluate[Ec[t] /. dynamicsmodel], {t, 0, 300}, PlotRange -> All, PlotLabel -> Ec],
Plot[Evaluate[Ic[t] /. dynamicsmodel], {t, 0, 300}, PlotRange -> All, PlotLabel -> Ic]},
{Plot[Evaluate[Vp[t] /. dynamicsmodel], {t, 0, 300}, PlotRange -> All, PlotLabel -> Vp],
Plot[Evaluate[ist[t] /. dynamicsmodel], {t, 0, 300}, PlotRange -> All, PlotLabel -> ist]
}}, ImageSize -> 600]

Mathematica graphics


Kind of underwhelming though -- maybe different parameters would be more interesting!



As a side note, I like to omit the [t] in the list of dependent variables.


Comments