Uncertainty Propagation in Calculations#

In science and engineering, uncertainties and errors are a fact of life. This post is a study of how uncertainties can be used in calculations. More importantly, this post explores how uncertainty is propagated to subsequent calculations. That is, given a series of calculations that build on top of one another, what happens to the uncertainty?

Purpose/Problem Statement#

The idea is to explore how uncertainty in numerical values can propagate within equations and how those uncertainties interact. I’ll explore some rock mechanics and calculate rock properties from other measured properties. We’ll assume that enough testing was done to calculate averages and standard deviations from the collected data. Given a density, \(\rho\), p-wave velocity, \(V_p\), and s-wave velocity \(V_s\), with an amount of uncertainty tied to each variable. What kind of uncertainty will Young’s modulus, \(E\), Shear modulus, \(G\), and Poisson’s ratio, \(\nu\), have if they are calculated from these values?

Uncertainties Library#

We are going to use the uncertainties library. We don’t want to use this library to replace all of our calculations. It is rather slow and doesn’t have support for everything we need. It is excellent for understanding what is/what could be happening. The idea is to explore the relationship between measured quantities with uncertainties. How the calculated/estimated values are affected, and the magnitude of the propagated uncertainty is.

Linear Approximation of Error Propagation Theory#

The uncertainties library uses linear approximations from error propagation theory. See this wikipedia article for details. The linear approximation derivations can be found here. It is a good overview and introduction to the programming complexity that the library deals with.

Basic Example - Uncertainties Library#

This section explores some of the basics from the uncertainties library to demonstrate its usage.

import uncertainties
from uncertainties import ufloat

x1 = ufloat(10, 1)
x2 = ufloat(10, 0.5)

print(f'{x1=}')
print(f'{x2=}')

    x1=10.0+/-1.0
    x2=10.0+/-0.5

print(f'x1 + x2 = {x1 + x2:.3u}')
print(f'x1 - x2 = {x1 - x2:.3u}')

    x1 + x2 = 20.00+/-1.12
    x1 - x2 = 0.00+/-1.12

Uncertainty propagates addition and subtraction by a quadrature:

(1)#\[\sigma_f = \sqrt{a^2\sigma_A^2 + b^2\sigma_B^2 \pm 2ab\sigma_{AB}}\]

Note

The covariance, \(\sigma_{AB} = 0\) as the variables are uncorrelated simplifying the quadrature too:

(2)#\[\sigma_f = \sqrt{a^2\sigma_A^2 + b^2\sigma_B^2}\]
print(f'x1*x2 = {x1*x2:.3u}')
print(f'x1/x2 = {x1/x2:.3u}')
print(f'x1/x1 = {x1/x1:.3u}')

    x1*x2 = 100.0+/-11.2
    x1/x2 = 1.000+/-0.112
    x1/x1 = 1.0+/-0

Uncertainty propagates multiplication and division by a quadrature with relative uncertainties:

(3)#\[\sigma_f \approx \left| f \right| \sqrt{ \left(\frac{\sigma_A}{A}\right)^2 + \left(\frac{\sigma_B}{B}\right)^2 \pm 2\frac{\sigma_{AB}}{AB} }\]

Note

The covariance, \(\sigma_{AB} = 0\) as the variables are uncorrelated simplifying the quadrature too:

(4)#\[\sigma_f \approx \left| f \right| \sqrt{ \left(\frac{\sigma_A}{A}\right)^2 + \left(\frac{\sigma_B}{B}\right)^2 }\]
print(f'x1**2 = {x1**2:.3u}')
print(f'x1*x1 = {x1*x1:.3u}')

    x1**2 = 100.0+/-20.0
    x1*x1 = 100.0+/-20.0

Uncertainty propagates exponentiation by:

(5)#\[\sigma_f \approx \left| {a}{b}{A}^{b-1}{\sigma_A} \right| = \left| \frac{{f}{b}{\sigma_A}}{A} \right|\]

How does uncertainty Propagate with Complex Calculations?#

This section explores how uncertainty is propagated in complex calculations. We’ll use the uncertainties library to model the calculations and Pint to handle the units. We’ll analyze a rock with a given density, p-wave and s-wave velocity. From these values, we’ll calculate:

To understand how the uncertainties can accumulate from the different calculations, we’ll use the idea of the coefficient of variation:

(6)#\[c_v = \frac{\sigma}{\mu}\]

If we want to apply the coefficient of variation to a density, say \(3.25 \pm 0.2 \frac{g}{cc}\). we’ll assume that \(\mu = 3.25\) and \(\sigma = 0.2\). Thus \(c_v = \frac{\sigma}{\mu}\) and we have a measure to understand how the uncertainty propagates.

Code Preamble#

The following sections are a mixture of text and code. We will use Pint to handle the units.

import pint
u = pint.UnitRegistry()

# define the units we'll use.
density_unit = u.g/u.cm**3
velocity_unit = u.m/u.s

Density#

The density of the rock, \(\rho\), is a fundamental property used in estimating all other elastic modulus values and acoustic properties. Typically it is measured in terms of mass per unit volume, i.e.:

(7)#\[\rho = \frac{g}{\text{cm}^3} \text{or} \frac{kg}{m^3}\]
rho = ufloat(2.2, 0.1)*density_unit

print(f'rho = {rho:.2f~P}')
print(f'rho = {rho.to(u.kg/u.m**3):.2f~P}')
print()
print(f'cv = {rho.std_dev}/{rho.nominal_value} = {rho.std_dev/rho.nominal_value:.1%}')

    rho = 2.20+/-0.10 g/cm^3
    rho = 2200.00+/-100.00 kg/m^3

    cv = 0.1/2.2 = 4.5%

P-Wave Velocity and S-Wave Velocity#

The acoustic properties of the rock, the p-wave and s-wave velocities, are vital for calculating the dynamic modulus values. They can be measured in situ or estimated from the modulus values. They are usually measured in terms of distance over time, i.e.:

(8)#\[v = \frac{m}{s}\]
Vp = ufloat(6000, 325)*velocity_unit
Vs = ufloat(3800, 289)*velocity_unit

print(f'Vp = {Vp:.2f~P} ({Vp.std_dev/Vp.nominal_value:.1%}) (m/s)')
print(f'cv = {Vp.std_dev:.2f}/{Vp.nominal_value:.2f} = {Vp.std_dev/Vp.nominal_value:.1%}')
print()
print(f'Vs = {Vs:.2f~P} ({Vs.std_dev/Vs.nominal_value:.1%}) (m/s)')
print(f'cv = {Vs.std_dev:.2f}/{Vs.nominal_value:.2f} = {Vs.std_dev/Vs.nominal_value:.1%}')

    Vp = 6000.00+/-325.00 m/s (5.4%) (m/s)
    cv = 325.00/6000.00 = 5.4%

    Vs = 3800.00+/-289.00 m/s (7.6%) (m/s)
    cv = 289.00/3800.00 = 7.6%

Shear Modulus#

Shear modulus \(\left( G \right)\):

(9)#\[G = \rho V_s^2\]

Where:

  • \(G\) - Shear Modulus \(\left(Pa \right)\)

  • \(\rho\) - density Shear Modulus \(\left(\frac{g}{cc} \right)\)

  • \(V_s\) - S-Wave Velocity \(\left(\frac{m}{s} \right)\)

G = rho*Vs**2

print('G = rho*Vs**2')
print(f'G = {G.to(u.Pa):.2f~P}')
print(f'G = {G.to(u.GPa):.2f~P}')
print()
print(f'cv = {G.std_dev:.2f}/{G.nominal_value:.2f} = {G.std_dev/G.nominal_value:.1%}')

    G = rho*Vs**2
    G = 31768000000.00+/-5043226459.96 Pa
    G = 31.77+/-5.04 GPa

    cv = 5043226.46/31768000.00 = 15.9%

Poisson’s Ratio#

Poisson’s Ratio \(\left( \nu \right)\):

(10)#\[\nu = \frac{V_p^2 - 2 V_s^2}{2 \left(V_p^2 - V_s^2 \right)}\]

Where:

  • \(\nu\) - Poisson’s Ratio (Unitless)

  • \(V_p\) - P-Wave Velocity \(\left(\frac{m}{s} \right)\)

  • \(V_s\) - S-Wave Velocity \(\left(\frac{m}{s} \right)\)

nu = (Vp**2 - 2*Vs**2)/(2*(Vp**2 - Vs**2))

print('')
print('nu = (Vp**2 - 2*Vs**2)/(2*(Vp**2 - Vs**2))')
print(f'nu  = {nu:.4u}')
print()
print(f'cv = {nu.std_dev:.3f}/{nu.nominal_value:.3f} = {nu.std_dev/nu.nominal_value:.1%}')

    nu = (Vp**2 - 2*Vs**2)/(2*(Vp**2 - Vs**2))
    nu  = 0.1651+/-0.1044 dimensionless

    cv = 0.104/0.165 = 63.2%

Young’s Modulus#

Young’s modulus, \(E\):

(11)#\[E = 2 G \left(1 + \nu \right)\]
E = 2*G*(1 + nu)

print('E = 2*G*(1 + nu)')
print()

print(f'E  = {E.to(u.Pa):.2f~P}')
print(f'E  = {E.to(u.GPa):.2f~P}')
print()

print(f'cv = {E.std_dev:.3f}/{E.nominal_value:.3f} = {E.std_dev/E.nominal_value:.1%}')

    E = 2*G*(1 + nu)

    E  = 74027102040.82+/-7773579654.28 Pa
    E  = 74.03+/-7.77 GPa

    cv = 7773579.654/74027102.041 = 10.5%

Summary#

The calculations and formulae for the modulus values are relatively straightforward. From the results, we can note some interesting observations. The following table summarizes the input data:

Variable

Unit

Nominal Value

Standard Deviation

Coefficient of Variation (%)

Density

g/cc

2.2

0.1

4.5

P-wave

m/s

6000

325

5.4

S-Wave

m/s

3800

289

7.6

And the following table summarizes the calculations:

Variable

Unit

Nominal Value

Standard Deviation

Coefficient of Variation (%)

Shear Modulus

GPa

31.77

5.04

15.9

Young’s Modulus

GPa

74.03

7.77

10.5

Poisson’s Ratio

0.1651

0.1044

63.2

Notice that the calculated modulus values have a reasonable coefficient of variation (at least, I think so). Their standard deviations are relatively small. For example, the Shear modulus involves the density (4.5%) and the s-wave velocity (7.6%), so a 15.9% coefficient seems reasonable. This is the problem with uncertainties. They don’t combine and interact linearly. Intuitively thinking that summing the coefficients yields anything meaningful is false. We need to do the calculations to avoid being led to a false conclusion that low standard deviations (low errors/uncertainties) are what we are striving for.

For Poisson’s Ratio, the uncertainty is far larger than I would expect (63.2%). That is quite an uncertainty! If we look at the equation for Poisson’s ratio, it is clear that the uncertainty will interact and evolve non-linearly. It is quite complex compared to the other equations as it involves division, multiplication and exponentiation, and each of these contributes to the uncertainty in a non-linear way. Physically, calculating Poisson’s Ratio is less reliable than either of the modulus values from the measured data. Unfortunately, if that is the data we have, we have no choice but to use it. At least we can understand the risks associated with it and understand the modelling results may need to be more accurate because we provided our model with good measured data.

Note

Observe that all of these equations have at most two independent variables.

This large uncertainty for Poisson’s Ratio would mean that we should explore other, less complex, relationships that can estimate Poisson’s Ratio or physically measure the value from the system to reduce noise and eliminate it as a source of error. I would push hard to physically measure Poisson’s Ratio from rock core samples, along with Young’s modulus and static compressive strength values. Poisson’s Ratio and Young’s modulus are typically key values recorded from rock core samples. They would provide good input sources for other calculations if carefully collected and measured.

Another interesting observation is Young’s Modulus. It is calculated from the Shear modulus and Poisson’s Ratio, which, in turn, were calculated from the measured data. I would have assumed Young’s modulus would have had a much higher uncertainty than Poisson’s ratio. But it doesn’t—another counter-intuitive result.

This has been an interesting and eye-opening set of calculations and can shed tremendous light on sources of errors in your models. We should be able to apply this analysis to other variables and determine how uncertainty plays a key role. We cannot assume to understand how uncertainties will propagate without doing the calculations.

References#

Fuel Tracker   Windows Terminal and Cmder