5
$\begingroup$

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}}
][[1]];

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?

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

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

Mathematica graphics

Plot[slice[x], {x, 0, 100}] // AbsoluteTiming

Mathematica graphics

$\endgroup$
  • $\begingroup$ Thanks, this is great! I used this as the basis for my own answer below. $\endgroup$ – Chris K Oct 15 '16 at 15:11
3
$\begingroup$

@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}}][[1]];

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

Plot[Evaluate[n[x, 15] /. sol], {x, 0, l}]

Mathematica graphics

Plot[Evaluate[n[20, t] /. sol], {t, 0, tmax}]

Mathematica graphics

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!

$\endgroup$
  • 1
    $\begingroup$ 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! $\endgroup$ – Chris K Oct 31 '16 at 23:43

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.