Estimate errors during integration of interpolated data

Question

I am interested into calculating the errors associated with an integration over interpolated data. Suppose I have something like this:

inc=Pi/10;
ord=1;
data=Table[N[{x,Sin[x]}],{x,0,Pi,inc}];
Integrate[Interpolation[data,InterpolationOrder->ord][x],{x,0,Pi/2}]

The output is 0.991762, so that is some error associated with the inc size and the interpolation order. Making ord=3, I get 0.999891. Is there a way to calculate such error? I looked up this question, but it's for NIntegrate and I'm more interested on the interpolation effects on the final result.

Edit: The data that I am trying to use comes from an experiment, so there's no analytic expression to go by. I just need to estimate how accurate is the result of Integrate[Interpolation[data,InterpolationOrder->ord][x],{x,a,b}] given the length of data and the interpolation order.

Edit 2: Maybe it could be done by doing a Spline fit, finding the errors for each point and then making a numerical integral propagating those errors?

Answer 1

Here is a method for estimating the error of integration. By sampling at two different rates, one can compare the integrals and estimate the error from the degree of precision of the integral (equals the InterpolationOrder in the cases at hand). If we interpolate $f(x)$ by a polynomial $p(x)$, then the error $\int_a^b f(x)\,dx - \int_a^b p(x) \, dx$ is given by $$ K \, (b-a)^{2+d} f^{(1+d)}(\xi) \, , $$ where $K$ is a constant depending on $p$, $\xi$ is in the interval containing the interpolation nodes, $a$ and $b$, and $d$ is the degree of precision (DOP, sometimes also called the order). The DOP $d$ is at least the interpolatory order ord. In the cases at hand, $d$ equals ord.

In the case of piecewise polynomial interpolation over uniformly spaced nodes, the error is the sum of such errors above. It will probably have a form like the following (I don't know for a fact in all cases): $$ K^* \, (b-a)\,h^{1+d} f^{(1+d)}(\xi) \, , $$ where $h$ is the spacing between the nodes, the same as inc in the OP's code. Here $K^*$ will depend on $h$ because of the way Interpolation works. A question I didn't bother to investigate is whether there is a $\xi$ such that the formula for the error is exactly equal to the error of integration. It may be well known, but I don't know it. In most cases either there is a $\xi$ or there is one that gets close enough.

If $f^{(1+d)}(\xi)$ does not vary in magnitude too much, then the error for a spacing of $2h$ will have an error approximately $2^{1+d}$ times the error for a spacing of $h$. The difference of the integrals therefore would be $2^{1+d}-1$ times the error of the $h$ integral.

The condition, "$f^{(1+d)}(\xi)$ does not vary in magnitude too much," is important. If the sampling rate is too low ($h$ is too big), then as it increases and $h$ decreases, the error will bounce around for while. After that, we enter a "convergent phase," in which the error decreases as $h^{1+d}$. For the estimate to be accurate we should be in the convergent phase. The bouncing around is explained by the formula in terms of the derivative factor varying greatly in magnitude. Often it can also be explained by offsetting errors or symmetry in the integrals over the subintervals. (This is a well-known feature of the trapezoidal rule on periodic functions.) In the pre-convergent phase, the estimate below may be unreliable.

Example that works

inc = Pi/10;
ord = 1;
data = Table[N[{x, Sin[x]}], {x, 0, Pi, inc}];
ifn = Interpolation[data, InterpolationOrder -> ord];
ifn2 = Interpolation[data[[;; ;; 2]], InterpolationOrder -> ord];
(Integrate[ifn2[x], {x, 0, Pi/2}] - 
    Integrate[ifn[x], {x, 0, Pi/2}]) /
  (2^(1 + ord) - 1) // ScientificForm
(* for comparison purposes *)
Integrate[Interpolation[ifn[x], {x, 0, Pi/2}] - 1 // ScientificForm

Pre-convergent-phase example

Change the order, change the interpolation, and you change the convergence behavior. The spacing inc = Pi/10 cannot really be considered small, and perhaps the previous example was lucky. It's not good enough for order 3.

inc = Pi/10;
ord = 3;
ifn = Interpolation[data, InterpolationOrder -> ord];
ifn2 = Interpolation[data[[;; ;; 2]], InterpolationOrder -> ord];
(Integrate[ifn2[x], {x, 0, Pi/2}] - 
    Integrate[ifn[x], {x, 0, Pi/2}]) /
  (2^(1 + ord) - 1) // ScientificForm
(* for comparison purposes *)
Integrate[Interpolation[ifn[x], {x, 0, Pi/2}] - 1 // ScientificForm

The estimate is a couple of orders of magnitude less than the actual error.

It's a very symmetric example. Maybe if we skip one of the interpolation nodes, we won't get the cancellation:

(Integrate[ifn2[x], {x, data[[3, 1]], Pi/2}] - 
    Integrate[ifn[x], {x, data[[3, 1]], Pi/2}]) /
  (2^(1 + ord) - 1) // ScientificForm
(* for comparison purposes *)
Integrate[ifn[x], {x, data[[3, 1]], Pi/2}] - 
  Integrate[Sin[x], {x, data[[3, 1]], Pi/2}] // ScientificForm

That's a much closer estimate (but it estimates the error of a slightly different integral). I assume one cannot (easily) change the sampling rate of experimental data. So one might drop the first or the last point and compute an error for each, and take the maximum to the overall estimate.

Dependence on the sampling rate

Just to show the effect of the sampling rate, we'll double it:

inc = Pi/20;
ord = 3;
ifn = Interpolation[data, InterpolationOrder -> ord];
ifn2 = Interpolation[data[[;; ;; 2]], InterpolationOrder -> ord];
(Integrate[ifn2[x], {x, 0, Pi/2}] - 
    Integrate[ifn[x], {x, 0, Pi/2}]) /
  (2^(1 + ord) - 1) // ScientificForm
(* for comparison purposes *)
Integrate[Interpolation[ifn[x], {x, 0, Pi/2}] - 1 // ScientificForm

That estimate is a much closer to the error.

Other considerations

• The higher the interpolatory order, the less stable. See the Runge phenomenon and Newton-Cotes rules.

• Likewise, higher-order derivatives are harder to estimate from data.

• Probably the most common choices are orders 1 and 3 and Method -> "Spline" (cubic splines).

• Noise in experimental data might be more significant than interpolation or integration error. Integration is generally less sensitive to perturbation than differentiation. But I'm unfamiliar with any theory about the error of numerical integration for perturbed data.