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 to 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 the 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 wants to show something, here is some sample 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}
              },
              r,
              {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
3  
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
1  
Stefan's answer to the schroedinger question is quite relevant –  gpap Aug 26 '13 at 8:50

3 Answers 3

up vote 8 down vote accepted

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

G = 4 Pi^2 // N;

stateData = 
 First[
    NDSolve`ProcessEquations[
        {  r''[t] == -G r[t]/Norm[r[t]]^3, 
           r[0] == {1, 0, 0},
           r'[0] == {0, 2 Pi // N, 0}
           },
        r,
        {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]];
  rAndDerivativesValues,
  {i, 1, 365, 7 (* 7 = each week*)}
   ] ;

rAndDerivativesValues is the value of r[t] and its derivative at each week.
Example : example

Here is one revolution of the Earth around the Sun :

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

weekly Earth positions

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

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

There is actually a very simple way to do this: instead of specifying a list of functions in the second argument, specify an empty list instead. This now begs the question of how one can obtain results. The solution is to use the event location functionality of Mathematica, by trapping the required events with Sow[]/Reap[].

In the current Mathematica, one will of course use WhenEvent[]. As I do not have a newer version on the machine I am currently borrowing, I will use the older Method setting "EventLocator"; translation to the modern form ought to be straightforward, however.

First, the easy case of returning only the value at the endpoint:

G = 4 Pi^2 // N; 
Reap[NDSolve[{r''[t] == -G r[t]/Norm[r[t]]^3,
              r[0] == {1, 0, 0}, r'[0] == {0, 2 Pi, 0}},
             {}, (* empty! *) {t, 1}, 
             Method -> {"EventLocator", "Event" -> (t == 1), 
                        "EventAction" :> Sow[r[t]],
                        Method -> "ExplicitRungeKutta"}, 
             MaxStepSize -> 1/365]][[-1, 1, 1]] // Quiet

yields {0.9999999995331317, 3.8260633123538*^-8, 0.} on this machine.

To borrow andre's example of sampling the Keplerian orbit at weekly intervals, we can do this:

weeks = Reap[NDSolve[{r''[t] == -G r[t]/Norm[r[t]]^3,
                      r[0] == {1, 0, 0}, r'[0] == {0, 2 Pi, 0}},
                     {}, {t, 1}, 
                     Method -> {"EventLocator",
                                "Event" -> {Sin[π t/(7./365)], t == 1}, 
                                "EventAction" :> Sow[r[t]],
                                Method -> "ExplicitRungeKutta"}, 
                     MaxStepSize -> 1/365]][[-1, 1]] // Quiet;

Graphics3D[Point[weeks]]

Keplerian orbit, sampled weekly

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.