Can we construct our own NDSolveStateData?

NDSolve can be broken into three stages:

1. NDSolveProcessEquations processes the equations and sets up an NDSolveStateData object
2. NDSolveIterate iterates the differential equations
3. NDSolveProcessSolutions processes the solutions into InterpolatingFunctions

What is inside an NDSolveStateData object? Can we create our own valid NDSolveStateData object to bypass NDSolveProcessEquations? Can we modify an existing NDSolveStateData object?

Knowing the answer to these fundamental questions might help address other questions such as these:

This is a partial answer to the first two questions (what is inside an NDSolveStateData object? can we create our own valid NDSolveStateData object to bypass NDSolveProcessEquations?). It's only a partial answer, because NDSolve has different modes for different kinds of problems (ordinary differential equations vs differential-algebraic equations vs partial differential equations). Hopefully others will add answers that address these other modes.

First, how can we look inside an NDSolveStateData object created by NDSolveProcessEquations to reverse engineer it? This is apparently version dependent. In versions 10.3 and 11.2, we can just take parts of an NDSolveStateData object:

s = NDSolveProcessEquations[{x'[t] == 13 x[t], x[0] == 73}, x, t][[1]]
s[[1]]
(* NDSolveStateData["<" 0. ">"] *)
(* {5, 256, {NDSolveProcessEquations, None, NDSolveProcessEquations,
NDSolveProcessEquations}} *)


Unfortunately this fails in versions 11.3 and 12.0. If you know a way around this, please comment. However, we can still construct valid NDSolveStateData objects in these later versions, so this is only an issue when trying to reverse-engineer the internals of NDSolveStateData.

Changing the Method->{EquationSimplification} option alters s[[1, 2]]:

s = NDSolveProcessEquations[{x'[t] == 13 x[t], x[0] == 73}, x, t,
Method -> {EquationSimplification -> MassMatrix}][[1]];
s[[1]]
(* {5, 257, {NDSolveProcessEquations, None, NDSolveProcessEquations,
NDSolveProcessEquations}} *)

s = NDSolveProcessEquations[{x'[t] == 13 x[t], x[0] == 73}, {x}, t,
Method -> {EquationSimplification -> Residual}][[1]];
s[[1]]
(* {5, 258, {NDSolveProcessEquations, None, NDSolveProcessEquations,
NDSolveProcessEquations}} *)


Evidently s[[1, 2]] == 256 corresponds to ODEs and s[[1, 2]] == 257 and s[[1, 2]] == 258 to two different methods for solving DAEs. I'm sure other modes exist for PDEs and who knows what else. For this answer, I'll focus only on systems of first-order ODEs with s[[1, 2]] == 256.

Returning to my first example, we see that NDSolveStateData has eleven parts:

Length[s]
(* 11 *)


Taking a look at them:

Do[Print[i,":"]; Print[s[[i]]], {i, 11}]


It's kind of tedious, but by using a few well-chosen calls to NDSolveProcessEquations as probes, we can figure out what goes where. The number of equations is a common element, as are the dependent variables, right-hand sides, initial conditions and initial derivatives.

Feynmann wrote, "what I cannot create, I do not understand." Without claiming to actually understand all of these internal parts, perhaps the easiest way to describe them is to write a function to create our own mode==256 NDSolveStateData object (no WhenEvents, no ParametricSensitivity, just first-order ODEs).

ProcessFirstOrderODEs[vars_List, rhs_List, icsin_List, t0in_?NumericQ,
opts___?OptionQ] := Block[{jacobian, neq, xvars, toxvars, fromxvars, uvars, uxss,
t0, ics, ids, part, parts, mon, mons, str, res},

jacobian = Evaluate[Jacobian /. Flatten[{opts, Options[ProcessFirstOrderODEs]}]];

If[debug, Print["calculating neq..."]];
neq = Length[vars]; (* # of eqns *)

(* if there are any non-Symbol vars, make TemporaryVariables in xvars
and Dispatches to convert *)

If[debug, Print["checking vars for non-Symbols..."]];
xvars = vars;
toxvars = fromxvars = {}
,
If[debug, Print["making xvars..."]];
xvars = Table[Unique[TemporaryVariable], neq];
If[debug, Print["making toxvars..."]];
If[debug, Print["making fromxvars..."]];
];

(* add $number to vars to stand in for derivatives in Functions *) If[debug, Print["making uvars..."]]; uvars = Unique[xvars]; If[debug, Print["making uxss..."]]; uxss = Table[Unique[NDSolvexs], neq]; If[debug, Print["making t0..."]]; t0 = N[t0in]; (* initial time *) If[debug, Print["making ics..."]]; ics = N[icsin]; (* initial conditions *) (* part[1] -- ?? part[1,2] = Mode (256=first-order ODEs) *) If[debug, Print["part[1]..."]]; part[1] = {5, 256, {NDSolveProcessEquations, None, NDSolveProcessEquations, NDSolveProcessEquations}}; (* part[2] -- NDSolveProcessEquations Options? *) If[debug, Print["part[2]..."]]; part[2] = {"TimeIntegration" :> Automatic, "BoundaryValues" :> Automatic, "DiscontinuityProcessing" :> Automatic, "EquationSimplification" :> Automatic, "IndexReduction" :> None, "DAEInitialization" :> Automatic, "PDEDiscretization" :> Automatic, "ParametricCaching" :> Automatic, "ParametricSensitivity" :> Automatic}; (* part[3] -- ExperimentalNumericalFunction with RHS *) If[debug, Print["part[3,1]..."]]; part[3, 1] = {Function[Evaluate[Join[{t}, xvars]], Evaluate[rhs /. toxvars]], Apply}; If[debug, Print["part[3,2]..."]]; part[3, 2] = {0, Join[{{{}, 1, 0, 0, 0, 0}}, Table[{{}, 2, i - 1, 0, 0, 0}, {i, neq}]]}; If[debug, Print["part[3,3]..."]]; part[3, 3] = {{{1, 1, 818}, {{}, {}}}, {{3, neq, 817}, {{jacobian, Automatic, None, 1, Automatic}}}}; If[debug, Print["part[3,4]..."]]; part[3, 4] = {0, 3, {neq}, 0}; If[debug, Print["part[3,5]..."]]; part[3, 5] = {8236, MachinePrecision, {{Automatic}, Automatic}, True, {{Automatic, "CleanUpRegisters" -> False, "WarningMessages" -> False, "EvaluateSymbolically" -> False, "RuntimeErrorHandler" -> ($Failed &)}, {}, Automatic, "WVM"},
NDSolveProcessEquations, Join[{t}, Table[var[t], {var, vars}]], None};

If[debug, Print["part[3,6]..."]];

(* by @MichaelE2 <https://mathematica.stackexchange.com/a/
202891> *)

mon = Unique[NDSolveMonitor];
mons = Table[Unique[mon], {neq + 1}];

part[3, 6, 1] = With[{code =
Join[Hold[{#1}, #2, #3],(*first args of Function and
InheritedBlock*)
Unset /@ Hold @@ #3,(*beginning of body*)
Set @@@ Hold @@ Transpose@{Prepend[Through[Rest[#3][First[#3]]],
First[#3]], #2}, Hold[#1]]},
Replace[code,
Hold[m1_, m2_, v_, body__] :>
Function[m1, Function[m2, InternalInheritedBlock[v, CompoundExpression[body]]]]]]
&[mon, mons, Prepend[vars, t]];

part[3, 6] = {part[3, 6, 1], None, None};

(*part[3,6]={#&,None,None};*)

part[3] = ExperimentalNumericalFunction[part[3, 1], part[3, 2], part[3, 3],
part[3, 4], part[3, 5], part[3, 6]];

(* part[4] -- ?? *)

If[debug, Print["part[4]..."]];
part[4, 1] = {{neq, 1, 0, neq, 0, 0, 0, 0, 0}, {0, 1, 1, neq + 1,
neq + 1, neq + 1, neq + 1, neq + 1, neq + 1}};

part[4, 2] = {0, {#1 /. toxvars &, #1 &, #1 /. fromxvars &},
{1, {t}}, {xvars, xvars, vars}};

part[4, 3] = part[4, 4] = None;

part[4, 5, 1] = {0, 1, 1, neq + 1, neq + 1, neq + 1, neq + 1, neq + 1, neq + 1};
part[4, 5, 2] = {0, Join[{{{}, 1, 0, 0, 0, 0}},
Table[{{}, 2, i - 1, 0, 0, 0}, {i, neq}]]};
part[4, 5, 3] = Function[Evaluate[Join[{t}, xvars, uvars]],
Evaluate[{t, {}, xvars, uvars, {}, {}, {}, {}}]];
part[4, 5] = Table[part[4, 5, i], {i, 3}];

part[4, 6] = Table[{var, var'}, {var, vars}];

part[4] = Table[part[4, i], {i, 6}];

(* part[5] -- Initial Conditions *)

If[debug, Print["making ids..."]];
ids = part[3][0, ics];

If[debug, Print["part[5]..."]];
part[5, 2] = {{t0, None, ics, ids, {}, {}, {}, {}}, 0, Automatic, None, None, True};
part[5] = {None, part[5, 2], None};

(* part[6] -- Results Store *)

If[debug, Print["part[6]..."]];
part[6, 2] = {neq, 1, 0, neq, 0, 0, 0, 0, 0};

part[6, 3] = Function[Evaluate[uxss], Evaluate[Thread[vars -> uxss]]];

part[6, 5] = {Range[neq], Table[1, neq], Table[0, neq],
{Table[0, 9], {}}, {{0, 0, 0, neq, neq, neq, neq, neq, neq},
Range[0, neq - 1]}, Range[neq]};

(* see <https://mathematica.stackexchange.com/questions/202869/> *)

With[{tcl = SystemOptions["CompileOptions" -> "TableCompileLength"]},
InternalWithLocalSettings[
SetSystemOptions["CompileOptions" -> {"TableCompileLength" -> \[Infinity]}],
part[6, 6] = {InternalBag[t0], {}, Table[InternalBag[{ics[[i]], ids[[i]]}], {i, neq}],
{}, {}, {}, {}, {}, {}},
SetSystemOptions[tcl]]
];

part[6, 7] = {{}, Table[InternalBag[], {4}]};

part[6] = {1, part[6, 2], part[6, 3], Automatic, part[6, 5], part[6, 6], part[6, 7]};

(* part[7] -- Options *)

If[debug, Print["part[7]..."]];
part[7] = {0, Automatic, {NDSolveScaledVectorNorm[2, {1.0536712127723497*^-8, 1.0536712127723497*^-8},
NDSolveProcessEquations], {Automatic, \[Infinity], 1/10}, t},
{Automatic, Automatic,

(* merge opts and default opts - <https://
mathematica.stackexchange.com/a/135242/> *)

GatherBy[
Flatten[Join[{opts}, {AccuracyGoal -> Automatic,  PrecisionGoal -> Automatic,
WorkingPrecision -> MachinePrecision, InterpolationPrecision -> Automatic,
Compiled -> Automatic, Jacobian -> Automatic,
Method -> {"TimeIntegration" :> Automatic,  "BoundaryValues" :> Automatic,
"DiscontinuityProcessing" :> Automatic,
"EquationSimplification" :> Automatic,
"IndexReduction" :> None,
"DAEInitialization" :> Automatic,
"PDEDiscretization" :> Automatic,
"ParametricCaching" :> Automatic,
"ParametricSensitivity" :> Automatic},
"StoppingTest" -> None, "Events" -> None,
InterpolationOrder -> Automatic, MaxSteps -> Automatic,
StartingStepSize -> Automatic, MaxStepSize -> \[Infinity],
MaxStepFraction -> 1/10, "MaxRelativeStepSize" -> 1/10,
NormFunction -> Automatic, DependentVariables -> Automatic,
DiscreteVariables -> {}, SolveDelayed -> Automatic,
"CompensatedSummation" -> Automatic,
EvaluationMonitor -> None, StepMonitor -> None,
"MethodMonitor" -> None, "ExtrapolationHandler" -> Automatic,
"MinSamplingPeriod" -> Automatic,
"Caller" -> NDSolveProcessEquations}]], First][[All, 1]]
}, None, None, None};

(* part[8] -- Initial Conditions *)

If[debug, Print["part[8]..."]];
part[8] = {{0, 0}, Thread[xvars == icsin], {}, All, {}};

(* parts[9-11] -- Nothing *)

If[debug, Print["parts[9-11]..."]];
part[9] = part[10] = part[11] = {};

(* put together *)
parts = Table[part[i], {i, 11}];

(*Do[Print["part ",i]; Print[part[i]], {i,11}];*)

If[debug, Print["res..."]];
ClearAttributes[NDSolveStateData, HoldAllComplete];
res = NDSolveStateData[Sequence @@ parts];
SetAttributes[NDSolveStateData, HoldAllComplete];

Return[res]

];

Options[ProcessFirstOrderODEs] = {Jacobian -> Automatic};


Hope there aren't too many transcription errors there!

In use:

s = ProcessFirstOrderODEs[{x}, {13 x}, {73}, 0]
(* NDSolveStateData["<" 0. ">"] *)
NDSolveIterate[s, 1]
sol = NDSolveProcessSolutions[s]
(* {x->InterpolatingFunction[Domain: {{0.,1.}}
Output: scalar]} *)


Multiple equations:

s = ProcessFirstOrderODEs[{x, y, z}, {13 x, 17 y, 19 x}, {73, 89, 101}, 0];


Indexed equations:

nmax = 10000;
vars = Table[p[i], {i, nmax}];
rhs = Table[p[i] (1 - p[i]/i), {i, nmax}];
ics = ConstantArray[1, nmax];
s = ProcessFirstOrderODEs[vars, rhs, ics, 0];


RepeatedTiming of the last one is 0.417 second, where the equivalent NDSolveProcessEquations takes 1.1. That's the overhead saved by dealing with only one kind of system.

A few notes:

• the ExperimentalNumericalFunction in part[3] doesn't seem to have the same format as one made by ExperimentalCreateNumericalFunction as described here, so it had to be made manually
• not so confident about my Option handling in part[7]
• using indexed variables like p[1], p[2] incurs a cost because they need to be changing into TemporaryVariable\$num in the NumericalFunction, then changed back at the end.

In general, there are probably many ways this code could be improved, which I hope you all will provide. My actual problem that initiated this investigation deep into the internals of NDSolveStateData remains unsolved, but at least there's some hope for improvement still!

edit 7/31/19 - now calculate initial derivatives with part[3]'s NumericalFunction

edit 8/1/19 - added Jacobian option to pass to NumericalFunction

• Concerning the change of NDSolve`State in version 11.3. Internally that object has a different representation now (more like a SparseArray or Graph and not like a normal expression). This has the advantage that it can be passed around more efficiently. That being said, it's not a bug and will not change back to the pre. 11.3 behavior. Jul 30, 2019 at 13:35
• Can you remind me, what exactly was the initial problem you wanted to solve. Perhaps I can squeeze out some time to look at that. Jul 30, 2019 at 13:36
• One further point: Such an interface will be very brittle; once the internal of that data object need to be changed for whatever reason this interface will no longer work. Not to say not to do it, but to keep this in mind. Jul 30, 2019 at 13:37
• @user21 Thanks for the offer! I haven't actually posted my original problem, because it's complex and I thought I might figure it out on my own. But I'll tag you if I give up and post it. Jul 30, 2019 at 13:38
• Sure. I am going to go on a looooong (well not that long) vacation soon, it will have to wait then until I get back. Jul 30, 2019 at 13:40