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

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

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more