# Plot 1D slice of 2D InterpolatingFunction

I'm numerically solving a PDE (time+1D space) with NDSolve and want to plot some spatial profiles. Plot works, but is very slow compared to ListPlot. Here's an example (Fisher-KPP equation):

l = 100;
tmax = 30;
pts = 1000;

sol = NDSolve[{
D[n[x, t], t] == n[x, t] (1 - n[x, t]) + D[n[x, t], {x, 2}],
n[x, 0] == If[45 < x < 55, 1, 0], n[0, t] == n[l, t]}, {n}, {t, 0, tmax}, {x, 0, l},
Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> pts}}
][];

First@Timing[Plot[Evaluate[n[x, 15] /. sol], {x, 0, l}]]
(* 13.0701 *)

First@Timing[ListLinePlot[Table[{x, n[x, 15] /. sol}, {x, 0, l}]]]
(* 0.064737 *)


Can Plot be sped up? Maybe there is a nice way to extract a 1D InterpolatingFunction slice from a 2D InterpolatingFunction?

• If you don't need too much accuracy, FunctionInterpolation[] is usable for slicing. – J. M. is away Oct 14 '16 at 21:19

## 2 Answers

Extract the appropriate values and re-interpolate. If we extract the derivative values (optional), you get the normal cubic Hermite interpolation type produced by NDSolve on a first-order ODE. If t is a grid point, then the derivatives at the values of x on the grid should be stored (as computed from the PDE) in the InterpolatingFunction. Between grid points, these will be interpolated. It makes little difference to the plot, so you can omit and reduce the time a little more.

(if = n /. First@sol;
grid = if["Grid"];
slice = Interpolation@Transpose@{
grid[[All, 1, {1}]],                (* extract x-grid *)
if[grid[[All, 1, 1]], 15.],         (* extract n-values on x-grid at t == 15 *)
Derivative[1, 0][if][grid[[All, 1, 1]], 15.] (* derivative values *)
}) // AbsoluteTiming Plot[slice[x], {x, 0, 100}] // AbsoluteTiming • Thanks, this is great! I used this as the basis for my own answer below. – Chris K Oct 15 '16 at 15:11

@MichaelE2's slice-taking solution works great, but it's a bit cumbersome to use. Here's an idea based on it, which overloads the definition of InterpolatingFunction that makes it transparent to use, and also work on the other dimension.

Unprotect[InterpolatingFunction];

InterpolatingFunction[stuff___][var_Symbol, num_?NumericQ] := Module[{if, grid},
if = InterpolatingFunction[stuff];
grid = (InterpolatingFunction[stuff])["Grid"];

Return[(Interpolation@Transpose@{
grid[[All, 1, {1}]], (* extract x-grid *)
if[grid[[All, 1, 1]], num]
})[var]]
];

InterpolatingFunction[stuff___][num_?NumericQ, var_Symbol] := Module[{if, grid},
if = InterpolatingFunction[stuff];
grid = (InterpolatingFunction[stuff])["Grid"];

Return[(Interpolation@Transpose@{
grid[[1, All, {2}]], (* extract t-grid *)
if[num, grid[[1, All, 2]]]
})[var]]
];

Protect[InterpolatingFunction];


In action:

l = 100;
tmax = 30;
pts = 1000; (* # grid points *)

sol = NDSolve[{
D[n[x, t], t] == n[x, t] (1 - n[x, t]) + D[n[x, t], {x, 2}],
n[x, 0] == If[45 < x < 55, 1, 0], n[0, t] == n[l, t]}, {n}, {t, 0, tmax}, {x, 0, l},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> pts}}][];

First@Timing[Plot[Evaluate[n[x, 15] /. sol], {x, 0, l}]]
(* 0.3865 *)

Plot[Evaluate[n[x, 15] /. sol], {x, 0, l}] Plot[Evaluate[n[20, t] /. sol], {t, 0, tmax}] This is my first attempt at modifying built-in functions, so let me know if it might be improved or if it might break anything!

• I've since been warned that modifying built-in commands is generally bad practice. This solution worked nicely for me, but use it at your own risk! – Chris K Oct 31 '16 at 23:43