StreamPlot with interpolating functions

Question

While working my way through the Wolfram Blog post on simulating fluid flow I started playing with StreamPlot because I preferred quick and ugly visualizations over the ListIntegralConvolutionPlot which simply tested my patience. When using StreamPlot with the interpolated velocity vectors, however, the function would throw an error:

InterpolatingFunction::dmval: Input value {1.05051,1.05051} lies outside the range of data in the interpolating function. Extrapolation will be used. >>

An example that reproduces this error is:

u = Interpolation[
  Flatten[Table[{{x, y}, -1 - x^2 + y}, {x, 0, 1, 0.1}, {y, 0, 1, 
     0.1}], 1]]
v = Interpolation[
  Flatten[Table[{{x, y}, 1 + x - y^2}, {x, 0, 1, 0.1}, {y, 0, 1, 
     0.1}], 1]]
StreamPlot[{u[x, y], v[x, y]}, {x, 0, 1}, {y, 0, 1}]

It looks like a similar issue was observed in earlier versions of Mathematica with LogLinearPlot. Following a similar troubleshooting process:

Reap@StreamPlot[{Sow[x], Sow[y]}, {x, 0, 1}, {y, 0, 1}];
%[[2]]

shows (I think, I'm not terribly comfortable with Reap/Sow) that StreamPlot is polling the function at 1.05051.

{{1.05051, 1.05051, 1.05051, 1.05051, 0.,...}

Furthermore, the polling seems independent of the range for x or y supplied in the function.

In the end, the plot that is displayed with interpolated values is identical to the one made with the original equations (in this case). My concern is that this behavior may impact a subset of interpolated functions that do behave badly outside of their defined region and it would be good to be able to tell StreamPlot that. Is it possible or is this a bug that we'll be stuck with until version 10 comes out?

Answer 1

What StreamPlot is doing with the values at points outside of the domain seems to have to do with finding a good arrangement of stream lines. The values at points outside the domain do not appear in the output, so your concern might turn out to be unwarranted. Bad behavior due to extrapolation might result only in an ugly but correct plot. I'm not sure why it has to look outside the domain to make a good plot.

A classic way of preventing extrapolation is to limit the definition of a function:

Clear[vel];
vel[x_, y_] /; 0 <= x <= 1 && 0 <= y <= 1 := {u[x, y], v[x, y]};
velsp = StreamPlot[vel[x, y], {x, 0, 1}, {y, 0, 1}, StreamStyle -> ColorData[1][3]];

This results in a stream plot with an uneven distribution of stream lines compare to the original. (See below.)

We can test what happens if the function has odd values outside the domain by explicitly extending the vector field discontinuously.

Clear[disc];
disc[x_, y_] := If[0 <= x <= 1 && 0 <= y <= 1, {u[x, y], v[x, y]}, {-100., -100.}];
discsp = StreamPlot[disc[x, y], {x, 0, 1}, {y, 0, 1}, StreamStyle -> ColorData[1][2]]

This produces a very similar plot, but it is different. The difference is perhaps worrying, but in this case we can see that the difference is negligible.

sp = StreamPlot[{u[x, y], v[x, y]}, {x, 0, 1}, {y, 0, 1}]; (* OP's *)
GraphicsRow[{sp, discsp, velsp}, ImageSize -> 600]

Comparing the discp with sp directly, the differences are almost imperceptible (especially in the image below). A few streamlines seem to be shifted over just barely.

Show[sp, discsp]

The similarity of the two plots supports the idea that the values of the vector field in the plot domain determine the streamlines. The rightmost plot velsp, while a bit ragged, still shows the flow of the vector field (which may be confirmed with Show[sp, velsp]). The ragged look suggests that one ought not to worry about StreamPlot going outside the interpolation domain.

Aside

I don't get the dmval message when I copy/paste/execute your example (Mma V9.0.1, MacOSX 10.8.4), but if I track the function evaluations (using Reap and the option EvaluationMonitor :> Sow[{x, y}]), values of x and y that are 0.02... outside of [0, 1] are being plugged into the function. I'm not sure why mine is Quiet and yours isn't, nor why yours are 0.05 outside the domain and mine 0.02. (A clue might be that the full value 0.0217563... corresponds exactly with the PlotRange of the plot.) I suspect that the Window's version is meant to be quiet, but for some, perhaps unintentional, reason isn't.