How to tell NDSolve to suppress creation of interpolation function on successful exit?

Question

Update: I tried to work around DependentVariables, as suggested (thanks !) Here is the example from documentation for an ODE system, where InterpolatingFunction is created only for the variable x: xsol = NDSolve[{x'[t] == y[t], y'[t] == -Sin[x[t]], x[0] == 3.1, y[0] == 0}, x, {t, 0, 20}, DependentVariables -> {x, y}]

And here is a simple PDE, where same idea does not work (note commented out u):

 NDSolve[{D[u[x, y, t], t] == 
 D[u[x, y, t], x, x] + D[u[x, y, t], y, y], 
 u[0, y, t] == u[1, y, t], u[x, 0, t] == u[x, 1, t], 
 u[x, y, 0] == 
 10/Sqrt[2 \[Pi]] Exp[-50 ((x - 1/2)^2 + (y - 1/2)^2)]}, 
 (*u,*){x, 0,
 1}, {y, 0, 1}, {t, 0, 1}, DependentVariables -> u[x, y, t], 
 AccuracyGoal -> MachinePrecision/3, 
 PrecisionGoal -> 1 + MachinePrecision/3, InterpolationOrder 
 -> All, 
 MaxSteps -> 10000, Compiled -> False, 
 Method -> {"MethodOfLines", 
 "SpatialDiscretization" -> {"TensorProductGrid", 
 "DifferenceOrder" -> 4, MinPoints -> 110, MaxPoints -> 110, 
 AccuracyGoal -> MachinePrecision/3, 
 PrecisionGoal -> 1 + MachinePrecision/3}, Method -> "BDF"}]

and even if I write: DependentVariables -> u it still does not work. Again, the goal is to get rid of InterpolatingFunction as the answer to the problem, or create InterpolationFunction only at t=1, i.e. at the endpoint of integration.

---

This may seem to many as stupid question to ask, but I would like to know if creation of interpolation function in NDSolve can be eliminated. The reason is, I call NDSolve in a loop, and when a particular tend is reached I dump the solution on a grid to a file, then read the file, interpolate data, and use this data as the initial condition for next call to NDSolve. I don't need the interpolation function creation at the end of each step in a loop, and I want to get rid of it, because it takes around one hour just to create this interpolating function.

My question is similar to the discussion here, but it is different.

Here is my code, except I do not state the right-hand side of a PDE, since it is too long and complicated.

tini = 0; \[CapitalDelta]t = 1;

For[i = 0, i <= 100, i++,

AbsoluteTiming[
NDSolve[{D[Subscript[C, B][x, y, t], t] == RHSCeqB, 
 Subscript[C, B][0, y, t] == Subscript[C, B][domainlengthX, y, t],
  Subscript[C, B][x, 0, t] == 
  Subscript[C, B][x, domainlengthY, t], 
 Subscript[C, B][x, y, tini] == ControlShapeCB[x, y], 
 WhenEvent[
  Mod[t, tini + \[CapitalDelta]t] == 
   0, {Print["t=", t, "   ", "i=", i], 
   Export["c:\\file-out-CB" <> ToString[i] <> ".csv", 
    Flatten[
     Table[{N[x], N[y], Subscript[C, B][x, y, t]}, {x, 0, 
       domainlengthX, 
       domainlengthX/200}, {y, 0, domainlengthY, 
       domainlengthY/200}], 1]]}]}, Subscript[C, 
   B], {x, 0, domainlengthX}, {y, 0, domainlengthY}, {t, tini, 
   tini + \[CapitalDelta]t}, AccuracyGoal -> MachinePrecision/3, 
   PrecisionGoal -> 1 + MachinePrecision/3, 
   InterpolationOrder -> All, MaxSteps -> 10000, Compiled -> False, 
   Method -> {"MethodOfLines", 
   "SpatialDiscretization" -> {"TensorProductGrid", 
    "DifferenceOrder" -> 4, MinPoints -> 200, MaxPoints -> 200, 
    AccuracyGoal -> MachinePrecision/3, 
    PrecisionGoal -> 1 + MachinePrecision/3}, 
   Method -> {"ImplicitRungeKutta", 
    "Coefficients" -> "ImplicitRungeKuttaGaussCoefficients", 
    "DifferenceOrder" -> 1}}]] // Echo;

   dataC = Import["c:\\file-out-CB" <> ToString[i] <> ".csv"]; 
   ControlShapeCB = Evaluate@Interpolation[dataC, Method -> "Spline"];

   tini = tini + \[CapitalDelta]t;

   ]

Answer 1

You can tell NDSolve to only return the final value instead of an interpolation over the time interval, as noted in the NDSolve documentation (2nd line of Details section).

For example:

NDSolve[{n'[t] == n[t] (1 - n[t]), n[0] == 0.1}, n[1], {t, 0, 1}]

which gives

{{n[1] -> 0.231969}}

Or more directly using NDSolveValue:

NDSolveValue[{n'[t] == n[t] (1 - n[t]), n[0] == 0.1}, n[1], {t, 0, 1}]

Use a list to get values at specified times, e.g.,

NDSolveValue[{n'[t] == n[t] (1 - n[t]), n[0] == 0.1}, {n[0.5], n[1]}, {t, 0, 1}]

to get

{0.154828, 0.231969}

Answer 2

Just replace the list of variables being solved for with an empty list. E.g.

NDSolve[{n'[t] == n[t] (1 - n[t]), n[0] == 0.1, 
  WhenEvent[n[t] == 0.5, Print[t]]}, n, {t, 0, 10}]
(* 2.19722 *)

NDSolve[{n'[t] == n[t] (1 - n[t]), n[0] == 0.1, 
  WhenEvent[n[t] == 0.5, Print[t]]}, {}, {t, 0, 10}]
(* 2.19722 *)
(* {{}} *)

Answer 3

I think this is what you want. In the docs somewhere it says NDSolve with {t, a, b} and an initial condition at t == c with c < a won't save the steps until t == a. If a == b, then only the last step is saved. This feature exists specifically to save memory and the time it takes to process the data generated.

Quit[] (* restart kernel *)

foo = MaxMemoryUsed[]

(* 90504064 *)

sol = NDSolve[{D[u[x, y, t], t] == 
    D[u[x, y, t], x, x] + D[u[x, y, t], y, y], 
   u[0, y, t] == u[1, y, t], u[x, 0, t] == u[x, 1, t], 
   u[x, y, 0] == 
    10/Sqrt[2  \[Pi]]  Exp[-50  ((x - 1/2)^2 + (y - 1/2)^2)]}
  , u, {x, 0, 1}, {y, 0, 1}, {t, 10, 10},
  AccuracyGoal -> MachinePrecision/3, 
  PrecisionGoal -> 1 + MachinePrecision/3,(*InterpolationOrder->All,*)
  MaxSteps -> 10000, Compiled -> False, Method -> {"MethodOfLines",
    "SpatialDiscretization" -> {"TensorProductGrid"
      , "DifferenceOrder" -> 4, MinPoints -> 110, MaxPoints -> 110, 
      AccuracyGoal -> MachinePrecision/3, 
      PrecisionGoal -> 1 + MachinePrecision/3}
    , Method -> "BDF"}]

Note the domain: Only t == 10. is returned. The following shows the difference in memory usage. The first is from the above. The second is when {t, 0, 10} specifies the domain in NDSolve instead of {t, 10, 10}:

u /. sol // ByteCount
MaxMemoryUsed[] - foo

(* 479016 *)
(* 119816624 *)

u /. sol // ByteCount
MaxMemoryUsed[] - foo

(* 57508792 *)
(* 141724816 *)

Also note that Interpolation -> All when the end points are the same in the time integration, that is {t, 10, 10} in this case, seems to cause the kernel to crash. A bug probably. OTOH, with only one interpolation node, Interpolation -> All doesn't make sense anyway.