6
$\begingroup$

I want to solve an ode system with pulses, for example

x'[t] == b*(1 - x[t]) - a*x[t]*y[t], 
y'[t] == a*x[t]*y[t] - d*y[t] 
in {t,0,1}

and then at the end of each period, a pulse occurs, for example, in the next period in {t,1,2},

x[0]=x[1],y[0]=y[1]*(1+x[1]),...,

and so on. I made the following code but it seems not to work.

IDESolver[{x0_, y0_}, {a_, d_, b_, T_}, 
  MaxPeriod_] := Module[{sol, resulttable},
  resulttable = {};
  Do[
    sol = 
     NDSolve[{x'[t] == b*(1 - x[t]) - a*x[t]*y[t], 
       y'[t] == a*x[t]*y[t] - d*y[t], x[0] == x0, 
       y[0] == y0}, {x, y}, {t, 0, T}];
    AppendTo[resulttable, sol];
    x[0] = Evaluate[x[T] /. sol];
    y[0] = Evaluate[(y[T]*(1+x[T])) /. sol] ;, {i, 1, MaxPeriod}
    ]
   resulttable
  ]

IDESolver[{0.2, 0.2}, {0.3, 0.1, 0.2, 1}, 10]

Can anybody help me with this issue? Thanks in advance.

$\endgroup$
  • 1
    $\begingroup$ "It seems not to work..." What about it doesn't work? Does it spit out errors? Does it spit out a wrong solution (if so, what is the solution you expect and how does it differ from the wrong solution)? Have you narrowed down where the problem is? Etc. $\endgroup$ – march Dec 18 '15 at 16:59
  • $\begingroup$ It does not work and gives errors such as "$RecursionLimit::reclim: Recursion depth of 256 exceeded. >>" $\endgroup$ – X Bruno Dec 18 '15 at 17:03
5
$\begingroup$

I made a small change and i think it is now working.

   IDESolver[{x0_, y0_}, {a_, d_, b_, T_}, MaxPeriod_] := 
      Block[{X0 = x0, Y0 = y0, sol, resulttable}, resulttable = {};
       Do[
        sol = 
         NDSolve[{x'[t] == b*(1 - x[t]) - a*x[t]*y[t], 
           y'[t] == a*x[t]*y[t] - d*y[t], x[0] == X0, y[0] == Y0}, {x, 
           y}, {t, 0, T}];
        AppendTo[resulttable, sol];

        X0 = Evaluate[x[T] /. sol];
        Y0 = Evaluate[(y[T]*(1 + x[T])) /. sol];, {i, 1, MaxPeriod}] ;
       Flatten[resulttable]
       ];
    listsols = IDESolver[{0.2, 0.2}, {0.3, 0.1, 0.2, 1}, 10];
    Table[Plot[listsols[[i]][[2]][x], {x, 0, 10}], {i, 1, 
      Length[listsols]}]
$\endgroup$
  • $\begingroup$ Thank you very much! It still gives errors such as "$RecursionLimit::reclim: ""Recursion depth of 256 exceeded." PS. I am using Mathematica 7 on Mac. $\endgroup$ – X Bruno Dec 18 '15 at 18:05
  • $\begingroup$ I am using Mathematica 10 on Windows 8 and its is working fine... $\endgroup$ – Diogo Dec 18 '15 at 18:08
  • $\begingroup$ I restart the Kernel and it is working! Thanks a lot! $\endgroup$ – X Bruno Dec 18 '15 at 18:16
  • $\begingroup$ You are welcome! $\endgroup$ – Diogo Dec 18 '15 at 18:18
7
$\begingroup$

In the case where you actually want a continuous time variable, the function WhenEvent will be extremely useful. I've modified your code to:

Clear[IDESolver, x, y, sol]
IDESolver[{x0_, y0_}, {a_, d_, b_, period_}, tMax_] :=
  NDSolve[{
    x'[t] == b*(1 - x[t]) - a*x[t]*y[t]
    , y'[t] == a*x[t]*y[t] - d*y[t]
    , x[0] == x0, y[0] == y0
    , WhenEvent[Mod[t, period] == 0, y[t] -> y[t] (1 + x[t])]
   }
   , {x, y}
   , {t, 0, tMax}]

Then,

sol = First@IDESolver[{0.2, 0.2}, {0.3, 0.1, 0.2}, 1, 10];
Plot[Evaluate[# /. sol], {t, 0, 10}] & /@ {x[t], y[t]} // GraphicsRow

enter image description here

$\endgroup$
  • $\begingroup$ @SteveX. You're welcome! Don't forget to accept the answer that you find most helpful by clicking the grey checkmark next to that answer! $\endgroup$ – march Dec 21 '15 at 19:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.