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

Question

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. — Michael E2, Jan 07 '21 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]? — Michael E2, Jan 07 '21 at 21:32
@MichaelE2 : I don't know the symbolic description and the formula. — John Taylor, Jan 11 '21 at 09:15

Michael E2 · Accepted Answer · 2021-01-07T18:57:14.360

4

You can use Integrate to accumulate an InterpolatingFunction:

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

Addendum:

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)

edited Jan 07 '21 at 18:57

answered Jan 07 '21 at 18:33

Michael E2

235,386
17
334
747

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"] ] – Michael E2 Jan 07 '21 at 19:21

score 3 · Answer 2 · answered Jan 07 '21 at 18:26

Find the domain over which the values are not zero:

minmaxes = MinMax /@ Transpose@Select[DataInt, #[[3]] != 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 <= #[[1]] <= 2.1 && 0 <= #[[2]] <= 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. — Michael E2, Jan 07 '21 at 18:42

score 3 · Answer 3 · edited Jan 07 '21 at 20:10

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.

score 2 · Answer 4 · answered Jan 07 '21 at 21:26

$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.) — Michael E2, Jan 07 '21 at 21:29

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

4 Answers4