# Choosing points of integration in NIntegrate (integrating a function given by a list without interpolation)

When using NIntegrate, say in two-dimensions, is it possible to specify the points of the grid Mathematica will use? For example, if the integrand was a function f[x,y] given as a list at specific points only of x and y and we don't want to Interpolate it over other points, is there a way to make Mathematica do something like

NIntegrate[f[x,y],{x,0,10},{y,0,10}]


using only the given known values of f on that list?

How does NIntegrate work in such cases? And if specifying grid points for it is possible, is there a common approach or rule of thumb to know when the grid is refined enough to give correct answer (within some error criterion, say 0.1%)? I suspect one wouldn't stop immediately after it first passes the criterion, but may want that to repeat a few times to be sure it converges. Any common practice for this?

The motivation behind this question is to learn how to numerically integrate a function given by a finite list, without interpolating it.

• Look for method "TrapezoidalRule" or use a linear interpolation of the gridpoints. Commented Nov 17, 2020 at 7:26
• Maybe the article "tutorial/NIntegrateOverview" in MMA's help will help you? Commented Nov 17, 2020 at 9:56
• You can use the function values to create an interpolation function and then take the primitive of that interpolation function without going through NIntegrate. This effectively applies the correct quadrature rules. See this answer, for example. Commented Nov 17, 2020 at 10:19
• @SjoerdSmit Thanks for this great comment! Commented Nov 17, 2020 at 22:34

Taking Alexei's example, here's a simple method using the Interpolation framework. I know you requested to not actually interpolate the function, but that will not be necessary. For example:

lst = Table[{x, x*Exp[-x] + 1}, {x, 0, 2, 0.2}];
primitive = Derivative[-1] @ Interpolation[lst]


The primitive computed this way essentially has all the quadrature information consistent with the InterpolationFunction generated from the list of points. You can directly use it to integrate the underlying function:

primitive[2] - primitive[0]
NIntegrate[x*Exp[-x] + 1, {x, 0, 2}]


2.59402

2.59399

Unfortunately, though, it seems like this only works with Hermite interpolation and not with Spline interpolation.

• Thanks for this advice. Can I ask: (1) why do you call it 'primitive'? (2) What is the difference between Hermite and Spline interpolation in practice and how can I tell during coding some function (any brief notes on this)? (3) Finally, I noticed that interpolation tends to always be fast for a well-behaved function (for a given size of a set of points), am I correct? For example, if I had a loop that basically has just one 'NIntegrate' and one Interpolation commands and it was slow, would NIntegrate be the usual suspect or can interpolation also be slow? If not so, how would you check? Commented Nov 17, 2020 at 22:12
• @user135626 Re (2): "In calculus, an antiderivative, inverse derivative, primitive function, primitive integral or indefinite integral..." -- synonyms. Commented Nov 18, 2020 at 3:02
• I don't have time right now for all questions, but for details on the interpolation methods I suggest looking here and here. The Method option of Interpolate controls the type of interpolation used (see documentation). The default is Hermite interpolation. As for your question about speed: interpolation is very fast in Mathematica. I've never seen a slow interpolation function. Commented Nov 18, 2020 at 9:17
• As for the question about terminology: personally I prefer the word "integral" to mean "area under a curve" and "primitive" as "a new function that gives the original one when you take its derivative". A primitive can be used for integration, but the two concepts are very different otherwise. Commented Nov 18, 2020 at 9:23
• @user135626 A search for "InterpolationPointsSubdivision" and "MaxSubregions" shows some Q&A that deal with problems NIntegrate sometimes has with interpolating functions. Commented Nov 18, 2020 at 16:28
1. I think Trapezoidal strategy is the closest NIntegrate built-in functionality to what you describe.

• You can modify the function TrapStep in the example Trapezoidal strategy implementation (in the linked advanced documentation section) to use points that you want.
2. You can also make and plug-in custom integration rules and integration strategies.

• If you want to use NIntegrate's framework...

Here are code and plot that illustrate point 1:

Needs["IntegrationNIntegrateUtilities"]

aPoints =
Table[pg ->
NIntegrateSamplingPoints[
NIntegrate[1/Sqrt[x + y + 1], {x, 0, 10}, {y, 0, 10},
Method -> {"Trapezoidal", "SymbolicProcessing" -> False},
PrecisionGoal -> pg, MaxRecursion -> 4]][[1, All, 1]], {pg, 1, 3}];

ListPlot[Values[aPoints],
PlotStyle ->
Table[PointSize[(Length[aPoints] - i + 1)*0.01], {i, 1,
Length[aPoints]}], PlotTheme -> "Detailed",
AspectRatio -> Automatic, PlotLegends -> Automatic,
ImageSize -> Large]


I do not really understand, why do you need to do as you described. However, one way would be to use the table to define a polygon, then transform it into Region and apply Area. This way, however, introduces errors related to the approximation of the intervals between the points by straight lines.

Let me give an example. Since you posted no table I make a simple one:

lst = Table[{x, x*Exp[-x] + 1}, {x, 0, 2, 0.2}]

(*  {{0., 1.}, {0.2, 1.16375}, {0.4, 1.26813}, {0.6, 1.32929}, {0.8,
1.35946}, {1., 1.36788}, {1.2, 1.36143}, {1.4, 1.34524}, {1.6,
1.32303}, {1.8, 1.29754}, {2., 1.27067}}  *)


looking as follows:

Now let us supplement it by the points lying on the axis to make a polygon:

lst1 = Append[
Append[Prepend[lst, {First[lst][[1]], 0}], {Last[lst][[1]],
0}], {First[lst][[1]], 0}]

(*  {{0., 0}, {0., 1.}, {0.2, 1.16375}, {0.4, 1.26813}, {0.6,
1.32929}, {0.8, 1.35946}, {1., 1.36788}, {1.2, 1.36143}, {1.4,
1.34524}, {1.6, 1.32303}, {1.8, 1.29754}, {2., 1.27067}, {2.,
0}, {0., 0}}  *)


Let us draw it:

Graphics[{Red, Polygon[lst1]}]


Now let us define it as a region and calculate its area:

r = Region[Polygon[lst1]];
Area[r]


(* 2.59022 *)

Let us now check, what the direct integration gives:

NIntegrate[x*Exp[-x] + 1, {x, 0, 2}]

(*  2.59399  *)


The difference is in the third position after point. This is, however, because the function was "good" and the points in the list lst were taken rather densely.

Have fun!

Direct programming of the trapezoidal rule:

Block[{f, xList, yList, fPoints},

(* Set up as indicated in OP with random data *)
f[x_, y_] := Sin[x^2/20] Exp[y/2];
xList = Prepend[0.]@Append[10.]@
Accumulate[
10 RandomVariate[DirichletDistribution[ConstantArray[1, 2000]]]];
yList = Prepend[0.]@Append[10.]@
Accumulate[
10 RandomVariate[DirichletDistribution[ConstantArray[1, 1500]]]];
fPoints =
f @@ Transpose[Outer[List, xList, yList], {2, 3, 1}]; (* 2d array *)

(* Trapezoidal rule *)
Fold[#2.(Most[#] + Rest[#])/2 &,
fPoints,
Differences /@ {xList, yList}]
]

(* 769.973  *)

NIntegrate[Sin[x^2/20] Exp[y/2], {x, 0, 10}, {y, 0, 10}]

(*  769.968  *)