# Estimate errors during integration of interpolated data

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?

• If you approximate $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 between $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. For instance, Simpson's Rule has interpolatory order $2$ and DOP $3$. A Gaussian-Legendre rule with interpolatory order $n$ has DOP $2n+1$.... Commented Aug 5, 2020 at 13:22
• ...For Interpolation, you're using a different polynomial for each subinterval, so $K$ and $f^{(1+d)}(\xi)$ can vary. The error is just the sum over the subintervals of the errors given above. $K$ really just depends on the interpolation nodes that determine $p(x)$, so it is often the same for regular interpolation schemes. For the trapezoidal rule (your ord = 1), the DOP is $1$ and $K = -1/12$ (on all subintervals). Often the different derivative values $f^{(1+d)}(\xi)$ can be replaced by single one with $\xi$ somewhere in the integration interval.... Commented Aug 5, 2020 at 13:29
• ...This can be done for standard integration rules like the trapezoidal rule, for instance. For ord = 1, you have different values of $K$ for the first and last subinterval than for the interior subintervals. This is because Interpolation has different interpolation schemes for the exterior subintervals. Probably the derivative values $f^{(1+d)}(\xi)$ can be replaced by single one, but the lack of a consistent interpolation scheme probably makes the analysis more difficult. I have neither done it nor seen it done. In any case this seems more a mathematics question than a Mathematica Q. Commented Aug 5, 2020 at 13:37
• I was careless in my first comment: $p(x)$ should be an interpolatory polynomial, interpolating $f(x)$ at some nodes $x_0,\dots,x_n$. Then $\xi$ in $f^{(1+d)}(\xi)$ will lie in the interval containing the nodes $x_j$ and the end points $a$ and $b$. Commented Aug 5, 2020 at 13:44
• I should point out that 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. Commented Aug 5, 2020 at 13:56

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.

• If the data is noisy anyway, the OP has no business trying out high-order polynomial interpolation or anything derived from it, as it will only make things numerically worse. Commented Aug 5, 2020 at 23:51