How to integrate numerically a function with error bars?

Question

Typically, the function that one wants to integrate numerically, $f$, is given, i.e. its values for various points $\{x_i\}$ are known precisely. The resulting error is due to the fact that we chose a finite number of points. Let's call this error $numerical$.

I have a different situation. I know the values of $f$, $\{f(x_i)\}_{i=1}^N$, but with corresponding variances, $\{var(f(x_i))\}$. If I use, for example, the trapezoidal rule: $$F(b)-F(a)=\int_a^b f(x) dx \approx \frac{\Delta x}{2} [f(x_1) + 2 f(x_2) + \ldots 2f(x_{N-1}) + f(x_N) ] $$ (assuming that $x_i$ are equally separated by $\Delta x$), I calculate the variance of this estimate as: $$\left(\frac{\Delta x}{2}\right)^2 \left[ var(f(x_1)) + 4\cdot var(f(x_2)) + \ldots + 4\cdot var(f(x_{N-1})) + var(f(x_N)) \right] $$ but what about the $numerical$ error?

How can I combine these two errors? What is the "final" error of my estimate?

EDIT:

I've calculated the variance the same way I would do it for a sum of $N$ independent random variables, $X_i$, with coefficients $a_i$. That is, if: $$ S :=\sum_{i=1}^{N}a_i X_i,\qquad \forall_{i\neq j}Cov(X_i,X_j)=0 $$ then $$ Var(S)=\sum_{i=1}^{N} a_i^2 Var(X_i). $$ If this error analysis is not correct, please let me know.

Also, I noticed a problem with the variance of partial sums, $S_k:=\sum_{i=1}^k a_i X_i$. It seems that the error propagates, and therefore $S_N$ has the largest variance, while $S_2$ -- the lowest. But if I were to reverse the order in the sum (which would correspond to estimating a reversed integral, $\int_b^a f dx$), the error would again propagate, but in the opposite direction. I would get two different estimates of "partial integrals", depending on the direction of integration.

Is there a way of combining information from numerical integration in both ways?

Answer 1

Assuming the integration is deterministic (that is, ignoring effects like bit flips and hardware issues that are generally negligible at scale), you have two types of uncertainties:

• the uncertainty in calculating the integral approximately with numerical methods instead of exactly; this error is epistemic, and with enough analysis, could be calculated or estimated

• the uncertainty in each of the function values sampled; unless there's additional relevant information, this uncertainty appears to be aleatoric

Since you have variances, one approach to quantifying the aleatoric uncertainty that is compatible with the information you've provided is polynomial chaos. You could integrate the mean and variance separately to obtain the mean and variance of the integral. This approach is different than the one you propose above; how you arrived at the formula above for the variance is unclear to me.

Another approach would use interval arithmetic, but this approach is better at capturing worst-case upper and lower limits, and it's not clear that such information would be valuable to you.

As for the remaining epistemic uncertainty, you could ignore it, if it's sufficiently small relative to the aleatoric uncertainty in the integral (and in your data points). If you feel you do need to quantify it, you'd need to estimate or bound the error somehow (estimates via embedded methods, arbitrary precision integration, or bounds from numerical analysis textbooks/research). As a first approximation, you could model the epistemic error in the integral of the mean as a uniform random variable at each point in space, calculate the resulting variance, and combine that variance information with the integral of the variance calculated from the aleatoric uncertainties, using a polynomial chaos approach. The error in integrating the variances could then be incorporated into calculations of higher-order moments of the random variable describing the integral you wish to calculate. Although this process could be continued with progressively higher-order moments, at some point, the moment expansion of this variable will need to be truncated in calculations.