How can I prevent a Plot expression from evaluating points forbidden by its RegionFunction?

Question

I just ran into this and it had me scratching my head for a good while. Say I have a function which may produce error messages for certain inputs, like e.g.

f[x_, y_] := (Sin[x] - Sin[y])/(x - y)

which will return Power::infy and Infinity::indet when called on the diagonal x==y. If I want to avoid that singular input in a plot, the standard way of telling Mathematica to not evaluate it there is to use a RegionFunction command, like

Plot3D[f[x, y]
 , {x, -10, 10}, {y, -10, 10}
 , RegionFunction -> Function[{x, y, f}, Abs[x - y] > 5]
 ]

This will successfully generate a plot with the diagonal taken off, but it will still produce the same error messages as f[1,1]. This behaviour remains even if I add the qualifiers f[x_: NumericQ, y_: NumericQ] to the definition, so some singular inputs are definitely passed to the function even though they are away from where I'm telling it to plot.

So: why is this? Is this actually the correct way to specify regions to stay away from? If not, what's the best practice for that?

Answer 1

One option would be to restrict the function from funky regions with a Condition liks this

f[x_, y_] /; Abs[x - y] > 5 := (Sin[x] - Sin[y])/(x - y);
Plot3D[f[x, y], {x, -10, 10}, {y, -10, 10}]

Out:

One can easily see that Plot3D will also sample points in the region which is "forbidden", the points are just not drawn due to the RegionFunction.

ClearAll[f];
f[x_, y_] /; Abs[x - y] > 5 := (Sin[x] - Sin[y])/(x - y);
xyValues = {};
Plot3D[xyValues = Append[xyValues, {x, y}];
  f[x, y], {x, -10, 10}, {y, -10, 10}, Evaluated -> False, 
  RegionFunction -> Function[{x, y, f}, Abs[x - y] > 5]];
ListPlot[xyValues]

Out:

thus mathematica will print some errors once the points close to zero are evaluated, if the function is not restricted. No errors will occure, if the function is restricted, because in that case just the Input will be returned:

ClearAll[f];
f[x_, y_] /; Abs[x - y] > 5 := (Sin[x] - Sin[y])/(x - y);
f[0, 0]

Out:

f[0, 0]

Another possible solution would be to solve the 0/0 Problem analyticaly like this:

Limit[(Sin[x] - Sin[y])/(x - y), x -> y]

Out:

Cos[y]

and thus redefine our function to account for the x==y situation.

ClearAll[f];
f[x_, y_] := (Sin[x] - Sin[y])/(x - y)
f[x_, y_] /; x == y := Cos[y]
Plot3D[f[x, y], {x, -10, 10}, {y, -10, 10}]

which gives the function over the whole area.

edit: added some pictures and a little more Info.