# How to integrate the interpolated data over the region in which it is non-zero?

Consider some toy-function which I then interpolate:

Function1[x_, y_] = If[4<x^3 + y^4 - x^2 < 5, Exp[-x^2 + y^2], 0];
DataInt =
Flatten[ParallelTable[{x, y, Function1[x, y]}, {x, 0, 100,
0.05}, {y, 0, 100, 0.05}], {2, 1}];
FunctionInt[x_, y_] =
Interpolation[DataInt, InterpolationOrder -> 1][x, y];


Note that there is a lot of points in which the interpolated function is zero. I want then to integrate this function:

NIntegrate[FunctionInt[x, y], {x, 0, 100}, {y, 0, 100}]


The integral gives zero due to the fact that the domain in which FunctionInt is non-zero is very narrow. An brute-force way to improve this is to increase the integration accuracy. However, it may be more useful to define somehow the domain in which the function is non-zero, i.e. to get a function y = y[x] which defines this region.

Is it possible to do having the data?

• Here's a quickish Riemann sum estimate: SparseArray[Head[FunctionInt[x, y]]["ValuesOnGrid"]]["NonzeroValues"] // Total[#] 0.05^2 & -- tends to overestimate the trapezoidal rule on a nonnegative function. Jan 7, 2021 at 18:49
• Do you know the symbolic description of the region such as 4<x^3 + y^4 - x^2 < 5 for your actual use case? Do you know the formula for the data being interpolated such as Exp[-x^2 + y^2]? Jan 7, 2021 at 21:32
• @MichaelE2 : I don't know the symbolic description and the formula. Jan 11, 2021 at 9:15

You can use Integrate to accumulate an InterpolatingFunction:

Integrate[FunctionInt[x, y], x, y] /. {x -> 100, y -> 100}
(*  0.601376  *)


If the interpolation grid is regular as in the example (spacing = 0.05), here's a quick way using a manual trapezoidal rule:

With[{fvals = Head[FunctionInt[x, y]]["ValuesOnGrid"]},
Nest[Total[Most[#] + Rest[#]] 0.05/2 &, fvals, 2]
]
(*  0.601376  *)


A more general trapezoidal rule is implemented here: Choosing points of integration in NIntegrate (integrating a function given by a list without interpolation)

• The general trapezoidal rule refactored for this problem: ifn = Head[FunctionInt[x, y]]; With[{fvals = ifn@"ValuesOnGrid"}, Fold[Function[{y, dx}, dx . (Most[y] + Rest[y])/2], fvals, Differences /@ ifn@"Coordinates"] ] Jan 7, 2021 at 19:21

Find the domain over which the values are not zero:

minmaxes = MinMax /@ Transpose@Select[DataInt, #[] != 0 &]

(* Out: {{0., 2.1}, {0., 1.5}, {0.0121552, 8.67114}} *)


Select a square region that includes that domain and interpolate over that:

squareint =
Interpolation[
Select[DataInt, 0 <= #[] <= 2.1 && 0 <= #[] <= 1.5 &]
];


Integrate over that interpolated region numerically:

NIntegrate[
squareint[x, y], {x, y} ∈ Rectangle @@ Transpose@squareint["Domain"]
]

(* Out:  0.554679 *)


Notice that you will get a NIntegrate::slwcon warning probably because of the perceived discontinuity over this domain.

This result compares well with the true value:

NIntegrate[
Boole[4 < x^3 + y^4 - x^2 < 5] Exp[-x^2 + y^2],
{x, 0, 100}, {y, 0, 100}
]

(* Out: 0.597656 *)

• Keep in mind that the OP's interpolating function will perform a linear extrapolation (InterpolationOrder -> 1) beyond the boundary. Jan 7, 2021 at 18:42

With the data:

Function1[x_, y_] = If[4 < x^3 + y^4 - x^2 < 5, Exp[-x^2 + y^2], 0];
DataInt =
Flatten[ParallelTable[{x, y, Function1[x, y]},
{x, 0, 100, 0.05}, {y, 0, 100, 0.05}], {2, 1}];
FunctionInt[x_, y_] = Interpolation[DataInt, InterpolationOrder -> 1][x, y];


You could try:

ir = ImplicitRegion[{4 < x^3 + y^4 - x^2 < 5, 0 < x < 100, 0 < y < 100}, {x, y}];
Region@ir NIntegrate[FunctionInt[x, y], {x, y} ∈ ir]
(*0.468457*)


However I get the warning message:

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

Something like this is not too surprising for a numerical integration routine and a discontinuous function.

\$Version

(* "12.2.0 for Mac OS X x86 (64-bit) (December 12, 2020)" *)

Clear["Global*"]


Use Piecewise rather than If for numeric functions

Function1[x_, y_] = Piecewise[{{Exp[-x^2 + y^2], 4 < x^3 + y^4 - x^2 < 5}}];

int1 = NIntegrate[Function1[x, y], {x, 0, 100}, {y, 0, 100}]

(* 0.597656 *)

reg = ImplicitRegion[
4 < x^3 + y^4 - x^2 < 5 && 0 <= x <= 100 && 0 <= y <= 100, {x, y}];

int2 = NIntegrate[Exp[-x^2 + y^2], {x, y} \[Element] reg]

(* 0.597656 *)


The results are identical

int1 === int2

(* True *)
`
• You don't seem to integrate an interpolating function? Or did I miss it? (The actual use-case may involved interpolation of data as indicated in the question, not a symbolic function.) Jan 7, 2021 at 21:29