2022-02-06

# Uncertainty Propagation

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 derived variables.

# Purpose/Problem Statement

The idea is to explore how uncertainty in numerical values can propagate within equations and how do those uncertainties interact. I’ll explore some rock mechanics and calculate some 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 of the variables. What kind of uncertainty will the 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 would need. In my opinion 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

$\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:

$\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:

$\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:

$\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:

$\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?

In this section we’ll explore how the 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:

$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 we can understand how the uncertainty propagates.

## Code Preamble

The following sections are a mixture of text and code. We are going to 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 $$\left( \rho \right)$$

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

$\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³
rho = 2200.00+/-100.00 kg/m³

cv = 0.1/2.2 = 4.5%

## P-Wave Velocity$$\left( V_p \right)$$ and S-Wave Velocity$$\left( V_s \right)$$

The acoustic properties of the rock, the p-wave and s-wave velocity 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.:

$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$$\left( G \right)$$

$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$$\left( \nu \right)$$

$\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$$\left( E \right)$$

$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 formula for the modulus values are relatively straight forward. From the results we can note some interesting observations. Given the following table summarizing 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 summarizing 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 for the most part, 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 in a linear manor. Intuitively thinking that summing the coefficients yields anything meaningful is false. Without doing the calculations we may be lead 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 in a non-linear way. 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, it would mean that calculating Poisson’s Ratio is less reliably 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 results of modeling may not be accurate because we provided our model with good measured data.

NOTE: Observe that none of these equations have more than two independent variables.

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

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 that the 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 a tremendous amount of 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.