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.

Finally I started to play with differential equations in Mathematica.

And I have faced the problem, which seems too me so basic that I'm afraid this question is going to be closed soon. However, I've failed in finding solution here or in documentation center.

My question is: how to set that NDSolve will not save whole InterpolationFunction for the result?

I'm only interested in final coordinates or for example each 100th. Is there simple way to achieve that?

Anticipating questions:

  • I know I can do something like r /@ Range[0.1, 1, .1] /. sol at the end, but still, there is whole interpolating function in memory. I want to avoid it because my final goal is to do N-Body simulation where N is huge and I will run out of the memory quite fast. What is important to me is only the set of coordinates as far in the future as it can be, not intermediate values.

  • I can write something using Do or Nest but I want to avoid it since NDSolve allows us to implement different solving methods in handy way.

  • I saw WolframDemonstrations/CollidingGalaxies and it seems there is an explicit code with Do :-/

  • Another idea would be to put NDSolve into the loop but this seems to be not efficient. Could it be even compilable?

Just in case someone want to show something here is the sample of code to play with:

G = 4 Pi^2 // N ;

sol = NDSolve[{
               r''[t] == -G r[t]/Norm[r[t]]^3,
               r[0] == {1, 0, 0},
               r'[0] == {0, 2 Pi // N, 0}
              {t, 0, 1}, Method -> "ExplicitRungeKutta", MaxStepSize -> (1/365 // N)

 ParametricPlot3D[Evaluate[r[t] /. sol], {t, 0, 1}]

  (* Earth orbiting the Sun. Units: Year/AstronomicalUnit/SunMass 
     in order to express it simply*)
share|improve this question
I did something similar by periodically reinitializing the system with current solution. See tutorial/NDSolveStateData for details. That allowed me to solve for large tmax and get out solution for [tmax - chunksize, tmax] without running out of memory –  ssch Aug 25 '13 at 17:13
@ssch looks like something what I can use but let me read it. Thanks. –  Kuba Aug 25 '13 at 17:19
Stefan's answer to the schroedinger question is quite relevant –  gpap Aug 26 '13 at 8:50
@gpap you're right, thanks. I've seen this but I was not paying attention to method then. –  Kuba Aug 26 '13 at 8:54

1 Answer 1

up vote 7 down vote accepted

Here is a solution inspired from tutorial/NDSolveStateData (Mathematica 8)

G = 4 Pi^2 // N;

stateData = 
        {  r''[t] == -G r[t]/Norm[r[t]]^3, 
           r[0] == {1, 0, 0},
           r'[0] == {0, 2 Pi // N, 0}
        {t, 0, 1}, 
         Method -> "ExplicitRungeKutta",
         MaxStepSize -> (1/365 // N)]]

res = Table[
  NDSolve`Iterate[stateData, i/365];
  rAndDerivativesValues = NDSolve`ProcessSolutions[stateData, "Forward"];
  stateData = First @ NDSolve`Reinitialize[stateData, Equal @@@ Most[rAndDerivativesValues]];
  {i, 1, 365, 7 (* 7 = each week*)}
   ] ;

rAndDerivativesValues is the value of r[t] and his derivative at each week.
Example : enter image description here

Here is one revolution of the Earth around the Sun :

ListPointPlot3D[#[[1, 2]] & /@ res ]

enter image description here

share|improve this answer

Your Answer


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.