We will use a simple example of a formula for a physical quantitity $g$ which depends on $x$, $y$, and $z$:
$$ f(x,y,z) = \frac{xy}{z^2}. $$Let's assume that we have measured the following values:
import numpy as np
from scipy import stats
def f(x,y,z):
'''Model function for error propagation examples'''
return x*y/z**2
x = 1.
alpha_x = 0.05
y = 2.
alpha_y = 0.07
z = 3.
alpha_z = 0.12
best = f(x,y,z)
print(best)
u to hold the values $(\alpha_f)_x$, $(\alpha_f)_x$, and $(\alpha_f)_z$.¶u = np.zeros(3)
u[0] = f(x+alpha_x, y, z)
u[1] = f(x, y+alpha_y, z)
u[2] = f(x, y, z + alpha_z)
u = u-best
unc = np.sqrt(u@u)
print(unc)
print("Fractional uncertainties", u/best)
When we write something like $x = 1.00 \pm 0.05$ it means that we believe that $x$ is a random variable selected from a normal distribution with a mean of $1.00$ and a standard deviation of $0.05$. In this technique we simulate very many experiments with values drawn from the appropriate distributions.
nEx = 10000 #Number of simulated experiments
x_sim = stats.norm.rvs(x,alpha_x,nEx) # array for simulated values of x
y_sim = stats.norm.rvs(y,alpha_y,nEx) # array for simulated values of y
z_sim = stats.norm.rvs(z,alpha_z,nEx) # array for simulated values of z
f_sim = f(x_sim, y_sim, z_sim)
print(np.mean(f_sim), np.std(f_sim))
In this example the partial derivatives are easy:
$$ \frac{\partial f}{\partial x} = \frac{y}{z^2}\quad\quad \frac{\partial f}{\partial y} = \frac{x}{z^2}\quad\quad \frac{\partial f}{\partial y} = -\frac{2xy}{z^3}. $$alpha_f_x = alpha_x*y/z**2
alpha_f_y = alpha_y*x/z**2
alpha_f_z = -2*alpha_z*x*y/z**3
alpha_f = np.sqrt(alpha_f_x**2 + alpha_f_x**2 + alpha_f_z**2 )
print(alpha_f)
version_information is from J.R. Johansson (jrjohansson at gmail.com)
See Introduction to scientific computing with Python:
http://nbviewer.jupyter.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-0-Scientific-Computing-with-Python.ipynb
for more information and instructions for package installation.
If version_information has been installed system wide (as it has been on Bucknell linux computers with shared file systems), continue with next cell as written. If not, comment out top line in next cell and uncomment the second line.
%load_ext version_information
%version_information numpy, scipy