5
$\begingroup$

I am using NDSolve for a Langevin dynamics problem. I want to the know long term behaviour of my system ($t>1$) but it has to be simulated with very small time steps ($dt\sim 10^{-9}$). An example code:

R[t_Real]:= RandomVariate[NormalDistribution[0,1]]

NDSolve[
  {(-x''[t] - k1*x'[t] + k2*R[t])==0, x[0]==0, x'[0]==0} //.values,
  x,
  {t, 0, 10},
  StartingStepSize-> 10^-9,
  Method->{"FixedStep",Method->"ExplicitEuler"},
  MaxSteps->\[Infinity]
]

The Problem: My computer runs out memory when trying to store $10^{10}$ data points necessary for this computation. Is there a way to sample and store only a small subset of all the integration points?

$\endgroup$
  • $\begingroup$ Have you tried [Mma's RandomVariate[...]](reference.wolfram.com/mathematica/ref/RandomVariate.html)? What does values do? Is there a reason you are using ExplicitEuler? You would benefit by letting Mathematica choose the solver automatically or use perhaps the magic-wand Method->"LSODA". Also, why not just let t range from 0 to a very large number (for long term behavior)? $\endgroup$ – dearN Apr 24 '14 at 13:46
  • $\begingroup$ Hi, the 'values' variable is used to store parameters k1 and k2. I do not trust the automatic solvers of Mathematica because 'R[t]' is not a continuos function. For me $t \sim 1$ is a long time behaviour. The problem presented here has a damping of momentum (second term) that happens on the timescale of $10^{-9}$. The momentum has to relax from thermal kicks obtained from 'R[t]' term. For this reason I am using constant time sampling. $\endgroup$ – Karolis Apr 24 '14 at 22:32
  • $\begingroup$ May I ask if this molecular dynamics or some such? What is the physical phenomena that you are trying to model? $\endgroup$ – dearN Apr 24 '14 at 23:23
  • $\begingroup$ It is a simple 1D Molecular Dynamics model for free diffusion $\endgroup$ – Karolis Apr 24 '14 at 23:40
  • $\begingroup$ I performed simulations on MD for my Master's degree and I am curious, if you don't mind. Are you simulating Argon? With periodic BCs? Is the eqn you are solving the convection-diffusion equation? In 2012 I had a similar question on SE and I am wondering if it would help you. Good luck! $\endgroup$ – dearN Apr 25 '14 at 11:56
4
$\begingroup$

You can have NDSolve do the integration but not save the results (set second argument to {}). Then use EvaluationMonitor to save points at whatever interval dt you please. This uses much less memory; in fact for dt = 10^-2, the memory use is negligible for the settings below.

Since parameters were not given by the OP, I used the simplest choices. Also, waiting for ten billion steps seemed silly for a proof-of-concept trial, so I lengthened the step size.

MaxMemoryUsed[]

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,
    {}, {t, 0, 10}, 
    StartingStepSize -> 10^-5, 
    Method -> {"FixedStep", Method -> "ExplicitEuler"}, MaxSteps -> ∞,
    EvaluationMonitor :> If[t >= lastt + dt, lastt = t; Sow[{t, x[t]}]]
    ];
sol = Interpolation[pts];

MaxMemoryUsed[]

Plot[sol[t], {t, 0, 10}]

(* 37243360 *)

(* 37243360 *)

Mathematica graphics

When the integral x is requested, the memory use by NDSolve is significantly higher:

MaxMemoryUsed[]

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]}]]
    ];
sol = Interpolation[pts];

MaxMemoryUsed[]

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

(* 37243360 *)

(* 127874624 *)

Mathematica graphics

$\endgroup$

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.