1
$\begingroup$

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?

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – Michael E2
    Jan 7, 2021 at 18:49
  • 1
    $\begingroup$ 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]? $\endgroup$
    – Michael E2
    Jan 7, 2021 at 21:32
  • $\begingroup$ @MichaelE2 : I don't know the symbolic description and the formula. $\endgroup$ Jan 11, 2021 at 9:15

4 Answers 4

4
$\begingroup$

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)

$\endgroup$
1
  • $\begingroup$ 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"] ] $\endgroup$
    – Michael E2
    Jan 7, 2021 at 19:21
3
$\begingroup$

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 *)
$\endgroup$
1
  • 2
    $\begingroup$ Keep in mind that the OP's interpolating function will perform a linear extrapolation (InterpolationOrder -> 1) beyond the boundary. $\endgroup$
    – Michael E2
    Jan 7, 2021 at 18:42
3
$\begingroup$

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

plot of the region where function is not zero

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.

$\endgroup$
2
$\begingroup$
$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 *)
$\endgroup$
1
  • 2
    $\begingroup$ 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.) $\endgroup$
    – Michael E2
    Jan 7, 2021 at 21:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.