ParametricNDSolve fails when I *don't* print the results

Question

Some updates added to end. Truly minimal exampled added at top.

Minimal example:

tstfn[debug_] := Module[{tst}, 
  tst = ParametricNDSolve[{a'[t] == -k a[t], b'[t] == k a[t], 
                           a[0] == 1, b[0] == 0}, 
                          {a[t], b[t]}, {t, 0, 1}, {k}]; 
  If[debug, Print[tst]]; 
  (# -> ((# /. tst) @@ {k})) & /@ {a[t], b[t]}]

DynamicModule[{tst2, dbg}, 
  tst2 = Null; 
  dbg = False;
  Column[{Row[{"Include Print? ", Checkbox[Dynamic[dbg]]}], 
          Button["Recalc", tst2 = tstfn[dbg]], 
          Dynamic[Plot[Evaluate[({a[t], b[t]} /. tst2) /. 
                       {k -> 5, a[0] -> 1}], {t, 0, 1}]]}]]

The problem seems to be with the a[0]->1 substitution, even though it isn't really being used. I am guessing that this is for some reason messing up the internal representation of the ParametricFunction. I am still baffled as to why Printing the solution fixes the problem, but at least having narrowed down where the problem is, I think I can find another way to accomplish what I am trying to accomplish. I am leaving this question open in case someone has suggestions as to what is going on with the Print.

Original Post: Ok, I am thoroughly baffled. I have a function I have written that contains a ParametricNDSolve. Under certain conditions, it fails. Adding a Print statement (for debugging purposes) causes it to succeed. I have tried to narrow down my code to a much smaller code set; I hope it is not still too complicated to qualify as a minimal example.

The main idea is that I have a set of differential equations, and want to be able to interactively choose which parameters (constants in the equations and/or initial conditions) should be the variable parameters within ParametricNDSolve.

Here is a quick demonstration of the problem. First, when I call ParametricNDSolve with no variable parameters (i.e., the list of parameters is {}):

When I put one of the parameters in as variable (in this case, k1), it fails:

But when I Print the solution before leaving the function, it succeeds!

Does anyone have any idea why it is behaving this way?

Here is the routine I have put together to do the calculation (a walkthrough will follow):

Kineticstest[rateEqs_, rateConstants_, conc_, {t_, t0_, tf_}, 
             rateConstFixed_, initConcFixed_] := 
  Module[{solution, initialConditions, rateConstReplacements, tDepVars,
          variableParameters, internalSubstitutions}, 
    internalSubstitutions = MapIndexed[(#1[0] -> ic[#2[[1]]]) &, 
      Pick[conc[[All, 1]], ! # & /@ initConcFixed]]; 
    initialConditions = (#1[0] == ((#1[0] /. internalSubstitutions) /. 
      (#1[0] -> #2))) & @@@ conc; 
    rateConstReplacements = Pick[(#1 -> #2) & @@@ rateConstants, 
      rateConstFixed]; 
    tDepVars = #1[t] & @@@ conc; 
    variableParameters = Join[
      Pick[rateConstants[[All, 1]], ! # & /@ rateConstFixed], 
      Values[internalSubstitutions]]; 
    solution = ParametricNDSolve[
      Join[rateEqs, initialConditions] /. rateConstReplacements, 
      tDepVars, {t, t0, tf}, variableParameters];
    (*Print["solution = ", solution];*)
    (# -> ((# /. solution) @@ (variableParameters /. 
          (internalSubstitutions /. (a_ -> b_) -> (b -> a))))) & 
           /@ tDepVars]

So the idea is to pass:

• a list of differential rate equations (rateEqs),
• a list of rate constants like {{k1,5},{k2,5}} where the first entry in each pair is the rate constant symbol, and the second is its value (rateConstants),
• a list of concentrations like {{a,1},{b,0}} where the first entry in each pair is the symbol for the concentration (a[t] for the time-dependent form or a[0] for an initial value) and the second is the initial value (conc),
• the time variable (t) and the initial and final times (t0 and tf), and
• optionally a list of True/False values for each rate constant and each initial concentration to determined whether they should be considered fixed or variable.

I then do some manipulations to get this set of input information into the right form for ParametricNDSolve. We have:

1. the list of equations: Join[rateEqs,initialConditions]/.rateConstReplacements. The rate equations are left unchanged. The initial conditions are either specific value assignments (e.g., a[0]==1) or an assignment to an internal variable if the initial condition is one of the variable parameters (e.g., a[0]==ic2; here ic2 is the first initial concentration); the assignment of these internal variables is handled by the internalSubstitutions local variable. The rateConstReplacements are only done for rate constants that are not variable parameters.
2. tDepVars gives the time-dependent variables. In this case, it is just {a[t],b[t]}.
3. {t,t0,tf} (hopefully self-explanatory)
4. variableParameters, which is simply the list of the parameters that are being allowed to float.

After calling ParametricNDSolve, I need to manipulate the resulting solution into a form that is useful outside of this function. Specifically, I do (#->(#/.solution))&/@tDepVars (which by itself just gives me back what I started with, but which I will modify by working on the #/.solution part). I then do @@variableParameters so that I can plug in the values of the parameters using replacement rules rather than as function calls. But I need to modify the variable parameters so that the initial condition parameters are in the original form (e.g., a[0]) rather than the internal variable form (e.g., ic2), so my /.(a_->b_)->(b->a) bit does that.

Notice the commented-out Print statement. This is the place I have an issue.

But first we need some test calls to work with. Here we go:

DynamicModule[{result, eqs, tDeps, tfinal, rateConst, rateConstValues,
               parFixed, initialConc, ICFixed}, 
  eqs = {a'[t] == -k1 a[t] - k2 a[t]^2, b'[t] == k1 a[t] + k2 a[t]^2}; 
  initialConc = {1, 0}; 
  tDeps = {a, b}; 
  tfinal = 1;
  rateConst = {k1, k2}; 
  rateConstValues = {5, 5}; 
  parFixed = {True, True}; 
  ICFixed = {True, True}; 
  result = Kineticstest[eqs, {rateConst, rateConstValues}\[Transpose], 
           {tDeps, initialConc}\[Transpose], {t, 0, tfinal}, 
           parFixed, ICFixed]; 
  Deploy@Row[{Panel[
    Column[{Row[{"Fix k1? ", Checkbox[Dynamic[parFixed[[1]]]]}], 
            Row[{"Fix k2? ", Checkbox[Dynamic[parFixed[[2]]]]}], 
            Row[{"Fix a[0]? ", Checkbox[Dynamic[ICFixed[[1]]]]}], 
            Row[{"Fix b[0]? ", Checkbox[Dynamic[ICFixed[[2]]]]}], 
            Button["Recalculate", result = 
              Kineticstest[eqs, {rateConst, rateConstValues}\[Transpose], 
              {tDeps, initialConc}\[Transpose], {t, 0, tfinal}, 
              parFixed, ICFixed]]}]], 
    Panel[Dynamic[Plot[
      Evaluate[(#[t] & /@ tDeps) /. result /. ((#1 -> #2) & @@@ 
       Join[{rateConst, rateConstValues}\[Transpose], {#[0] & /@ tDeps, 
          initialConc}\[Transpose]])], {t, 0, tfinal}]]]}]]

The way to use this is to execute both blocks of code, uncheck the k1 box, and click Recalculate. Then uncomment the print statement, reexecute, and follow the same procedure.

Update: Print[solution//FullForm]; does not have the same effect (meaning that with FullForm, the solution doesn't display). Which gave me an idea... print it out FullForm before and after a normal print, as follows:

Print[solution//FullForm];
Print[solution];
Print[solution//FullForm];

That should tell me what has changed, right? Except the two FullForm versions that get printed are character-for-character identical. This is really confusing!

Update 2: Exactly the same behavior is seen if the Print is inside the Button.

Update 3: So if you uncheck more than one box, it works fine. It is only when there is a single variable parameter that there is a problem. This is making me wonder if somehow, when there is a single variable parameter, MMA is getting confused about whether the bracketed parameter is the parameter or is t. This would mean that the issue is with this line:

(# -> 
  ((# /. solution) @@
    (variableParameters /. 
       (internalSubstitutions /. (a_ -> b_) -> (b -> a))
     )
   )
 ) & /@ tDepVars

Maybe the @@ bit correctly identifies what follows as the parameter list when there are zero or two-or-more, but the @@ bit plugs in for t when there is only one parameter?

Update 4: Major progress in tracking down exactly what line is the problem. I was building up from a direct call of ParametricNDSolve with hard-coded parameters, adding in levels of sophistication, and checking at each level that the code still worked. I had no problems until I added in this part of the Plot statement, replacing a hard-coded {k1->5}:

((#1 -> #2) & @@@ Join[{rateConst, rateConstValues}\[Transpose], 
                       {#[0] & /@ tDeps, initialConc}\[Transpose]])

And lo, it fails without a Print and succeeds when there is one.

I am now working to remove all of the other bits I added to get a truly minimal example.

Answer 1

Your issue is that Dynamic doesn't know it has to refresh the value of tst2 inside of the Dynamic object. When you print the value of tst, then Mathematica realizes that there is a dependency in the plot on the printed value of tst (which contains the symbol k), and so it updates the plot dynamically.

One idea is to include k as a DynamicModule variable, and change the definition of tstfn so that it explicitly depends on k:

tstfn[v_] := Module[{tst,r},
    tst = ParametricNDSolve[
        {a'[t] == -k a[t], b'[t] == k a[t], a[0] == 1, b[0] == 0},
        {a[t], b[t]},
        {t, 0, 1},
        {k}
    ];
    (#->((#/.tst)@@{v}))&/@{a[t],b[t]}
]

DynamicModule[{k=5},
    Column[{
        Slider[Dynamic[k], {1, 10}],
        Dynamic[Plot[Evaluate[({a[t],b[t]} /. tstfn[k])], {t, 0, 1}]]
    }]
]

Also, I would avoid the contortions you use to modify the output of ParametricNDSolve, and just change to using ParametricNDSolveValue instead, something like:

DynamicModule[{pf, k},
    pf = ParametricNDSolveValue[
        {a'[t] == -k a[t], b'[t] == k a[t], a[0] == 1, b[0] == 0},
        {a[t], b[t]},
        {t, 0, 1},
        {k}
    ];
    Column[{
        Slider[Dynamic[k], {1, 10}],
        Dynamic[Plot[Evaluate[pf[k]], {t, 0, 1}]]
    }]
]