# Multidimensional NIntegrate with Interpolating

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. 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@"NDSolveFEM";
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$$.

• If I specify rectangle Rectangle[{-1.04065, -1.04155}, {0.001, -0.148707}], ElementMeshInterpolation gives error about interpolation on unstructured grids. How can it happen? I assume that rectangle is structured grid.. Jun 29, 2020 at 9:38
• @ArtemAlexandrov I get no such error with this: Needs@"NDSolveFEM"; emesh = ToElementMesh[ Rectangle[{-1.04065, -1.04155}, {0.001, -0.148707}]]; fIFN = ElementMeshInterpolation[{emesh}, f @@@ emesh["Coordinates"]]; NIntegrate[fIFN[x, y], {x, y} \[Element] emesh, Method -> "FiniteElement"]. Did you use Interpolation instead of ElementMeshInterpolation? Jun 29, 2020 at 13:25
• Can this error occurs due to nature of my function f or should this approach work for any function? Jun 29, 2020 at 16:04
• @ArtemAlexandrov I don't think ElementMeshInterpolation gives such an error. Interpolation does give such an error. It doesn't have anything to do with the function f. If you can't share the code you used, then I'm guessing at what you're doing and I may be wrong. Jun 29, 2020 at 16:28
• Sorryy! I have made the silliest mistake. Everything works perfectly, thank you! Jun 29, 2020 at 18:05

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.

• It seems that argument of Interpolation should be given in form { {x,y}, f } ,not {x, y, f}` Jun 29, 2020 at 7:21