# How to build a grid of integrand points and numerically integrate?

If I have some function I know numerically only, say f(x) and each point $x$ takes significant time to compute so I have them all stored in some file as f(1)=0.232423, f(1.1)=0.3243432,....Then is it possible to tell NIntegrate only sample on the grid x=1,1.1,1.2,1.3...... and try to integrate this function?

Then is the best method to keep refining my grid (gathering more and more data in my file) until the result of integration stops changing?

• you could try to ListInterpolate
– user21
Oct 12, 2012 at 10:24
• ...or you could directly implement the trapezoidal rule on your function. Oct 12, 2012 at 11:22
• Or perhaps you could guess a good fit. Can you upload a significant sample of your points? Oct 12, 2012 at 11:23

The question was if we can tell NIntegrate to only use discrete points on a grid, and I believe the answer is no. But that's not a very useful answer, so one should look for an alternative, which will either involve interpolation, or summation instead of integration.

Although the comments point to some possible answers, I'll mention one approach here that is especially useful if you know the function on a large number of equally spaced grid points and moreover have some derivatives of the function at the end points of the grid only.

Then (assuming the actual function can be Taylor expanded), one can estimate the error of a discrete approximation to the integral with the help of the BernoulliB numbers, going arbitrarily far beyond the trapezoidal approximation depending only on how many derivatives you can calculate. Since the latter is only needed at the end points, this effort may be small compared to the cost of evaluating the function on the rest of the grid.

Here I'll use a Bessel function as an example to generate the discrete data and compare the integral to the sum approximation:

xMax = 15.;
nMax = 100;
\[CapitalDelta]x = xMax/nMax;
Clear[n]

ftest[n_] = BesselJ[3, n \[CapitalDelta]x]

(* ==> BesselJ[3, 0.15 n] *)

points = Map[ftest, Range[0, nMax]];

Show[Plot[ftest[n], {n, 0, nMax}],
ListPlot[Transpose[{Range[0, nMax], points}]]] The points indicate the discrete values I'll use in the sum, and the horizontal axis has been scaled so that the integral or sum will extend from n = 0 to n = 100 (the sum will be in steps of 1).

Next is the application of the BernoulliB numbers. I'll assume we can calculate the function ftest and its first three derivatives at n = 0 and n = 100. The variable r will then contain a general formula for the error of the sum compared to the integral:

Clear[f];
With[{numberOfDerivatives = 3},
r = Sum[D[Integrate[f[n], n], {n, k}] BernoulliB[k]/k!, {k, 0,
numberOfDerivatives + 1}] - \[Integral]f[n] \[DifferentialD]n]


$-\frac{f(n)}{2}+\frac{f'(n)}{12}-\frac{1}{720} f^{(3)}(n)$

error = With[{rest = r /. f -> ftest},
rest /. n -> nMax - rest /. n -> 0]

(* ==> 0.0980142 *)


The variabele error is the specific error for our function ftest and choice of interval for n.

Now we just have to compare the exact and approximate ways of calculating the integrals:

NIntegrate[ftest[n], {n, 0, nMax}]

(* ==> 6.20721 *)

Sum[ftest[n], {n, 0, nMax}]

(* ==> 6.1112 *)

% + error

(* ==> 6.20922 *)


The last line means that by adding the error term to the sum (previous line), we get very close to the integral (the line above that). Or in other words, the integral is approximated by

Sum[ftest[n], {n, 0, nMax}] + error


and the correction error is determined purely by ftest and its derivatives at n = 0 and n = 100.