Trapezoid rule-esque integration with unevenly spaced n-D grid

Question

I have a function (or rather many functions) defined over an n-dimensional grid. The grid points will often but not always be evenly spaced. I want to compute integrals of these functions numerically.

My question now is what the best approach to this is. In 1-D I can always works with the trapezoid rule, e.g.:

chapeaux[
   grid_List?(VectorQ[#, Internal`RealValuedNumericQ] &), 
   f_List?(VectorQ[#, Internal`RealValuedNumericQ] &)
   ] :=
  Module[
   {
    diffs = Differences[grid],
    means = MovingAverage[f, 2]
    },
   diffs.means
   ];

And here are some simple tests:

gigi = Range[-10., 10., .1];
grigri = Sort@BlockRandom[RandomReal[{-10. , 10.}, 100*Length@gigi]];
sinf = Sin@gigi;
sinfr = Sin@grigri;

NIntegrate[Sin[x], {x, -10, 10}] // RepeatedTiming

{0.0049, 0.}

chapeaux[gigi, sinf] // RepeatedTiming

{0.000066, -8.32667*10^-17}

chapeaux[grigri, sinfr] // RepeatedTiming

{0.00070, 0.00164231}

harm = gigi^2;
harmr = grigri^2;

NIntegrate[x^2, {x, -10, 10}] // RepeatedTiming

{0.0034, 666.667}

chapeaux[gigi, harm] // RepeatedTiming

{0.000068, 666.7}

chapeaux[grigri, harmr] // RepeatedTiming

{0.000694, 666.281}

And sure in the unstructured case I need about 100X as many points to even get okay-ish accuracy, but that's not the end of the world, so now the question is how best to do this or something similar in n-dimensions. Of course, the functions I'm using have no analytic analog and so the best I can do is NIntegrate the Interpolation, but I'm hoping to do this faster and cleaner.

At some level this is mostly a math question, of course, but I figured it's worth asking here before implementing as the accumulated math knowledge here is much greater than my own.

Answer 1

I would suggest Tai's method! ;)

Okay, there are two issues here with the trapezoidal rule:

1. It is crucial that the interval endpoints are included in the grid.
2. Trapezoidal rule with uniform subdivision is exact of all odd functions just because the quadrature is symmetric. Unfortunately, you tested with Sin which is an odd function. This made the trapezoidal rule look way more potent than it is.

This is a fairer comparison:

f = Function[x, Sin[x] + x^2 + x/2]
gigi = Range[-10., 10., .1];
grigri = Join[{-10.}, Sort@BlockRandom[RandomReal[{-10., 10.}, Length@gigi - 2]], {10.}];
sinf = f@gigi;
sinfr = f@grigri;

a = NIntegrate[f[x], {x, -10, 10}];
b = chapeaux[gigi, sinf];
c = chapeaux[grigri, sinfr];

Abs[a - b]/Abs[a]
Abs[a - c]/Abs[a]

0.00005

0.000259904

The error of the trapezoidal ruled for the grid points $x_0,\dotsc,x_n$ should be be bounded by

$$\sup_{x\in I} |f''(x)| \, \cdot \sum_{i=1}^n |x_{i}-x_{i-1}|^2.$$

So uniform grids optimize this bound by minimizing the sum of squares of the subintervals.