1
$\begingroup$

I want to simulate the movement of a damped single particle which vibrates due to Brownian motion.

$m \ddot{x} + \gamma \dot{x} - \xi(t)=0$

where $\gamma$ is the friction constant and $\xi(t)$ a "random" or "fluctuating force".

In principle I would like to solve a similar equation as in here but I would like to get out the (x,y) coordinate pairs and how they vary with time.

So, how can I extend the 1d Langevin equation so that NDSolve produces the 2d coordinates?

This is the answer from Michael E2 for the 1d problem:

SeedRandom[1];
R[t_Real] := RandomVariate[NormalDistribution[0, 1]];
values = {k1 -> 1, k2 -> 1};

lastt = -1;
dt = 10^-2;

{foo, {pts}} = 
  Reap@NDSolve[{(-x''[t] - k1*x'[t] + k2*R[t]) == 0, x[0] == 0, x'[0] == 0} /. values,
    x, {t, 0, 10}, 
    StartingStepSize -> 10^-5, 
    Method -> {"FixedStep", Method -> "ExplicitEuler"}, MaxSteps -> ∞,
    EvaluationMonitor :> If[t >= lastt + dt, lastt = t; Sow[{t, x[t]}]]
    ];

Plot[x[t] /. First[foo] // Evaluate, {t, 0, 10}]

Mathematica graphics

$\endgroup$
2
$\begingroup$

This extension of your answer seems to work. I have taken out the evaluation monitor since it is not used in the plot.

SeedRandom[1];
R[t_Real] := RandomVariate[NormalDistribution[0, 1]];
values = {k1x -> 1, k2x -> 1,k1y -> 1, k2y -> 1};

foo=NDSolve[{
    (-x''[t] - k1x*x'[t] + k2x*R[t]) == 0, x[0] == 0, x'[0] == 0,
    (-y''[t] - k1y*y'[t] + k2y*R[t]) == 0, y[0] == 0, y'[0] == 1} /.values,
    {x[t],y[t]}, {t, 0, 10}, 
    StartingStepSize -> 10^-3, 
    Method -> {"FixedStep", Method -> "ExplicitEuler"}, MaxSteps -> Infinity
];

(* plot in x-y plane  *)
ParametricPlot[Evaluate[{x[t],y[t]}/.foo[[1]]],{t,0,10}]
(* plot with time on z axis  *)
ParametricPlot3D[Evaluate[{x[t],y[t],t}/.foo[[1]]],{t,0,10},BoxRatios->{1,1,1}]
$\endgroup$
  • $\begingroup$ Why is the effective displacement much larger in x than in y? That should not be in average. The movement is somehow directed, what is unclear to me. $\endgroup$ – mrz Mar 9 '18 at 13:41
  • 1
    $\begingroup$ @Lenoil Actually the displacement is much larger in y, that's because ha loa chooses x'[0] == 0 and y'[0] == 1 as the i.c.. And I believe $F$ should also be a vector if you need undirected result. $\endgroup$ – xzczd Mar 26 '18 at 14:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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