# First order ODE with random force

I need to solve using first order Euler method some simple equation $\frac{dc}{dt} +3 c = f(t)$

where $f(t)$ is random delta-correlated force: $<f(t_1)f(t_2)>=2A \delta(t_1-t_2)$

Searching for similar problems I found that:

My code:

MyEuler[start_, end_, initialvalue_, nrOfsteps_] :=
Module[{a = start, b = end, j, m = nrOfsteps}, h = (b - a)/m;
T = Table[a + (j - 1) h, {j, 1, m + 1}];
rnd = RandomVariate[NormalDistribution[0, 1]];
Y = Table[initialvalue, {j, 1, m + 1}];
For[j = 1, j <= m, j++,
Y[[j + 1]] = Y[[j]] + h f[T[[j]], Y[[j]]]] + (2 rnd)/Sqrt[h];
Transpose@{T, Y}]
f[t_, x_] = -3 x;(*rhs of ODE without thermal noise*)
pts = MyEuler[0.0, 3.0, 10.0, 20];
ListLinePlot[pts, Mesh -> All, MeshStyle -> Red, Frame -> True]

But the results is just decaying amplitude without any random force indications in it. How should I modify this code to get correct answer? Namely $c(t)$ should relax not to zero value but to some level given by a random force amplitude I believe.

• Not sure about the meaning of $<f(t_1)f(t_2)>=2A \delta(t_1-t_2)$ but 1. + (2 rnd)/Sqrt[h] should be inside For; 2. rnd = should probably be rnd := – xzczd Dec 22 '17 at 13:35
• Thx for your suggestions. Puting random part inside of for helps, but resulting function is random, but not what I expect - just random walks around 0 value and not relaxing into some value – denkorw Dec 22 '17 at 13:51
• Have a look at ItoProcess in order to get a solution which might be a useful guide to your approach. – b.gates.you.know.what Dec 22 '17 at 14:22
• Someone can correct me if I'm wrong, but isn't this an Ornstein-Uhlenbeck Process? If so, there's a builtin OrnsteinUhlenbeckProcess function. Otherwise, as @b.gatessucks says, ItoProcess can handle more general stochastic differential equations. – Chris K Dec 22 '17 at 14:32
• If I understand the question right, you want to solve the ode c'[t]+3 c[t]==2 A DiracDelta[t-t1]? To use Euler etc. is your solution idea and not a requirement? If so, Mathematica is able to solve the problem c'[t]+3 c[t]==DiracDelta[t-t1] analytically. The solution is greensfunction which can be used to describe the general solution of your ode... – Ulrich Neumann Dec 22 '17 at 19:54

Taking into account the comments made xzczd and rewriting your code to both simplify it and make more robust, I come up with this.

myEuler[func_, start_, end_, initialvalue_, nrOfsteps_] :=
Module[{h, T, rnd, Y, j, m = nrOfsteps + 1},
rnd := RandomVariate[NormalDistribution[0, 1]];
h = (end - start)/nrOfsteps;
T = start + h (Range[m] - 1);
Y = ConstantArray[initialvalue, m];
For[j = 1, j < m, j++,
Y[[j + 1]] = Y[[j]] + h (func[T[[j]], Y[[j]]] + (2 rnd)/Sqrt[h])];
Transpose[{T, Y}]]

f[t_, x_] := -3 x + t

SeedRandom[42]
pts = myEuler[f, 0.0, 3.0, 10.0, 50];
ListLinePlot[pts,
PlotRange -> All,
Mesh -> All,
MeshStyle -> Red,
Frame -> True]

### Update

Here is a functional programming solution that greatly simplifies the function myEuler. It will produce exactly same result as the implementation using For. It may be a little harder to understand, but it is worth studying because it is very concise and much more efficient than any code building arrays and processing them with For.

myEuler[func_, start_, end_, initialvalue_, nrOfsteps_] :=
Module[{rnd, h},
rnd := RandomVariate[NormalDistribution[0, 1]];
h = (end - start)/nrOfsteps;
FoldPairList[
{{#2, #1}, #1 + h (func[#2, #1] + (2 rnd)/Sqrt[h])} &,
initialvalue,
start + h Range[0, nrOfsteps]]]
• this class of problem needs specialized methods - not my area of expertise for sure, but I do know enough that comments in the question itself are on the mark (Ito process, Ornstein-Uhlenbeck process). – Paul_A Dec 24 '17 at 2:14
• wolfram.com/broadcast/… for example (I am not affiliated with Wolfram or Pavlyk) – Paul_A Dec 24 '17 at 2:26
• @Paul_A. I am not addressing the underlying issue of whether the OP can hope to get reasonable results by employing Euler's method. I am only trying to address the programming issues his question raises. – m_goldberg Dec 24 '17 at 3:33
• I missed that, my apologies. I read many of your polished and detailed MMA programming solutions - am an admirer of yours in that regard. – Paul_A Dec 24 '17 at 4:54
• Thx for your answer, but why f[t_, x_] := -3 x + t and not just f[x]=-3x? Also from previous comments I'm still not sure how noise should be included: as h (2 rnd)/Sqrt[h] or h (2 rnd)/Sqrt[h]) – denkorw Dec 24 '17 at 12:55