# Solving Stochastic differential equation

How do I solve the following simple stochastic differential equation:

$$m x''[t] + \Gamma x'[t] + k x[t] = \sqrt[]{2 kb T/\Gamma)} \eta[t]$$

here $\eta[t]$ is brownian motion , i.e. Wiener process. The equation describes the stochastic motion of a particle in a harmonic potential.

I am new to stochastic calculus and almost know nothing. From fiddling around with Mathematica in the last two days I figured that in general the OrnsteinUhlenbeckProcess can be used to simulate $x[t]$.

The values are (converting everything in MKS) $$kb = 1.3806488 * 10 ^{-23}; (*Boltzmanns Constnat*)$$ $$T = 350.; (*K*);$$ $$r = 70.*10^{-9}; (*70 nm*);$$ $$\rho = (2./1000)*100^3; (*Density Of the Particle = 2 gm/cc*)$$ $$k = (50.*1000*2 \pi )^2 m; (* k = m \Omega^2. Resonant Frequency = 50KHz)$$ $$m = 4./3 \pi r^3 \rho; (*mass Of The Particle*)$$ $$\Gamma = 2 \pi 10^{-2};$$ these are the values I am using.

-
It would be important to know if this should be a Ito or a Stratonovich SDE. I guess it's the latter, so you probably have to use StratonovicProcess in @b.gatessucks answer below. –  sebhofer Aug 16 '13 at 9:04
However, I just realized that in your case the Ito and the Stratonovich form coincide. Still, you should be aware of the fact that there is a difference. –  sebhofer Aug 16 '13 at 10:55

I think it can be quite instructive to see how to integrate a stochastic differential equation (SDE) yourself. Of course there are different ways of doing that (a nice introduction is given in this paper). I chose the Euler-Maruyama method as it is the simplest one and is sufficient for this simple problem. Note that this assumes your SDE to be in Ito-form, which in your case coincides with the Stratonovic-form.

I write the equations of motion for the harmonic oscillator as a system of first order equations $$\dot{x}=\omega\, p,\\ \dot{p}=-\omega \,x -\gamma\, p + \xi,$$ which can easily be converted to the original equation. $\xi$ is a Wiener process which is basically just a rescaled version of $\eta$. We first sample the Wiener process from a Gaussian distribution

dt = .01; NT = 10000;
wn=Sqrt[dt] RandomVariate[NormalDistribution[0,1],NT];


and then define the update step of the Euler-Maruyama iteration

om = 1; ga = .1; n = 1;
update[x_,w_]:=(IdentityMatrix@2+{{0,om},{-om,-ga}}dt).x+Sqrt[n]{{0},{1}}w;


where n is the variance of the Wiener process. The actual integration is then just a matter of defining the initial condition and folding update over the Wiener process

x0 = {{0}, {20}};
xn = FoldList[update,x0,wn];


This yields a result similar to

ListLinePlot[{xn[[All, 1, 1]], xn[[All, 2, 1]]}, PlotRange -> All]


-
Any chance you can apply it to this ? –  b.gatessucks Aug 16 '13 at 17:29
@b.gatessucks I guess it can be, but I never worked with stochastic partial differential equations before. You probably use standard (non-stochastic) integration schemes for the spatial dimension and you use the extended (stochastic) ones for the time dimension. –  sebhofer Aug 17 '13 at 9:51

You can set it up as a system of coupled equations :

proc =
ItoProcess[{\[DifferentialD]x[t] == y[t] \[DifferentialD]t,
\[DifferentialD]y[t] == -(g/m y[t] + k/m x[t]) \[DifferentialD]t +
Sqrt[2 kb bigT/(m^2 g)] \[DifferentialD]w[t]},
x[t], {{x, y}, {x0, y0}}, {t, 0}, w \[Distributed] WienerProcess[]]


Now you can calculate different properties, like :

Mean[proc[t]]
(*
(1/(2 Sqrt[g^2 - 4 k m]))E^(-(((g + Sqrt[g^2 - 4 k m]) t)/(2 m)))
((-1 + E^((Sqrt[g^2 - 4 k m] t)/m)) g x0 + (1 +
E^((Sqrt[g^2 - 4 k m] t)/m)) Sqrt[g^2 - 4 k m] x0 + 2 (-1 + E^((Sqrt[g^2 - 4 k m] t)/m)) m y0)
*)


You can use numerical values for your parameters :

m = 1; g = 1; k = 1; kb = 1; bigT = 1;
x0 = 0; y0 = 1;
procNum =
ItoProcess[{\[DifferentialD]x[t] == y[t] \[DifferentialD]t,
\[DifferentialD]y[t] == -(g/m y[t] + k/m x[t]) \[DifferentialD]t +
Sqrt[2 kb bigT/(m^2 g)] \[DifferentialD]w[t]},
x[t], {{x, y}, {x0, y0}}, {t, 0}, w \[Distributed] WienerProcess[]]

Mean[procNum[t]]
(* (2 E^(-t/2) Sin[(Sqrt[3] t)/2])/Sqrt[3] *)

pdf[t_, x_] = Simplify[ComplexExpand[PDF[procNum[t], x], TargetFunctions -> {Re, Im}],
Assumptions -> {t > 0, x \[Element] Reals}]


Quick check :

Integrate[pdf[3, x], {x, -Infinity, Infinity}] // N
(* 1. *)

Plot[{pdf[1, x], pdf[5, x]}, {x, -5, 5}, PlotRange -> All]


With the actual values for the parameters :

g = 2 Pi 10^(-2); kb = 1.3806488 10^(-23); bigT = 350
r = 70 10^(-9); rho = (2/1000) 100^3; m = 4/3 Pi r^3 rho; k = (50 1000 2 Pi)^2 m
x0 = 0; y0 = 1;

a = Rationalize[g/m, 10^-3]
b = Rationalize[k/m, 10^-3]
c = Rationalize[Sqrt[2 kb bigT/(m^2 g)], 10^-3]

procNum2 =
ItoProcess[{\[DifferentialD]x[t] == y[t] \[DifferentialD]t,
\[DifferentialD]y[t] == -(a y[t] + b x[t]) \[DifferentialD]t + c \[DifferentialD]w[t]},
x[t], {{x, y}, {x0, y0}}, {t, 0}, w \[Distributed] WienerProcess[]]

Mean[procNum2[t]] // N
(* 4.57333*10^-17 2.71828^(-2.18659*10^16 t) (-1. 2.71828^(4.5137*10^-6 t) +
2.71828^(2.18659*10^16 t)) *)

-
Thanks. It works fine with your number. But with my numbers I get a complex infinity. Ugghh!!... Thanks very much –  user1188038 Aug 16 '13 at 16:40
Can you please edit your question and add the values for your constants ? –  b.gatessucks Aug 16 '13 at 16:52
Thanks. I just did –  user1188038 Aug 16 '13 at 17:11
Please see edit; I think it's just a numerical issue and you need to find a good compromise between accuracy and speed. –  b.gatessucks Aug 16 '13 at 17:55
I suspected the same. I didnot used this Rationalize[] function. At least now it's not throwing errors yet. Thanks. On the note: I have another question. I assume this can handle additional unitstep function or other time dependent force term as well? –  user1188038 Aug 16 '13 at 18:20