import scipy as sp
from scipy import stats
import sympy as sym # for symbolic differentiation
from sympy.interactive import init_printing # provides LaTex formatted cell output
sym.init_printing()
In PHYS 211 we did an experiment to determine "little $g$" from measurement of the length and period of a pendulum. The following cell defines a function giving $g$ for known values of $l$ and $T$; the next cell it gives typical data and associated uncertainties.
def g(l,t):
return 4*sp.pi**2*l/t**2
mean_l = 1.0 # Measured length of pendulum in meters
sigma_l = 0.004 # Uncertainty in length - assumed to be
#i.d. of normal distribution
mean_t = 2.0; # Measured period
sigma_t = 0.005 # Uncertainty in period - assumed to be
# s.d. of normal distribution
sigma_gl = g(mean_l+sigma_l,mean_t)-g(mean_l,mean_t)
sigma_gt = g(mean_l,mean_t+sigma_t)- g(mean_l,mean_t)
sigma_g = sp.sqrt(sigma_gl**2 + sigma_gt**2)
g(mean_l,mean_t),sigma_gl,sigma_gt,sigma_g
Result: $g = 9.87\pm 0.06$
Redefine $g(L,T)$ in terms of sympy floats so that we can do symbolic work. Will also use sym.sqrt() below.
NOTE: sympy floats != regular python floats See http://docs.sympy.org/dev/gotchas.html Either sym.sqrt(sym.pi), or sp.sqrt(float(sym.py)), or sym.sqrt(sym.Float(sp.pi)), or sym.sqrt(sym.syimpify(sp.pi))
def g(l,t):
return 4*l*sym.pi**2/t**2
In sympy you must declare symbolic variable explicitly.
The method below allows you to control assumptions, as in
n = sym.symbols('n',postive=True)
You can also import directly from a set of common symbols, e.g.,
from sympy.abc import w
or
sym.var('z')
l, t = sym.symbols('l, t')
sym.diff(g(l,t),t,1),sym.diff(g(l,t),l)
Evaluate the symbolic expressions at the values $\bar{l}$ and $\bar{t}$, and calculate the uncertainty.
sigma_gt = sym.diff(g(l,t),l).evalf(subs={l:mean_l,t:mean_t})*sigma_t
sigma_gl = sym.diff(g(l,t),l).evalf(subs={l:mean_l,t:mean_t})*sigma_l
sigma_g = sym.sqrt(sigma_gl**2 + sigma_gt**2).evalf()
g(mean_l,mean_t).evalf(),sigma_gl,sigma_gt,sigma_g
Result: $g = 9.87 \pm 0.06$
This agrees with previous result.
Here we are going to get the uncertainty in a different way, using simulated data. Our propagation of uncertainty formulas assume that the uncertainties are standard deviations of normal distribution. We can "redo" the experment, with new simulated "measurements" sampled from the assumed distributions.
First we return to our definition of $g(L,T)$ that returns a regular python float.
def g(l,t):
return 4*sp.pi**2*l/t**2
nEx = 10000 # Number of simulated experiments
l = sp.random.normal(mean_l,sigma_l,nEx) # Pick values of l from normal distribution
t = sp.random.normal(mean_t,sigma_t,nEx) # Pick values of t from normal distribution
gg = g(l,t) # Calculate nEx values of g using random l's and t,s
gg
sp.mean(gg),sp.std(gg)
Result: $g = 9.87 \pm 0.06$
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
#%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
%version_information numpy, sympy