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.:

   grid_List?(VectorQ[#, Internal`RealValuedNumericQ] &), 
   f_List?(VectorQ[#, Internal`RealValuedNumericQ] &)
   ] :=
    diffs = Differences[grid],
    means = MovingAverage[f, 2]

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.

  • $\begingroup$ Using NIntegrate over an ElementMesh as done here, does a simple integration rule without adaptive refinement. $\endgroup$
    – Michael E2
    Commented Jan 7, 2019 at 22:40
  • $\begingroup$ @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. $\endgroup$
    – b3m2a1
    Commented Jan 7, 2019 at 22:42
  • $\begingroup$ @MichaelE2 follow up: can that be done without Interpolation or is that basically a necessity in methods like this? $\endgroup$
    – b3m2a1
    Commented Jan 7, 2019 at 22:43
  • $\begingroup$ @b3m2a1 By "grid", do you mean a tensor product grid? $\endgroup$ Commented Jan 7, 2019 at 22:49
  • $\begingroup$ @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. $\endgroup$
    – b3m2a1
    Commented Jan 7, 2019 at 22:51

1 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]



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.

  • $\begingroup$ I see. Any hope if I can't have a uniform grid? $\endgroup$
    – b3m2a1
    Commented Jan 7, 2019 at 23:14
  • $\begingroup$ If you know that data belongs to a function with bounded 4th derivative, you may obtain better accuracy with Simpson's 3/8 rule. $\endgroup$ Commented Jan 7, 2019 at 23:17
  • $\begingroup$ 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. $\endgroup$ Commented Jan 7, 2019 at 23:33

Your Answer

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

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