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

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;

]
• Perhaps you could make use of an NDSolve state object? Feb 12 at 6:06

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}

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 *)
(* {{}} *)
• Chris, thanks for the suggestion. It may work for the simple ODE that you use as example, but it does not work for my PDE problem. On attempt to use empty list, there appears a warning message: NDSolve: There are more dependent variables {list} than equations, so the system is underdetermined. I can't even paste the message, because Mathematica goes into freeze and crashes.
– Mike
Feb 11 at 18:22
• I see. Hopefully someone with more info will chime in, but in the meanwhile you might try to make a more self-contained, minimal example to get more people interested in it. Feb 11 at 18:40
• BTW, my question can be stated differently, and the solution may be also useful and not only to me. My understanding is that the reason behind slow creation of interpolation function upon NDSolve completion in many problems is that it contains the entire time history, from tinitial to tend. Is there a way to tell NDSolve to create it at needed time instances only ? In my case, this would be at tend only.
– Mike
Feb 11 at 18:42
• You could try {t, tend, tend} but I think @tad's answer is more elegant. Feb 11 at 20:10
• Perhaps you could get it work by specifying "DependentVariables" Feb 12 at 6:07

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.

• Thanks ! This works just fine.
– Mike
Feb 17 at 23:11