# Solving a differential equation with colored Gaussian noise

I'm trying to reproduce the following equation:

With $n(t)$ defined as:

Assuming $g_A=g_B=1$ I tried using NDSolve but struggle with the last term representing the colored Gaussian noise. The best thing I can come up with is using a normally distributed white noise as the last term in equation A1 above, but this seems not to be correct:

  v[t_] := RandomVariate[NormalDistribution[]]
s = NDSolve[{10 r'[t] == -4 r[t] (r[t]^2 - 1) - 2*1 (r[t] - 1) -
2*1 (r[t] + 1) + n[t],
n'[t] == -(n[t]/100) + 0.7 Sqrt[2/100]*v[t], r[0] == 0,
n[0] == 0}, {r[t], n[t]}, {t, 0, 100}];


Can someone help me out here?

• You will have to re-cast your system in terms of a random process. Have you seen this mathematica.stackexchange.com/questions/30558/… ? – Lotus Sep 26 '17 at 10:19
• Thanks! I browsed through some of the answers involving random processes, but it seems I can use the Itoprocess function only with a Wiener process but not with an Ornstein-Uhlenbeck process. I'm not familiar with random process though, so it might still work, but don't know how. – holistic Sep 26 '17 at 10:36

As @Lotus suggested, I think RandomFunction and ItoProcess are what you need.

τ = 1;
ga = gb = 1;
σ = 0.7;
τs = 100;

sol = RandomFunction[ItoProcess[{
\[DifferentialD]r[t] == \[DifferentialD]t/τ*(-4 r[t] (r[t]^2 - 1) - 2 ga (r[t] - 1) - 2 gb (r[t] + 1) + n[t]),
\[DifferentialD]n[t] == -n[t]/τs \[DifferentialD]t + σ Sqrt[2/τs] \[DifferentialD]w[t]
}, {r[t], n[t]}, {{r, n}, {0, 0}}, {t, 0}, w \[Distributed] WienerProcess[]], {0, 100, 0.01}];

ListLinePlot[sol]


Does this look right? I haven't used these functions much myself, so this is a learning experience for me too!

• Thank you, that looks similar to the results of the article. I tried it like this myself and it seems it is only possible to use a Wiener process in this function, although the authors formulate some Ornstein-Uhlenbeck process, as in equation A1. I'm gonna wait a bit longer, maybe someone can figure it out, because I have no idea ;). – holistic Oct 4 '17 at 9:22
• @holistic The OU process for n given by eqn. A1 is included in the ItoProcess. Simultaneously solving for r and n means r is driven by a OU process (the gold trajectory), unless I'm missing something. – Chris K Oct 4 '17 at 12:09
• oh I'm sorry, I must have overlooked that, you are right! Then it's all fine, awesome! – holistic Oct 5 '17 at 9:10

I have changed the definition of v[t] and got this:

Clear[t]
v[t_] := PDF[NormalDistribution[0, 0.7], t]
s = NDSolve[{10 r'[t] == -4 r[t] (r[t]^2 - 1) - 2*1 (r[t] - 1) -
2*1 (r[t] + 1) + n[t],
n'[t] == -(n[t]/100) + 0.7 Sqrt[2/100]*v[t], r[0] == 0,
n[0] == 0}, {r[t], n[t]}, {t, 0, 100}];
Plot[Evaluate[r[t] /. s], {t, 0, 100}, PlotRange -> All,
Frame -> True, FrameLabel -> {"t", "r[t]"}]
Plot[Evaluate[n[t] /. s], {t, 0, 100}, PlotRange -> All,
Frame -> True, FrameLabel -> {"t", "n[t]"}]


• Thanks for the answer. I had similar results, but this does not seem to be correct, since the results should fluctuate over time and it seems NDSolve cannot handle this random part somehow :( – holistic Oct 4 '17 at 9:23