# How to tell Mathematica: continue NDSolve using values when you reach limit?

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?

• @MichaelE2 Well, OK, but I was hoping for a more general and automatic method for NDSolve in general, not just this specific example. – user1722 Dec 15 '16 at 1:22
• Perhaps, it is possible to use the information in NDSolve Components and Data Structure to accomplish your goal. – bbgodfrey Dec 15 '16 at 2:09
• @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. – bbgodfrey Dec 15 '16 at 2:25
• @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, – bbgodfrey Dec 15 '16 at 2:44
• Have a look at this – user21 Dec 15 '16 at 17:45

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.

Quit[]

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"}}]
]


### "LinearlyImplicitMidpoint"

The "LinearlyImplicitMidpoint" method is the fastest.

Quit[]

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


### "LinearlyImplicitModifiedMidpoint"

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

Quit[]

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


### "LinearlyImplicitEuler"

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.

Quit[]

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


### 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.

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