Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
1  
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
1  
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

2 Answers 2

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]

enter image description here

share|improve this answer
    
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]

plot

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)) *)
share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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