# Defining stochastic differential equation & simulating a system of three SDEs

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?

• 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! Nov 5 '15 at 16:22
• Also, what's wrong with directly substituting $(dX)_g$ into the other two equations? Am I missing something? 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[%]

• Good you figured it out, but you still haven't posted a sensible set of values, neither in the question nor in your answer! Nov 5 '15 at 20:51
• @ 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[%] ? Nov 5 '15 at 21:31
• 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]"}] Nov 5 '15 at 22:14
• 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! Nov 6 '15 at 10:28
• 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. Nov 6 '15 at 10:36