I'm trying to do a fairly complex integration. The details aren't relevant, but are at https://github.com/barrycarter/bcapps/blob/master/ASTRO/bc-integrate-init.m and https://github.com/barrycarter/bcapps/blob/master/ASTRO/bc-integrate.m for reference.

When I run the NDSolve command, this happens:

sol = NDSolve[{posvel,accels},planets,{t,-366*500,366*500}, MaxSteps->100000, AccuracyGoal -> 50];

NDSolve::mxst: Maximum number of 100000 steps reached at the point t == -65522.6.

NDSolve::mxst: Maximum number of 100000 steps reached at the point t == 65285.8.

How do I tell Mathematica: treat the data at t == -65522.6 and t == 65285.8 as initial conditions and continue solving the equation for the next 65000+ time units, and so on, until you've reached the interval I originally requested: -366*500,366*500 (or -183000, +183000).

I realize the results won't be as accurate as if I'd just used more steps, but I'm running out of memory, so I want to NDSolve this piecewise. Any thoughts?

  • $\begingroup$ @MichaelE2 Well, OK, but I was hoping for a more general and automatic method for NDSolve in general, not just this specific example. $\endgroup$
    – user1722
    Commented Dec 15, 2016 at 1:22
  • $\begingroup$ Perhaps, it is possible to use the information in NDSolve Components and Data Structure to accomplish your goal. $\endgroup$
    – bbgodfrey
    Commented Dec 15, 2016 at 2:09
  • 1
    $\begingroup$ @MichaelE2 Are you assuming that the OP only wants the solution at the two endpoints? I had assumed that he wanted the solution throughout the region but wished to compute and store it in pieces. $\endgroup$
    – bbgodfrey
    Commented Dec 15, 2016 at 2:25
  • 1
    $\begingroup$ @MichaelE2 I think he wants the entire solution on {-366*500, 366*500} but in bite-sized pieces that he can store somewhere. After each piece is stored, I believe he wants to run a new computation with the final result of the previous piece as the initial condition for the next piece. I am confident that this can be done in a straightforward manner using the NDSolve Components but would be interested in knowing a better way. In fact, it can be done using NDSolve in the usual manner too, $\endgroup$
    – bbgodfrey
    Commented Dec 15, 2016 at 2:44
  • 1
    $\begingroup$ Have a look at this $\endgroup$
    – user21
    Commented Dec 15, 2016 at 17:45

1 Answer 1


As an alternative, let me suggest the "Extrapolation" method, with one of the implicit methods as a base integrator. It's a bit slower, but it seems to do a fairly accurate job, at least compared to the NASA data provided in the Github files. And it completes the integration in 28,000 to 45000 steps, a considerable savings over the Automatic method.

General procedure

First, I quit the kernel to get a measure of the memory usage. The OP's initialization file was put in an initialization cell (not shown), so that it would be automatically initialized for each trial.

Second, set the base integration method in basemethod.

Third, execute the last code block, which calls NDSolve, collects some data, and outputs a couple of tables, one of time and space usage and another comparing the result to the NASA data.

I had to increase PrecisionGoal a bit, to get an accurate result. One might push it up a little more, but soon one will also have to increase WorkingPrecision, too, which will slow down the integration even more.


basemethod = < method goes here >

nasa = {1.764338998100163*^-01, -3.499459824403553*^-01, -2.055213183324884*^-01};
mem = MaxMemoryUsed[];
time = First@ AbsoluteTiming[
    sol = NDSolve[{posvel, accels}, Array[planet, 9], {t, -366*500, 366*500},
       Method -> {"Extrapolation", Method -> basemethod}, 
       PrecisionGoal -> 10, MaxSteps -> 100000, AccuracyGoal -> 50];];
mem = MaxMemoryUsed[] - mem;
steps = planet[1]["Grid"] /. First[sol] // Length;
size = ByteCount[sol];
TableForm[{{time, mem, steps, size}},
 TableHeadings -> {None, {"time", "memory", "steps", "size"}}]
With[{n = nasa, p = planet[1][360*20] /. First[sol]},
 TableForm[Transpose@{n, p, (p - n)/n},
   TableHeadings -> {{"x", "y", "z"}, {"NASA", "sol", "rel. error"}}]


The "LinearlyImplicitMidpoint" method is the fastest.


basemethod = "LinearlyImplicitMidpoint"
(* run the above *)

Mathematica graphics


The "LinearlyImplicitModifiedMidpoint" method takes the fewest steps and therefore the least memory.


basemethod = "LinearlyImplicitModifiedMidpoint"
(* run the above *)

Mathematica graphics


Well, I tried this one. There's not much to recommend "LinearlyImplicitEuler" over the other two. Still, I'll share the results to save others the trouble.


basemethod = "LinearlyImplicitEuler";
(* run the above *)

Mathematica graphics

The Automatic method (OP's)

The Automatic method used 200,000 steps for only one third of the interval in just over 14 sec.; it also used over 900MB of memory.

  • $\begingroup$ Thanks! I'd considered running every combination of method and submethod to see which combination was fastest, but felt it would be too painful. $\endgroup$
    – user1722
    Commented Dec 16, 2016 at 21:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.