5
$\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$
7
  • $\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
    Commented Dec 22, 2017 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
    Commented Dec 22, 2017 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$ Commented Dec 22, 2017 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
    Commented Dec 22, 2017 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$ Commented Dec 22, 2017 at 19:54

1 Answer 1

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$
6
  • $\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
    Commented Dec 24, 2017 at 2:14
  • $\begingroup$ wolfram.com/broadcast/… for example (I am not affiliated with Wolfram or Pavlyk) $\endgroup$
    – Paul_A
    Commented Dec 24, 2017 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
    Commented Dec 24, 2017 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
    Commented Dec 24, 2017 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
    Commented Dec 24, 2017 at 12:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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