Using NDSolve`StateData without processing it into rule

Question

Problem

I need to solve a differential equation system. Nothing fancy here so far with the issue of too less memory on my machine. The interpolation functions I get from NDSolve eat up my 16GB memory very fast and the machine starts to create large swap files then. This should be avoided. I almost fear this question will be closed as duplicate - but even with thorough reading I didn't find an answer to my specific questions below. Before you vote for close please continue with reading - I won't ask "How can I do this", my issue is a bit more complex as it involves both understanding the NDSolve`StateData (Question 1) and getting information from it very fast (Question 2).

Ideas and Questions

Having a look at the 2D Schrödinger equation or the official MMA documentation seems promising as we are able to split NDSolve into parts of

• 1st) preparing with NDSolve`ProcessEquations,
• 2nd) iterating with NDSolve`Iterate,
• 3rd) processing with NDSolve`ProcessSolutions.

This way only the NDSolve`StateData has to be kept in memory. A first early question that comes to my mind now:

---

Question 1: How does all this work? I read the linked documentation but am not able to fully understand it. Let us assume we have an explicit, one-step integration scheme of fourth order (like the well-known Runge-Kutta-4 approach). In university I learned it is only needed to keep information about the latest step $i$ in order to compute the $i+1$st step. Measuring the size of the underlying NDSolve`StateData object seems to prove this. It grows only with the first step (i.e. the first usage of NDSolve`Iterate) and then maintains a constant size.

I would naively assume it is now only possible to get information about the data values of the functions and first derivatives for the $i+1$st step. This is only partly true.

NDSolve`ProcessSolutions[stateData, "Forward"]

gives me the numerical values in the $i+1st$ step as expected but

NDSolve`ProcessSolutions[stateData]

gives me full InterpolatingFunction objects that I up to now have not considered possible because this would only be possible if all values in between would be stored as well. The constant ByteCount seems to say something different. So: How is MMA able to build InterplatingFunctions over the complete domain without storing all data for all done steps?

---

Let us now assume the differential equation system was built and iterated. To get results I now have the possibility to process the solutions (see 3rd step above) to get rules (either for a specific step or interpolation functions as mentioned in Question 1). I could for example repeatedly use

h[ti]/.NDSolve`ProcessSolutions[stateData, "Forward"]

to get the vector $\vec h$ for the latest step $t_i$. The replacement in each step is very time-consuming so I looked for alternatives in the linked Wolfram documentation. It seems it is possible to get the underlying raw solution very fast by means of

stateData@"SolutionData".

Part[#,3]& of it would be the sought after numerical list in my case. Unfortunately, NDSolve sorts the coefficients different (h1,h1',h10,h10',…) than I do (h1,h2,h3,…). This leads me to the

---

Question 2: Is there a fast way to get the function values (derivates not needed) correctly sorted from the NDSolve-Framework other than sorting it myself in each step? If the answer is yes: Is it faster than time and again using ReplaceAll as shown above?

---

Answer 1

[Update notice: I discovered that MaxStepFraction is applied locally in NDSolve`Iterate[] to the length of the interval being "iterated," not to the global interval specified in NDSolve`ProcessEquations[]. Also added interpolation methods that do not copy/unpack the array of function values.]

---

(1) How does this all work?

Since "all" seems rather broad, I'll give my interpretation for ODEs, for which I am much more familiar with the interpolation schemes than for PDEs. I suspect that the cumulative solution is stored in a private variable (which I failed to find) and not in state data structure returned by NDSolve`ProcessEquations. The state stores only the data for the current step. But it's clear the total memory use of Mathematica grows linearly with the number of steps NDSolve takes. Furthermore, the growth rate is optimal for the solution being computed. By optimal, I mean if the setting is InterpolationOrder -> Automatic, then the rate is. For an ODE of order $n$, NDSolve stores by default the values $$x_k,\ y_k,\ y'_k, \dots,\ y^{(n)}_k$$ plus step indices for the Developer`PackedArrayForm of the Hermite InterpolatingFunction of the solution (see here and here). That means that $n+3$ 64-bit numbers are stored for each step, and memory usage can be observed to grow at $8\,(n+3)$ bytes per step.

Note that the data is not stored if the NDSolve call is of the form

NDSolve[ode, {}, {x, x1, x2}]

This indeed saves memory and returns no output. Presumably one uses this in conjunction with StepMonitor or EvaluationMonitor to construct other data. I tried using StepMonitor to construct a table of just the function values, but I wasn't really able to beat NDSolve.

(2) What is a fast way to get the function values?

The reason that NDSolve`ProcessSolutions[] is slower than NDSolve`SolutionDataComponent[] is that NDSolve`ProcessSolutions[] restructures the stored data for the InterpolatingFunction that it creates, which entails copying the data. NDSolve`SolutionDataComponent[] simply copies the current values stored in state.

Here is a way to extract the values of the variables without knowing the structure of the solution data (from the docs): This could be used to get function values at prescribed values of the independent variables without accumulating all the data. It may be sufficient for some use-cases.

vars = {y, z};
t1 = 0; t2 = 20; dt = 1./8;
{state} = NDSolve`ProcessEquations[{y'[t] == z[t], z'[t] == -y[t], y[0] == 1, z[0] == 0},
   {}, {t, t1, t2}, 
   MaxStepFraction -> 1];   (* see usage note below *)
vpos = Position[
   NDSolve`SolutionDataComponent[state@"VariableData", "X"], 
   Alternatives @@ vars, 1];
data = Table[
   NDSolve`Iterate[state, tt];
   Extract[NDSolve`SolutionDataComponent[state@"SolutionData"["Forward"], "X"], vpos],
   {tt, t1, t2, dt}];
data // Developer`PackedArrayQ
(* False *)

ListLinePlot[Transpose@data, DataRange -> {t1, t2}]

One drawback is that data is not packed, no matter how long the table. One could generate packed data as follows:

comp[tt_] := (NDSolve`Iterate[state, tt]; 
   Extract[NDSolve`SolutionDataComponent[state@"SolutionData"["Forward"], "X"], vpos]);
gendata = Compile[{t1, t2, dt},
   Table[comp[tt], {tt, t1, t2, dt}],
   {{comp[_], _Real, 1}}];
data = gendata[t1, t2, dt];
data // Developer`PackedArrayQ
(* True *)

The compiled function consists of a single MainEvaluate call to comp[tt], so the only benefit of compiling is that it generates a packed array for output. For the test case t1 = 0; t2 = 10000; dt = 1./16, the compiled version took 10% longer but the packed data was 1/8 the size.

If interpolation is needed, one can avoid copying data by applying InterpolatingPolynomial to just the data values needed, which is relatively easy to do since the values are over a uniform dt grid. For data as above (each row is a vector of values of different functions), the following will calculate the piecewise cubic interpolation at a real number t of the function with index f in the rows of data. On the example above, it is also much faster than Interpolation, unless you want to do more than a few thousand individual interpolations. Also note that if data is not a packed array, the pattern check ?MatrixQ will slow down evaluation considerably.

Clear[interp];
interp[data_?MatrixQ, f_Integer, t1_Real, t2_Real, dt_Real, t_] := 
 With[{idx = Clip[Floor[(t - t1)/dt], {1, Length@data}]},
  InterpolatingPolynomial[
   Transpose@{t1 + Range[idx - 1, idx + 2] dt, data[[idx ;; idx + 3, f]]}, t]
  ]

On packed data, using a compiled function is much faster:

interpC = Compile[
  {{data, _Real, 2}, {f, _Integer}, {t1, _Real}, {t2, _Real}, {dt, _Real}, {t, _Real}},
  Module[{idx, vals, dd = {0., 0., 0., 0.}, ft},
   idx = Floor[(t - t1)/dt];
   If[idx < 1, idx = 1];
   If[idx > Length@data - 3, idx = Length@data - 3];
   dd = data[[idx ;; idx + 3, f]];           (* Newton divided differences *)
   Do[
    dd[[j ;; 4]] = Differences[dd[[j - 1 ;; 4]]]/((j - 1)*dt),
    {j, 2, 4}];
   ft = dd[[4]];                             (* interpolate *)
   Do[
    ft = dd[[j]] + (t - t1 - (idx + j - 2) dt) ft,
    {j, 3, 1, -1}];
   ft
   ],
  RuntimeAttributes -> {Listable}, Parallelization -> True
  ]

NDSolve usage note: The default setting of MaxStepFraction is 1/10. When using NDSolve`Iterate[], this is applied to the length of the interval being computed by Iterate[], i.e., dt in the example above; it is not applied to the length of the interval {t, t1, t2} of the original call. For a small dt, the default setting might lead to as many as ten times the steps as necessary. Setting it to larger value will save time. For the case t1 = 0; t2 = 10000; dt = 1./16, the setting MaxStepFraction -> 1 takes 1/3 less time.

---

Example:
This is a second order equation, so it should use 8*(2 + 3) == 40 bytes per step (marginal rate).

Quit[]

(* new cell *)
times = Flatten@N@Outer[Times, 10^Range[-1, 5], {1, 2, 5}];
{state} = 
  NDSolve`ProcessEquations[
   {y''[x] == -y[x], y[0] == 1, y'[0] == 0},
   y,  (* change to {} to see no growth in memory usage *)
   {x, 0, 1000000}, MaxSteps -> Infinity, MaxStepSize -> 1];
data = Table[NDSolve`Iterate[state, x1]; {x1, MemoryInUse[]}, {x1, times}];

ListLinePlot[data, PlotRange -> All]

Getting the bytes per step rates:

sol = NDSolve`ProcessSolutions[state];
grid = y["Grid"] /. sol // Flatten;
steps = LengthWhile[grid, x \[Function] x <= #] & /@ times
(* steps at the end of each NDSolve`Iterate:
  {15, 26, 40, 55, 69, 95, 137, 216, 448, 829, 1587, 3856, 7627, 15157, 
   37771, 75450, 150832, 377150, 754497, 1509044, 3743165}
*)

Differences@data[[All, 2]] / Differences@steps // N
(* Bytes per step, for each iteration interval
  {61.0909, 6.28571, 44.8, 94.2857, 26.4615, 45.7143, 39.6962, 34.8276, 
   38.9921, 39.905, 38.8365, 39.1811, 39.0704, 39.0806, 39.0522, 
   39.0707, 39.081, 39.0805, 39.077, 39.0725}
*)

(Note that due to auxiliary code using and releasing memory, exact results rarely occur.)

One can reduce the per step memory usage by reducing the order of the system:

{state} = NDSolve`ProcessEquations[
   {y'[x] == z[x], z'[x] == -y[x], y[0] == 1, z[0] == 0},
   y,  (* keep only data on  y, discard  z  data *)
   {x, 0, 1000000}, MaxSteps -> Infinity, MaxStepSize -> 1];

This has roughly a 32 byte per step marginal rate of memory usage (8 * (1 + 3) == 32).

Answer 2

Personal remark: I never answered my own question up to now and usually refrain from doing so for good reason as it looks like simple reputation increase. I therefore won't accept my own answer, maybe someone comes around with something more specific whose answer will be accepted then by me.

---

Possible fast solution: Assume sys is a DES that was built by NDSolve`ProcessEquations. Let h[x] be the sought after coefficient vector with the dependent variables of the DES in it.

First find the structure that maps NDSolve's internal data format to our coefficient vector by using

DESDataOrder[sourceList_, targetList_] := {Ordering@sourceList, Ordering@Ordering@targetList};
hOrder = DESDataOrder[NDSolve`SolutionDataComponent[sys@"Variables", "DependentVariables"], Head /@ h[x]];

and save the variable hOrder somewhere. Let's compare two approaches now:

SetAttributes[{DESStep, DESStepData}, HoldFirst];

First approach:

DESStep[stateData_, t_] := (NDSolve`Iterate[stateData, t]; NDSolve`ProcessSolutions[stateData, "Forward"]);

Second approach:

DESStepData[stateData_, t_] := (NDSolve`Iterate[stateData, t]; NDSolve`SolutionDataComponent[(stateData@"SolutionData")[[2]], "DependentVariables"][[hOrder[[1]]]][[hOrder[[2]]]]);

Comparison:

Note that DESStep builds replacement rules whereas DESStepData just sorts the underlying data from NDSolve in such a form that we can get the numerical coefficient vector directly. Due to the not needed replacement via rules the second approach is noticeably faster:

numerical = DESStepData[sys, 1]; // AbsoluteTiming
(* Out: {0.000052, Null}*)
ruleBased = h[1.] /. DESStep[sys, 1]; // AbsoluteTiming
(* Out: {0.04334, Null}*)

Caution: This approach heavily relies on knowledge about the underlying data structure. This can (but hopefully won't) be changed by WRI.