4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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 := $\endgroup$ – xzczd Dec 22 '17 at 13:35
  • $\begingroup$ 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 $\endgroup$ – denkorw Dec 22 '17 at 13:51
  • 3
    $\begingroup$ Have a look at ItoProcess in order to get a solution which might be a useful guide to your approach. $\endgroup$ – b.gates.you.know.what Dec 22 '17 at 14:22
  • 1
    $\begingroup$ 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. $\endgroup$ – Chris K Dec 22 '17 at 14:32
  • 3
    $\begingroup$ 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... $\endgroup$ – Ulrich Neumann Dec 22 '17 at 19:54
5
$\begingroup$

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]

plot

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]]]
$\endgroup$
  • $\begingroup$ 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). $\endgroup$ – Paul_A Dec 24 '17 at 2:14
  • $\begingroup$ wolfram.com/broadcast/… for example (I am not affiliated with Wolfram or Pavlyk) $\endgroup$ – Paul_A Dec 24 '17 at 2:26
  • $\begingroup$ @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. $\endgroup$ – m_goldberg Dec 24 '17 at 3:33
  • $\begingroup$ I missed that, my apologies. I read many of your polished and detailed MMA programming solutions - am an admirer of yours in that regard. $\endgroup$ – Paul_A Dec 24 '17 at 4:54
  • $\begingroup$ 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]) $\endgroup$ – denkorw Dec 24 '17 at 12:55

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.