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

Post a Comment

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...
Post a Comment
Read more