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.

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?

share|improve this question
    
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)? –  drN Apr 24 at 13:46
    
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. –  Karolis Apr 24 at 22:32
    
May I ask if this molecular dynamics or some such? What is the physical phenomena that you are trying to model? –  drN Apr 24 at 23:23
    
It is a simple 1D Molecular Dynamics model for free diffusion –  Karolis Apr 24 at 23:40
    
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! –  drN Apr 25 at 11:56

1 Answer 1

up vote 2 down vote accepted

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

share|improve this answer

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.