# How to make a parameter stochastic in a differential equation system with NDSolve?

I constructed the differential equation system below, which I solved using NDSolve. Now I need one parameter ($mu$) to be stochastic, e.g. Poisson distributed around a mean and changing slightly at every time step. This is my current model, which is now fully deterministic:

sol[fmaxnW_?NumericQ, fminnW_?NumericQ, fmaxpW_?NumericQ,fmaxnR_?
NumericQ,fminnR_?NumericQ, micW_?NumericQ, micR_?NumericQ,kW_?NumericQ, kR_?
NumericQ , K_?NumericQ, a_?NumericQ, b_?NumericQ,mu_?NumericQ, Amax_?
NumericQ,w_?NumericQ, pinit_?NumericQ, tinit_?NumericQ, propR_?NumericQ] :=
Module[{},
s = NDSolve[
{A'[t] == -w*A[t],
nW'[t] == (fmaxnW - ((fmaxnW - fminnW)*(A[t]/micW)^kW)/((A[t]/micW)^kW - fminnW/fmaxnW))*
Which[A[t] < micW, (1 - 1/K (nW[t] + pW[t] + nR[t])), True, 1]*nW[t] -
a*1/K (nW[t] + pW[t] + nR[t])*nW[t]  +  b*pW[t] -
w*nW[t] - mu*(fmaxnW - ((fmaxnW - fminnW)*(A[t]/micW)^kW)/((A[t]/micW)^
kW - fminnW/fmaxnW))*nW[t],

pW'[t] == fmaxpW*pW[t]  +  a*1/K (nW[t] + pW[t] + nR[t])*nW[t]  - b*pW[t] -
w*pW[t] - mu*fmaxpW*pW[t],

nR'[t] == (fmaxnR - ((fmaxnR - fminnR)*(A[t]/micR)^kR)/((A[t]/micR)^kR - fminnR/fmaxnR))*
Which[A[t] < micR, (1 - 1/K (nW[t] + pW[t] + nR[t])), True, 1]*nR[t] -
w*nR[t] + mu*(fmaxnW - ((fmaxnW - fminnW)*(A[t]/micW)^kW)/((A[t]/micW)^kW -
fminnW/fmaxnW))*nW[t] + mu*fmaxpW*pW[t],

A[0] == 0, WhenEvent[Mod[t, 24] == 0, A[t] -> Amax],
pW[0] == pinit*(1 - propR),
nW[0] == (tinit - pinit)*(1 - propR),  nR[0] == tinit*propR
},
{A, nW, pW, nR}, {t, 0, 240}];
Return[s]]

sol[1, -6, 0.00001, 0.9, -6, 0.01, 0.05, 1, 1,
5*10^9, 0.0001, 0.14, 8*10^(-8), 0.02, 0.231, 1, 2*10^5, 0]


How can I implement stochastic variation in the model parameter $mu$ in a simple way?

• Can you be more specific? Do you want mu to be a different (random) value in each call to sol? Or do you want mu to change at every time step in NDSolve? Or change at random times in NDSolve? – march Dec 14 '15 at 23:02
• "and changing slightly at every time step" isn't well defined because:1) Time steps can be adaptative and 2) "changing slightly" doesn't have a precise meaning – Dr. belisarius Dec 15 '15 at 0:43
• @belisariushassettled. When I phrased it in that way, I had in mind numerical solutions to stochastic differential equations where the time step is fixed (non-adaptive), and the mu might some be something like a Wiener increment. I suspect that's not actually what the OP wants, but, you know, I figured I'd throw it out there. – march Dec 15 '15 at 2:57
• Ok, voting to reopen. Here's a way: sol = NDSolve[{x'[t] == Sin[w[t] t], x[0] == 0, w[0] == 0, WhenEvent[Mod[t, .1] == 0, w[t] -> RandomVariate[NormalDistribution[]]]}, x, {t, 10}, DiscreteVariables -> {w}]; Plot[x[t] /. sol, {t, 0, 10}] – Dr. belisarius Dec 15 '15 at 14:25
• Or better: sol := NDSolveValue[{x'[t] == Sin[w[t] t], x[0] == 0, w[0] == 0, WhenEvent[Mod[t, .1] == 0, w[t] -> RandomVariate[NormalDistribution[]]]}, x, {t, 10}, DiscreteVariables -> {w}]; Plot[Through[Flatten@Table[sol, {6}][t]], {t, 0, 10}, Evaluated -> True] – Dr. belisarius Dec 15 '15 at 14:41