Multidimensional NIntegrate with Interpolating

Question

I have a numerically evaluated function f[x,y] (it is impossible to write down analytical epxression for the function f) and array of points {x,y} with constant steps in x and y directions. I try to perform suggested strategy (see this question) My data is not pretty and Interpolation says me that InterpolationOrder should be set to All or to 1. If I use ->1, it returns

Interpolation::fememtlq: The quality -1.21606*10^-15 of the underlying mesh is too low. 
The quality needs to be larger than 0.`.

Then, I obtain interpolation function and try to perform NIntegrate. However, when I try to do this, it seems that Mathematica kernel crashes.

To be honest, I know nothing about numerical integration but I think that Method specification in NIntegrate can help. I can provide my data and result of "naive" calculation (just replace integrate by sum over points).

I try to understand which method of NIntegrate should be used. I plot DensityPlot of interpolating function and (may be) it helps.

Answer 1

You could use a simpler test function for the sake of illustration. Let's take the following function over the implied rectangle, which can integrated exactly for the sake of comparison with the approximations:

f = Function[{x, y}, Exp[2 x - y]];

N@Integrate[f[x, y], {x, 0, 2}, {y, 0, 3}]
(*  25.4648  *)

One approach, depending on whether the fixed-step rectangle grid is itself important, is to use an ElementMesh:

Needs@"NDSolve`FEM`";
emesh = ToElementMesh[Rectangle[{0., 0.}, {2., 3.}]]
(*  ElementMesh[{{0., 2.}, {0., 3.}}, {QuadElement["<" 425 ">"]}]  *)

fIFN = ElementMeshInterpolation[{emesh}, f @@@ emesh["Coordinates"]];

NIntegrate[fIFN[x, y], {x, y} ∈ emesh, 
 Method -> "FiniteElement"]
(*  25.4647  *)

To increase accuracy, use a finer mesh:

emesh = ToElementMesh[Rectangle[{0., 0.}, {2., 3.}], 
  MaxCellMeasure -> "Length" -> 0.01]
(*  ElementMesh[{{0., 2.}, {0., 3.}}, {QuadElement["<" 60000 ">"]}]  *)

fIFN = ElementMeshInterpolation[{emesh}, f @@@ emesh["Coordinates"]];

NIntegrate[fIFN[x, y], {x, y} ∈ emesh, 
 Method -> "FiniteElement"]
(*  25.4648  *)

Another approach is to use a regular interpolation and Integrate:

fIFN = FunctionInterpolation[Exp[2 x - y], {x, 0, 2}, {y, 0, 3}];

Integrate[fIFN[x, y], x, y] /. {x -> 2, y -> 3}
(*  25.4618  *)

To improve the result, use more interpolation points:

fIFN = FunctionInterpolation[Exp[2 x - y], {x, 0., 2.}, {y, 0., 3.}, 
   InterpolationPoints -> 101];

Integrate[fIFN[x, y], x, y] /. {x -> 2, y -> 3}
(*  25.4648  *)

Integrate[InterpolatingFunction[...][x,y], x, y] returns $\int_a^x\int_b^y f(x,y) \, dx \, dy$, where the domain of the interpolating function is of the form $a \le x \le c$, $b \le y \le d$.

Answer 2

First, make sure your data is in the following form data = {{x1,y1,z1},{x2,y2,z2},...,{xn,yn,zn}} Then you can create your interpolating function, where the interpolation will be linear. Something to be careful of is that if your grid spacings, e.g. x1 - x2, are not uniform across your whole grid, then you may run into trouble. Moreover, if your grid isn't square, you'll have problems. Make a square grid, and any region you don't want to integrate just set it to zero. With this in mind, define your interpolating function fn = Interpolation[data, InterpolationOrder->1] and integrate answer = NIntegrate[ fn[x,y],{x,x1,xn},{y,y1,yn},Method->"QuasiMonteCarlo"] I like the method QuasiMonteCarlo, since monte carlo tends to be fast, and the quasi random sequence of points generated in the quasimontecarlo routine is more uniform than the points chosen from the PRNG in the MonteCarlo method, so you're less subject to sqrt(N) noise in your answer.