Solving a stochastic differential equation

How do I solve the following simple stochastic differential equation:

$$m x''[t] + \Gamma x'[t] + k x[t] = \sqrt[]{(2 k_{b} 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 know almost 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 I'm using are (converting everything in MKS):

$$k_{b} = 1.3806488 \times 10 ^{-23} - \text{ Boltzmann's constant}$$ $$T = 350 \text{ K}$$ $$r = 70 \times 10^{-9} \text{ m}$$ $$\rho = (2/1000)*100^3 - \text{ Density of particle} = 2 gm/cc$$ $$k = (50*1000*2 \pi )^2 m = m \Omega^2 - \text{ Resonant Frequency = 50 KHz}$$ $$m = 4/3 \pi r^3 \rho - \text{ Mass of particle }$$ $$\Gamma = 2 \pi 10^{-2}$$

• 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. Commented Aug 16, 2013 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. Commented Aug 16, 2013 at 10:55
• But if the above same system follows fractional Brownian motion, Is it possible to write code in Mathematica ? Commented Feb 22, 2017 at 11:53

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 ? Commented Aug 16, 2013 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. Commented Aug 17, 2013 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)
*)


$$e^{-\frac{g t}{2 m}} \left(\frac{(g \text{x0}+2 m \text{y0}) \sinh \left(\frac{t \sqrt{g^2-4 k m}}{2 m}\right)}{\sqrt{g^2-4 k m}}+\text{x0} \cosh \left(\frac{t \sqrt{g^2-4 k m}}{2 m}\right)\right)$$ 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 Commented Aug 16, 2013 at 16:40