8
$\begingroup$

NDSolve can be broken into three stages:

  1. NDSolve`ProcessEquations processes the equations and sets up an NDSolve`StateData object
  2. NDSolve`Iterate iterates the differential equations
  3. NDSolve`ProcessSolutions processes the solutions into InterpolatingFunctions

(see also this answer by @xzczd).

What is inside an NDSolve`StateData object? Can we create our own valid NDSolve`StateData object to bypass NDSolve`ProcessEquations? Can we modify an existing NDSolve`StateData object?

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

$\endgroup$
6
$\begingroup$

This is a partial answer to the first two questions (what is inside an NDSolve`StateData object? can we create our own valid NDSolve`StateData object to bypass NDSolve`ProcessEquations?). 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 NDSolve`StateData object created by NDSolve`ProcessEquations to reverse engineer it? This is apparently version dependent. In versions 10.3 and 11.2, we can just take parts of an NDSolve`StateData object:

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

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 NDSolve`StateData objects in these later versions, so this is only an issue when trying to reverse-engineer the internals of NDSolve`StateData.

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

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

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

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 NDSolve`StateData has eleven parts:

Length[s]
(* 11 *)

Taking a look at them:

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

Mathematica graphics

It's kind of tedious, but by using a few well-chosen calls to NDSolve`ProcessEquations 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 NDSolve`StateData 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..."]];
  If[VectorQ[vars, Head[#] == Symbol &],
    xvars = vars;
    toxvars = fromxvars = {}
  ,
    If[debug, Print["making xvars..."]]; 
    xvars = Table[Unique[TemporaryVariable], neq];
    If[debug, Print["making toxvars..."]]; 
    toxvars = Dispatch[Thread[vars -> xvars]];
    If[debug, Print["making fromxvars..."]]; 
    fromxvars = Dispatch[Thread[xvars -> vars]];
  ];

  (* 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[NDSolve`xs], 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, {NDSolve`ProcessEquations, None,
    NDSolve`ProcessEquations, NDSolve`ProcessEquations}};


  (* part[2] -- NDSolve`ProcessEquations 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] -- Experimental`NumericalFunction 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"},
    NDSolve`ProcessEquations, Join[{t}, Table[var[t], {var, vars}]], None};

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

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

  mon = Unique[NDSolve`Monitor];
  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, Internal`InheritedBlock[v, CompoundExpression[body]]]]]]
      &[mon, mons, Prepend[vars, t]];

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

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

  part[3] = Experimental`NumericalFunction[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"]},
    Internal`WithLocalSettings[
      SetSystemOptions["CompileOptions" -> {"TableCompileLength" -> \[Infinity]}], 
    part[6, 6] = {Internal`Bag[t0], {}, Table[Internal`Bag[{ics[[i]], ids[[i]]}], {i, neq}],
      {}, {}, {}, {}, {}, {}},
    SetSystemOptions[tcl]]
  ];

  part[6, 7] = {{}, Table[Internal`Bag[], {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, {NDSolve`ScaledVectorNorm[2, {1.0536712127723497`*^-8, 1.0536712127723497`*^-8},
    NDSolve`ProcessEquations], {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" -> NDSolve`ProcessEquations}]], 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[NDSolve`StateData, HoldAllComplete];
  res = NDSolve`StateData[Sequence @@ parts];
  SetAttributes[NDSolve`StateData, 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]
(* NDSolve`StateData["<" 0. ">"] *)
NDSolve`Iterate[s, 1]
sol = NDSolve`ProcessSolutions[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 NDSolve`ProcessEquations takes 1.1. That's the overhead saved by dealing with only one kind of system.

A few notes:

  • the Experimental`NumericalFunction in part[3] doesn't seem to have the same format as one made by Experimental`CreateNumericalFunction 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 NDSolve`StateData 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

$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – user21 Jul 30 at 13:35
  • $\begingroup$ 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. $\endgroup$ – user21 Jul 30 at 13:36
  • 1
    $\begingroup$ 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. $\endgroup$ – user21 Jul 30 at 13:37
  • $\begingroup$ @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. $\endgroup$ – Chris K Jul 30 at 13:38
  • $\begingroup$ 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. $\endgroup$ – user21 Jul 30 at 13:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.