# 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 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


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?#

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:

(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 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#

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.:

(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 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.:

(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 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.

## References#

Previous: Fuel Tracker