2
$\begingroup$

I am struggling with writing a stochastic version of Lotka-Volterra predator-prey model. This is as far as I have gotten:

a = 0.1;
b = 0.2;
c = 0.3;
d = 0.4;
proc = 
  RandomFunction[
    ItoProcess[
      {\[DifferentialD]X[t] == (a \[DifferentialD]t X[t] + \[DifferentialD]W[t]) - 
         b \[DifferentialD]t X[t] Y[t],
       \[DifferentialD]Y[t] == -c \[DifferentialD]t Y[t] + 
       d \[DifferentialD]t Y[t] X[t]}, 
      {X[t], Y[t]}, {{X, Y}, {0.2, 0.2}}, t, 
      W \[Distributed] WienerProcess[0, 0.1]], {0, 100, 0.01}, 
      Method -> "StochasticRungeKuttaScalarNoise"];
ListLinePlot[proc, PlotRange -> All]

For some reason, only one variable depends on the stochastic noise and it does not bear influence on the other part of the equation, I have tried to correct it, but no matter what I do, no luck.

First thing would be to get it working correctly, then I would worry about excluding "extinction Values > 0", since we are not interested in those.

Also, in regards to the graph, is there other way I would create a phase plot (X against Y)? My current way:

ListLinePlot[proc["Values"], PlotRange -> All]

plots X and Y, against time. How would I plot only X in time?

$\endgroup$
  • $\begingroup$ As to the 2nd question, you might try ListLinePlot[Transpose@proc["Values"], PlotRange -> All] $\endgroup$ – m_goldberg Jan 20 '18 at 7:12
  • $\begingroup$ Oh right, that plots amplitude of both functions (x and y) in time, on the same graph. Would there be a way to separate those functions, and create two graphs, first with only function X values in time, and second with only function Y values in time? @m_goldberg $\endgroup$ – Peter Jan 20 '18 at 7:42
2
$\begingroup$

There are a few questions wrapped up here. Let me try to take them one-by-one.

only one variable depends on the stochastic noise and it does not bear influence on the other part of the equation

I don't think that's actually correct. After running your code,

ListLinePlot[proc, PlotRange -> All]

Mathematica graphics

The predator (gold) seems to smooth over the variation in the prey (blue) but take a look around t=85: there is a random bump in prey which results in a little shoulder in the decline of the predators. So I think the simulation results are accurate for this model and parameter set.

then I would worry about excluding "extinction Values > 0", since we are not interested in those.

That prey can become negative is a problem with the model formulation. I don't think RandomFunction[ItoProcess has the equivalent of WhenEvent in NDSolve.

Two ideas: 1) replace your additive noise \[DifferentialD]W[t] with a multiplicative noise X[t] \[DifferentialD]W[t], which seems like it should prevent the prey population from becoming negative. 2) Instead of using a stochastic differential equation as your model, you could use a more basic continuous time stochastic process and simulate it with @IstvánZachar's GillespieSSA function from this answer (specifically, see his example 2).

is there other way I would create a phase plot (X against Y)?

Funny, when I run your code, it works as desired (v11.2):

ListLinePlot[proc["Values"], PlotRange -> All]

Mathematica graphics

How would I plot only X in time?

This recent answer by @kglr addresses this point, using an undocumented property "PathComponent":

ListLinePlot[proc["PathComponent", 1]]

Mathematica graphics

Anyone want to write a wrapper (call it SNDSolve) for RandomFunction[ItoProcess that mimics NDSolve in input and output? Seems like it'd be appreciated by many!

$\endgroup$
  • $\begingroup$ Or maybe NSDSolve makes more sense ;) $\endgroup$ – Chris K Jan 21 '18 at 16:40
  • $\begingroup$ Thank you very much for the explanation. In regards to the phase graph, I was just wondering whether mine was a correct way of doing it, since I was not able to find "Values" or "PathComponent" anywhere online... I will try to tweak it for the extinction event case. Once again, thank you! $\endgroup$ – Peter Jan 22 '18 at 1:48
  • $\begingroup$ If this answers your question, you can accept it by clicking on the check mark beside the answer to toggle it from greyed out to filled in. $\endgroup$ – Chris K Jan 22 '18 at 14:23
1
$\begingroup$

This solution will first define the ItoProcess and then simulate it; care is taken to initialize all the available parameters at their appropriate values (using With):

With[{a = 0.1, b = 0.2, c = 0.3, d = 0.4, Xo = 0.2, Yo = 0.2, wm = 0.,ws = 0.1},
 proc = ItoProcess[{
    \[DifferentialD]X[t] == (a - b Y[t]) X[t] \[DifferentialD]t + \[DifferentialD]W[t],
    \[DifferentialD]Y[t] == (-c + d X[t]) Y[t] \[DifferentialD]t
   }, {X[t], Y[t]}, {{X, Y}, {Xo, Yo}},
   t, W \[Distributed] WienerProcess[wm, ws]
  ]
 ];

Next, we simulate 3 paths for the process (proc):

td = RandomFunction[proc, {0, 100, 0.01}, 3, Method -> "StochasticRungeKuttaScalarNoise"]

Finally we plot the derived paths:

Grid[
 Partition[
  MapIndexed[

   With[{val = Part[#1, All, -1]},
     ListLinePlot[
      Transpose[val],
      PlotRange -> All,
      PlotLegends -> {X[t], Y[t]},
      PlotLabel -> Row[{Path, Null, #2[[-1]]}]
    ]] &, td["Paths"]], 2, 2, {1, 1}, Null]]

enter image description here

$\endgroup$

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.