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} $$

Answer

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]