How does uncertainty of dataset propagates through numerical integration?

Question

In the following text, uncertainty refers to standard deviation.

I have 500 time series which I use in a few equations and get averages and uncertainties. Through calculations for each time series, I end up with a data set with 500 values of the position z, the resistance $R(z)$, the absolute and the relative uncertainty: $$ z \,, R( z ) \,, \pm \Delta R(z) \,, \pm \frac{\Delta R(z)}{ R(z) }$$

What I want is to calculate the permeation $P$ through the following formula: $$ P=\frac{1}{\int_{z_{min}}^{z_{max}}R(z)dz} $$

therefore, I need to numerically integrate these 500 values. This is easy with a trapezoidal but the question is what happens to the uncertainty? How does it propagate through the integral? Essentially I want to report the average $ \left( P \pm \Delta P \right) \, [units]$

Answer 1

I found an answer to my question which in my opinion is correct (or at least has a point). Nevertheless it would be nice to have a validation from someone that might know more on the topic.

Since the average is calculated with a trapezoidal (which in reality is a summation) then the standard deviation can also be calculated by uncertainty propagation algebra for addition.

Based on Wikipedia for trapezoidal rule: $$ \int_c^d f(x) dx \approx (d-c) \left[ \frac{f(c)+f(d)}{2} \right] = \left( \frac{d-c}{2} \right) f(c) + \left( \frac{d-c}{2} \right) f(d)$$

which will follow the 2nd uncertainty propagation rule from Wikipedia: $$ f=aA+bB$$

with:

• $ a = b = \frac{d-c}{2}$
• $ A = f(c) $
• $ B = f(d) $

In terms of Python3 code that translates to:

import numpy as np

# FUNCTIONS
def uncAddition(A=1.,B=1.,sigmaA=0.,sigmaB=0.,a=1.,b=1.,sigmaAB=0.):
    f =  (a * A) + (b * B)
    sigmaf = np.sqrt( (a * sigmaA)**2 \
                    + (b * sigmaB)**2 \
                    + (2 * a * b * sigmaAB) )
    return f, sigmaf


# MAIN
# Data: is ab nx2 array that contains n number of f values 
#           and their absolute standard deviation
# x:    is the independent variable of the f values
# Obviously x and f should have the same number of rows

integralAve=0.
integralStD=0.
for i in range(1,len(data[:,0])):
    trapzPartAve,trapzPartStd = uncAddition(A=f[i-1,0],
                                            sigmaA=f[i-1,1],
                                            a=(x[i]-x[i-1])/2.0,
                                            B=f[i,0],
                                            sigmaB=f[i,1],
                                            b=(x[i]-x[i-1])/2.0) 
    integralAve += trapzPartAve
    integralStd += trapzPartStd

print(integralAve,integralStd)