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

2 Answers 2

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

There is actually a somewhat cryptic but very simple way to make NDSolve only return the solution at the end point:

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

As it is easy to oversee: the difference is that you set the third argument to {t,1,1}, NDSolve will in such a case still start to solve from t=0 as that is where the initial conditions are given but it only returns the solution at the end point. That idea is of course much more useful when solving e.g. partial differential equations. If you compare the ByteCount of that solution you will find that it is indeed much smaller than that of thefull solution. It probably is interesting to note that it is also possible to use something like {t,0.9,1} to only get the solution for 0.9<=t<=1.

If you need the solution at only some points, you can also do something like:

fullsol = Reap[NDSolveValue[{
    r''[t] == -G r[t]/Norm[r[t]]^3,
    r[0] == {1, 0, 0},
    r'[0] == {0, 2 Pi // N, 0},
    WhenEvent[Mod[10*t, 1] == 0, Sow[{t, r[t]}]]
  }, r, {t, 1, 1},
  Method -> "ExplicitRungeKutta", MaxStepSize -> (1/365 // N)


fullsol will then contain the accumulated solution points in a list which you can either convert to an interpolating function or use as it is for further investigations, e.g. visualizations:

Graphics3D[Point[fullsol[[2, 1, All, 2]]]]

You might want to test whether this or the NDSolve`Iterate approach does meet your needs better (it might well be that one is more performant than the other one for certain cases).

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.