1
$\begingroup$

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?

$\endgroup$
  • 2
    $\begingroup$ 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? $\endgroup$ – march Dec 14 '15 at 23:02
  • 4
    $\begingroup$ "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 $\endgroup$ – Dr. belisarius Dec 15 '15 at 0:43
  • 1
    $\begingroup$ @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. $\endgroup$ – march Dec 15 '15 at 2:57
  • 1
    $\begingroup$ 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}] $\endgroup$ – Dr. belisarius Dec 15 '15 at 14:25
  • 2
    $\begingroup$ 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] $\endgroup$ – Dr. belisarius Dec 15 '15 at 14:41

Your Answer

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

Browse other questions tagged or ask your own question.