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

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[#, InternalRealValuedNumericQ] &),
f_List?(VectorQ[#, InternalRealValuedNumericQ] &)
] :=
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.

• Using NIntegrate over an ElementMesh as done here, does a simple integration rule without adaptive refinement. – Michael E2 Jan 7 '19 at 22:40
• @MichaelE2 any idea how fast that can be? I'd be totally down to just let NIntegrate do it all for me but I might need to this to be fast in the future. – b3m2a1 Jan 7 '19 at 22:42
• @MichaelE2 follow up: can that be done without Interpolation or is that basically a necessity in methods like this? – b3m2a1 Jan 7 '19 at 22:43
• @b3m2a1 By "grid", do you mean a tensor product grid? – Henrik Schumacher Jan 7 '19 at 22:49
• @HenrikSchumacher unfortunately no. I have points in space that won't necessarily have come from a tensor product...although that is my dominant use case (although with different mesh spacings in different dimensions) so if you have something clean for that I'd be happy to see that too. – b3m2a1 Jan 7 '19 at 22:51

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.

• I see. Any hope if I can't have a uniform grid? – b3m2a1 Jan 7 '19 at 23:14
• If you know that data belongs to a function with bounded 4th derivative, you may obtain better accuracy with Simpson's 3/8 rule. – Henrik Schumacher Jan 7 '19 at 23:17
• For unstructured grid upto, 3 dimension, one may employ also a "trapezoidal rule" by (i) computing DelaunayMesh and employing the finite-element mass matrix for piecewise-linear functions on that mesh. But I am not aware of any integration rules with higher precision on totally unstructured grids. And for unstructured grids of dimension $d \geq 4$, I am totally clueless. – Henrik Schumacher Jan 7 '19 at 23:33