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

Question

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.

Answer 1

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.

Answer 2

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["Integration`NIntegrateUtilities`"]

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]

Answer 3

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!

Answer 4

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  *)