35

One of the most annoying "features" of Mathematica is that the Plot family does extrapolation on InterpolatingFunctions without any warning. I'm sure it was discussed to hell previously, but I cannot seem to find any reference. While I know how to simply overcome the problem by defining a global variable for the domain of the interpolation, from time to time I forget to do this and then I spend days figuring out where the numerical error originates. This could be avoided if Plot was to give a warning.

Consider the following example. An ODE system is defined and integrated for two different time ranges:

odes = {
   a'[t] == -a[t] - .2 a[t]^2 + 2. b[t], 
   b'[t] == a[t] + .1 a[t]^2 - 1.1 b[t], a[0] == 1, b[0] == 1
};
sol100 = First@NDSolve[odes, {a, b}, {t, 0, 100}];
sol500 = First@NDSolve[odes, {a, b}, {t, 0, 500}];

Now querying the function value for a point outside of the range correctly gives a warning:

(a /. sol100)[500]
InterpolatingFunction::dmval: Input value {500} lies outside
the range of data in the interpolating function. Extrapolation will be used. >>

651.034

The same is not done when we use the function in Plot:

Show[
 Plot[{a[t], b[t]} /. sol100, {t, 0, 400}, PlotStyle -> {Thick, Red}],
 Plot[{a[t], b[t]} /. sol500, {t, 0, 400}, PlotStyle -> {Thick, Blue}]
 ]

Mathematica graphics

I've tried to force a warning, with no avail. The following example won't give a warning.

On[InterpolatingFunction::dmval]
Check[Plot[{a[t], b[t]} /. sol100, {t, 0, 500}], "Error", 
 InterpolatingFunction::dmval]

Interestingly, one can be sure that InterpolatingFunction::dmval is NOT turned off at all inside the Plot family. In the following example, LogLinearPlot is able to drop a warning about sampling from below the domain (that can be ignored being unrelated, see this post, also it seems to be fixed in v9), but it does not give the same warning when sampling from above (> 100)!

LogLinearPlot[{a[t], b[t]} /. sol100, {t, 0.1, 500}]
InterpolatingFunction::dmval: Input value {-2.30241} lies outside the
range of data in the interpolating function. Extrapolation will be used. >>

It is even more disturbing to see that Plot checks the lower boundary but not the upper (thanks to J.M. for the comment):

Plot[{a[t], b[t]} /. sol100, {t, -1, 500}]
InterpolatingFunction::dmval: Input value {-0.989765} lies outside the
range of data in the interpolating function. Extrapolation will be used. >>

As Oleksandr has pointed out, it is not about lower vs. upper boundaries but first point vs. the rest.

Plot[{a[t], b[t]} /. sol100, {t, 101, 500}]
InterpolatingFunction::dmval: Input value {101.008} lies outside the
range of data in the interpolating function. Extrapolation will be used. >>

Questions

  1. Why Plot does not give a warning when extrapolating an InterpolatingFunction? Is there some higher-level consideration that justifies this behaviour, or is it a bug?
  2. How can one force Plot to give a warning? Is there any workaround that forces InterpolatingFunction::dmval not to be attenuated inside Plot?
István Zachar
  • 47,032
  • 20
  • 143
  • 291
  • as a remark Plot also ignores say complex values which is sometimes useful. – chris Nov 05 '12 at 12:20
  • Plot[Check[{a[t], b[t]} /. sol100, "Error"], {t, 0, 500}] ? – Dr. belisarius Nov 05 '12 at 12:21
  • @belisarius: It stops plotting because of Plot::accbend and not because of extrapolating outside of the range (InterpolatingFunction::dmval). It is not informative for the user, and I'm not sure they are related. – István Zachar Nov 05 '12 at 13:14
  • 1
    Very odd. OTOH, try Plot[{a[t], b[t]} /. sol100, {t, -1, 500}]. – J. M.'s missing motivation Nov 05 '12 at 13:35
  • Thanks @J.M., included, if you don't mind. – István Zachar Nov 05 '12 at 13:43
  • You're welcome; helping with the diagnosis is why I wrote that comment. :) If memory serves, older versions of Mathematica allowed InterpolatingFunction[] to spit out a warning when extrapolation is being done within Plot[]. I wonder what changed... – J. M.'s missing motivation Nov 05 '12 at 13:50
  • 1
    Strange indeed, and well spotted. Compare Plot[{a[t], b[t]} /. sol100, {t, 400, 500}]. It seems to be determined not by lower vs. upper bounds but rather somehow depends on the range chosen (first point vs. the rest?); apparently InterpolatingFunction::dmval is turned off, but only later in the plotting process. The strange thing is that hooking Quiet and Off show no relevant usage of the former in Plot and no instances at all of the latter. – Oleksandr R. Nov 05 '12 at 13:56
  • @J.M. I do remember encountering this odd behaviour as early as 2007 for the first time. Before that, I never really used NDSolve and interpolated functions. – István Zachar Nov 05 '12 at 13:57
  • Furthermore, hooking Message shows that rather a lot of InterpolatingFunction::dmval messages are produced--they just somehow don't get out of Plot. Weird, and in my estimation, almost certainly a bug... – Oleksandr R. Nov 05 '12 at 14:12
  • Tagging this as bug. @Oleksandr Perhaps there's a non-selective Quiet somewhere inside. – Szabolcs Nov 13 '12 at 20:02
  • @Szabolcs definitely a bug. As Mr. W's answer shows, Quiet isn't used; nor is Off. (I had hooked these two, plus Message, before I left my comment above. There is something going on with the internal message handling routines but I am not sure how these work and didn't have time to figure it out.) – Oleksandr R. Nov 14 '12 at 13:03

3 Answers3

8

It seems that indeed messages are suppressed somewhere in the evaluation of Plot. Interestingly this appears to happen some time after the first numeric evaluation. If we issue this message for each (numeric) point plotted we see that it is printed only once:

f[x_?NumericQ] := (Message[InterpolatingFunction::"dmval", {21.06`}]; Sin[x])

Plot[f[x], {x, 0, 10}]

If we make f to redefine itself such that it only issues the message after the first point is evaluated the message is never printed:

f[x_?NumericQ] :=
 (
  f[y_?NumericQ] := (Message[InterpolatingFunction::"dmval", {21.06`}]; Sin[y]);
  Sin[x]
 )

Plot[f[x], {x, 0, 10}]

This shows that the suppression of this message is not the result of some special optimization while plotting an InterpolatingFunction (which might be done outside the normal evaluation chain) but rather something more general.

It should be noted that the ::dmval message is not the only one suppressed, e.g.:

f[x_?NumericQ] := (Message[qrs::argx, 1, 2]; Sin[x])

Plot[f[x], {x, 0, 10}]

We can further demonstrate that Plot does not (appear to) use the high level function Quiet itself for this suppression, as this does not restore the messages:

Block[{Quiet = # &}, Plot[f[x], {x, 0, 10}] ]

One can get messages past Plot with something like this:

Unprotect[Message]

Message[x : InterpolatingFunction::"dmval", other___] :=
  Print @ Style[Row[{HoldForm[InterpolatingFunction::"dmval"], ":  ", 
      StringForm[x, other]}], "Message"]

f[x_?NumericQ] := (Message[InterpolatingFunction::"dmval", {21.06`}]; Sin[x])

Plot[f[x], {x, 0, 10}, PlotPoints -> 5, MaxRecursion -> 0]

But there is no automatic limit to the number of messages printed making this impractical to use as written.

I am continuing to explore this question.

Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371
5

As you can see with

if = Interpolation[Range[10]];
With[{mess = $Messages},
      Plot[
       Block[{$Messages = mess}, if[x]], {x, 5, 20}, Evaluated -> False]
  ]~Quiet~Message::msgl

the problem is that $Messages gets blocked somehow. I have so far no good recommendations as to how to fix it properly.

Rojo
  • 42,601
  • 7
  • 96
  • 188
4

As Rojo showed, the symbol $Messages appears to become unset somewhere inside the Plot internals. I wondered if it was possible to prevent this by setting the Protected and Locked attributes for $Messages, but it would seem that Plot has magic powers:

SetAttributes[$Messages, {Protected, Locked}];
f = Interpolation[Range[5]];
g[x_] := (Print[Attributes[$Messages], "  ", $Messages]; f[x]);
Plot[g[x], {x, 1, 6}]

enter image description here

You can see that after the first evaluation $Messages loses the Protected attribute and its value, though curiously it remains Locked throughout.

Since it seems impossible to stop $Messages losing its value, an imperfect workaround might be to reset it whenever the InterpolatingFunction::dmval message is generated.

Quit[];

Unprotect[Message];
With[{mess = $Messages},
 m : Message[InterpolatingFunction::dmval, __] :=
  Block[{$mIFdmval = True}, $Messages = mess; m] /; !TrueQ[$mIFdmval]]

With this approach the automatic limit of 3 messages kicks in if the messages are generated outside of Plot:

f = Interpolation[Range[5]];
Table[f[x], {x, 1, 10}]

enter image description here

Unfortunately, the automatic limit doesn't work with Plot (perhaps because Plot messes around with $MessageList as well), so you still get a long list of messages.

Simon Woods
  • 84,945
  • 8
  • 175
  • 324