I am no expert on SDE but I've been messing around with MMA's built in functions and it makes it quite easy to do some simple simulations. I bumped into this system of equations (below) in a paper and was trying to tweak the documentation center equations to simulate it.

It describes biological growth of biomass $X$, in a reactor that is fed substrate $S$, and the growth of bacteria is defined by $(dX)_g$

$$(dX)_g= \mu S/(K_j+S) *Xdt+\sigma dW,$$

$$dX=(dX)_g-Q_eXdt+(Q_rX_r-Q_oX)dt+\sigma_i dW$$

$$dS=-1/Y (dX)_g+Q_f(S_f-S)dt+\sigma_j dW $$

And $W, W_i, W_j$ are Wiener processes, $\sigma$s are volatility and $\mu$ is defined as the specific growth rate but represents drift.

Where the coefficients: $S_f, X_r, Y, Q_o,Q_f,Q_r,Q_e$ are scalars and constants (I can provide what they stand for is requested).

My approach was to define $(dX)_g$ first, simulate it and call it in the equations for $dX$ and $dS$: Qe = 10; Qr = 10; Xr = 30; Qo = 10; [Sigma] = 10; [Sigma]2 = 20; [Mu] = 0.25; Y = 0.25; Qf = 10; Sf = 300; Ks = 200;

a = ItoProcess[\[DifferentialD]xg[t] == \[Mu]*s[t]/(Ks + s[t])*x[t] \[DifferentialD]t + \[Sigma] \[DifferentialD]w[t], xg[t], {xg, 25}, t, w \[Distributed] WienerProcess[]];
aproc = RandomFunction[%, {0., 5., 0.01}]

And then call the function aproc I defined above in the two equations for $dX$ and $dS$:

b = ItoProcess[\[DifferentialD]s[t] == -1/Y (\[Mu]*s[t]/(Ks + s[t])*x[t] \[DifferentialD]t + \[Sigma] \[DifferentialD]w[t]) +Qf (Sf - s[t]) \[DifferentialD]t + \[Sigma] \[DifferentialD]w[t], {s[t], x[t]}, {s, 200}, t, w \[Distributed] WienerProcess[]];
bproc = RandomFunction[%, {0., 5., 0.01}];
c = ItoProcess[\[DifferentialD]x[t] == (\[Mu]*s[t]/(Ks + s[t])*x[t] \[DifferentialD]t + \[Sigma] \[DifferentialD]w[t]) - Qe*x[t] \[DifferentialD]t +(Qr*Xr - Qo*x[t]) \[DifferentialD]t + \[Sigma] \[DifferentialD]w[t], {x[t],s[t]}, {x, 100}, t, w \[Distributed] WienerProcess[]];
cproc = RandomFunction[%, {0., 5., 0.01}];

but I need to add the $(dX)_g$ into the equation before Qe*x[t] \[DifferentialD]t I tried defining its SDE separately and simulating it with RandomFunction then inserting it into the $dX$ equation above and clearly saw that's not the way to do this...

So my question is: How do I write the system of SDE above? And perhaps shed some insight as to whether the ItoProcess is what I want to use?

  • 2
    $\begingroup$ Please ALWAYS provide a COMPLETE code sample! Don't expect others to make up their own set of parameters if they want to help with your question! $\endgroup$ – sebhofer Nov 5 '15 at 16:22
  • 2
    $\begingroup$ Also, what's wrong with directly substituting $(dX)_g$ into the other two equations? Am I missing something? $\endgroup$ – sebhofer Nov 5 '15 at 16:24

Got it! Thanks @sebhofer.

My mistake was trying to simulate ItoProcess with RandomFunction and then calling it into another ItoProcess simulation... But the SDE simulations don't work that way. Instead I used the ItoProcess function for a vector process $x[t], s[t]$; inserted $(dW)_g$ directly into the equations and ended up with a 2-equation system of SDEs:

Qe = 10; Qr = 10; Xr = 30; Qo = 10; σ = 10; σ2 = 20; μ = 0.25; Y = 0.25; Qf = 10; Sf = 300; Ks = 200;

proc = ItoProcess[{\[DifferentialD]x[t] == (μ s[t]/(Ks + s[t]) x[t] \[DifferentialD]t + σ \[DifferentialD]w1[t]) - Qe x[t] \[DifferentialD]t + (Qr*Xr -Qo*x[t]) \[DifferentialD]t + σ2 \[DifferentialD]w1[t], \[DifferentialD]s[t] == -1/Y (μ s[t]/(Ks + s[t]) x[t] \[DifferentialD]t + σ \[DifferentialD]w2[t]) + Qf (Sf - s[t]) \[DifferentialD]t + σ2 \[DifferentialD]w2[t]},
{x[t], s[t]}, {{x, s}, {25, 800}}, t, {w1 \[Distributed] WienerProcess[], 
w2 \[Distributed] WienerProcess[]}];
RandomFunction[%, {0., 5., 0.001}];
ListLinePlot[%, PlotLegends -> {"x[t]", "s[t]"}]

This simulates both equations and then I can plot using ListLinePlot[%]

  • $\begingroup$ Good you figured it out, but you still haven't posted a sensible set of values, neither in the question nor in your answer! $\endgroup$ – sebhofer Nov 5 '15 at 20:51
  • $\begingroup$ @ sebhofer How come when I try to simulate 20 paths with RandomFunction[%,{0,5,0.01},20], I see TemporalData be generated by nothing appears when I try to plot it with ListLinePlot[%] ? $\endgroup$ – E3labs Nov 5 '15 at 21:31
  • $\begingroup$ If we set sample = RandomFunction[%, {0., 5., 0.001}]; then sample["ValueDimensions"] will return 2, hence ListLinePlot[%] above returns nothing. Try ListLinePlot[{sample["PathComponent", 1], sample["PathComponent", 2]}, PlotLegends -> {"x[t]", "s[t]"}] $\endgroup$ – Gosia Nov 5 '15 at 22:14
  • $\begingroup$ Yeah, this is not quite obvious. You should look at the TemporalData documentation for all the options. You can select one of the sampled trajectories by doing something like rf=RandomFunction[proc, {0., 5., 0.001}, 20]; rf["Path", 1], which will give you a List of raw data. Note, however, that this is not in the correct format to use with ListPlot directly. So this solution is rather clunky and @Gosia's solution is probably better! $\endgroup$ – sebhofer Nov 6 '15 at 10:28
  • 1
    $\begingroup$ I just realized that @Gosia's solution does not work for me. You can do rf=RandomFunction[proc, {0., 5., 0.001},20]; Show[ListPlot[rf["PathComponent", 1]], ListPlot[rf["PathComponent", 2]], PlotRange -> All] though. $\endgroup$ – sebhofer Nov 6 '15 at 10:36

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.