NIntegrate fails to give non-numerical value error, instead returns zero

Question

When I evaluate

NIntegrate[a(x - y)^2, {x, 0, 1}, {y, 0, 1}, AccuracyGoal -> 8]

Mathematica returns 0 without giving a non-numerical value error, despite that fact that there is an undefined a in the integrand. If I get rid of the AccuracyGoal specification then the expected error is raised. It is also raised if I multiply a by pretty much any function other than (x-y)^2.

I'm wondering why, in the particular instance describe above, NIntegrate does not warn me that I'm using incorrect syntax.

EDIT: Using evaluation monitor

Reap[NIntegrate[a (x - y)^2, {x, -1, 1}, {y, -1, 1}, AccuracyGoal -> 8, EvaluationMonitor :> Sow[{x, y}]]]

results in

{0., {{{0., 0.}, {0.948683, 0.948683}, {-0.948683, -0.948683}, {-0.688247, -0.688247}, {0.688247, 0.688247}}}}.

This explains why no error was given: NIntegrate evaluates the integrand only at points where it evaluates to the nonsymbolic value zero. I suppose my remaining question is why NIntegrate only evaluated the integrand at those points? If I do the same with an explicit value for a (say 1) I find

Reap[NIntegrate[1 (x - y)^2, {x, -1, 1}, {y, -1, 1}, AccuracyGoal -> 8, EvaluationMonitor :> Sow[{x, y}]]]

{2.66667, {{{0., 0.}, {0.358569, 0.}, {-0.358569, 0.}, {0.948683, 0.}, {-0.948683, 0.}, {0., 0.358569}, {0., -0.358569}, {0., 0.948683}, {0., -0.948683}, {0.948683, 0.948683}, {0.948683, -0.948683}, {-0.948683, -0.948683},{-0.948683, 0.948683}, {-0.688247, -0.688247}, {-0.688247,0.688247}, {0.688247, -0.688247}, {0.688247, 0.688247}}}}

in disagreement with the case where a is symbolic.

Answer 1

NIntegrate treats nonnumeric results as singularities and they are discarded. If there are numeric function values, it assumes that the singularities will be handled in the course of refinement. If there are no numeric function values at all sampling points, you get a NIntegrate::inumr error.

Try the following to see the integrand is actually evaluated throughout the domain.

ii[x_?NumericQ, y_?NumericQ] := (Sow[{x, y}, "All"]; a (x - y)^2);
Reap[
   NIntegrate[
    ii[x, y], {x, 0, 1}, {y, 0, 1},
    AccuracyGoal -> 8,
    EvaluationMonitor :> Sow[{x, y}, "Used"]],
   {"All", "Used"}
   ][[2, All, 1]] // 
 ListPlot[#, AspectRatio -> 1, PlotLegends -> {"All", "Used"}] &

Since only points with x == y are used, the integral and error are both 0.:

NIntegrate[a (x - y)^2, {x, 0, 1}, {y, 0, 1}, AccuracyGoal -> 8, 
  IntegrationMonitor :> (Sow[Through[#@"Error"]] &)] // Reap

(*  {0., {{{0.}}}}  *)

The OP observes that an error message does not appear until AccuracyGoal is at least 324. The actual limit is connected to the smallest (subnormal) machine-precision number, which is the same as Accuracy[0.], accuracy of the error:

ag = -Log10[$MinMachineNumber] + Log10[2^53]
Accuracy[0.]

(*
  323.607
  323.607
*)

NIntegrate[a (x - y)^2, {x, 0, 1}, {y, 0, 1}, 
 AccuracyGoal -> (1 - $MachineEpsilon/2) ag]  (* no error *)

(*  0.  *)

NIntegrate[a (x - y)^2, {x, 0, 1}, {y, 0, 1}, 
 AccuracyGoal -> ag]                          (*  error   *)

NIntegrate::inumr: The integrand a (x-y)^2 has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,0.5},{0,1}}. >>

Ignoring PrecisionGoal for this case, if you set AccuracyGoal higher than the Accuracy of the error, NIntegrate will subdivide the region with the largest error (which at the first step is the whole domain). As subdivision goes on, not all of the new regions will be symmetric with respect to the integrand and two things can happen: The integrand is nonnumeric on the entire subregion, or the error is at best 0. on a subregion. If the first happens, you get a NIntegrate::inumr error and integration stops. In the second, you get more subdivision.