In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.stats import chi2
14.01
In [2]:
x = np.random.uniform(0., 10., 100)
y = np.polyval([1, 2, -3], x) + np.random.normal(0., 10., 100)
e = np.random.uniform(5, 10, 100)
In [3]:
def poly1 (x,a0,a1):
return a0*x+a1
def poly2 (x,a0,a1,a2):
return a0*x**2+a1*x+a2
In [4]:
xp=np.linspace(0,10,100)
plt.errorbar(x,y,yerr=e, marker='o', markersize=3, capsize=3, linestyle='none')
popt,pcov = curve_fit (poly1,x,y,sigma=e, absolute_sigma=True)
plt.plot(xp,poly1(xp,popt[0],popt[1]))
popt,pcov = curve_fit (poly2,x,y,sigma=e, absolute_sigma=True)
plt.plot(xp,poly2(xp,*popt))
Out[4]:
[<matplotlib.lines.Line2D at 0x7ca2f3f765d0>]
14.02
In [5]:
def temp_seasonal(x,a0,a1,a2):
return a0*np.cos(2*np.pi*x+a1)+a2
In [6]:
# The following code reads in the file and removes bad values
import numpy as np
date, temperature = np.loadtxt('data/munich_temperatures_average_with_bad_data.txt', unpack=True)
keep = np.abs(temperature) < 90
date = date[keep]
temperature = temperature[keep]
popt,pconv = curve_fit(temp_seasonal,date,temperature)
plt.scatter(date,temperature,marker='.')
plt.plot(date,temp_seasonal(date,*popt),color='red',lw=3)
plt.xlim(2008,2012)
Out[6]:
(2008.0, 2012.0)
In [7]:
print ("temp=",popt[0],"*cos(2*pi*x+",popt[1],")+",popt[2],")")
temp= -9.955183583978071 *cos(2*pi*x+ 12.31340590740421 )+ 9.040845452201035 )
The average temperature is 9.0 degrees, the minimum/maximum are -3.3/21.4 degrees.
The meaning of popt[1]:
In [8]:
diff=(4*np.pi-popt[1])
print (365*diff/(2*np.pi))
14.695112991981393
This means that according to this model the coldest time of year is about 15 days after New Year.
15.01
In [9]:
from scipy.integrate import quad
from scipy.integrate import simpson
def gaussian(x,mu,sig):
return (1/(np.sqrt(2*np.pi)*sig))*np.exp(-(x-mu)**2/(2*sig**2))
In [10]:
res, _ = quad(gaussian,-1,1,args=(0,1)) # ignore uncertainty here
print (res)
0.682689492137086
In [11]:
x=np.linspace(-1,1,100)
y=gaussian(x,0,1)
simpson(y,x=x)
Out[11]:
np.float64(0.6826894964572258)
15.02
In [12]:
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
In [13]:
#critical damping
def diffeq(x, t, k, g):
#define/clear derivative
deriv = np.zeros(2)
deriv[0] = x[1]
deriv[1] = -k*x[0] - 2 * g * x[1]
return deriv
k=0.1
g=np.sqrt(k)
x = np.array([2.0,-8.2]) # initial position as before
time = np.linspace(0.0,100.0,100) # evaluation times
""" calling the solver """
solution = odeint(diffeq, x, time, args=(k,g))
"""plotting the result"""
plt.plot(time,solution[:,0])
Out[13]:
[<matplotlib.lines.Line2D at 0x7ca2f143d950>]
In [14]:
from scipy.integrate import odeint
# forced oscillation
def diffeq(x, t, k, g, a0, om0):
#define/clear derivative
deriv = np.zeros(2)
deriv[0] = x[1]
deriv[1] = -k*x[0] - 2* g*x[1] + a0*np.cos(om0*t)
return deriv
""" constants """
k=0.1
g=0.025
a0=2.0
om0=(k**2-g**2)**0.5
x = np.array([2.0,-0.8]) # initial position as before
time = np.linspace(0.0,300.0,100) # evaluation times
""" calling the solver """
solution = odeint(diffeq, x, time, args=(k,g,a0,om0))
"""plotting the result"""
plt.plot(time,solution[:,0])
Out[14]:
[<matplotlib.lines.Line2D at 0x7ca2f148f610>]
In [15]:
# quantum harmonic oscillator
def pot(z):
return z*z
# Schroedinger equation
def diffeq(psi, z, lam):
#define/clear derivative
deriv= np.zeros(2)
deriv[0] = psi[1]
deriv[1] = (pot(z) - lam) * psi[0]
return deriv
z = np.linspace(0,4,100) # evaluation times
""" calling the solver
boundary conditions at x=0 """
psi = np.array([1,0]) # even function: amplitude arbitrarily =1, slope = 0
solution1 = odeint(diffeq, psi, z, (1,))
psi = np.array([0,1]) # odd function: amplitude =0, slope arbitrarily=1
solution2 = odeint(diffeq, psi, z, (3,))
psi = np.array([1,0]) # even function
solution3 = odeint(diffeq, psi, z, (5,))
In [16]:
"""plotting the result"""
plt.plot(z,solution1[:,0],z,solution2[:,0],z,solution3[:,0])
plt.legend(['ground state','first state','second state'])
plt.title('psi')
plt.show()
plt.plot(z,solution1[:,0]**2,z,solution2[:,0]**2,z,solution3[:,0]**2)
plt.legend(['ground state','first state','second state'])
plt.title('|psi* psi|')
Out[16]:
Text(0.5, 1.0, '|psi* psi|')
Radioactive Decay¶
In [17]:
t,counts,err=np.loadtxt('data/new_decay_data.txt',unpack=True)
In [18]:
def chisquare(data,sigma,model):
return ((data-model)**2/sigma**2).sum()
In [19]:
plt.errorbar(t,counts,yerr=err, marker='o', markersize=3, capsize=3, linestyle='none')
plt.xlabel('time [d]')
plt.ylabel('counts / sec')
Out[19]:
Text(0, 0.5, 'counts / sec')
Define a function that calculates the probability to obtain a greater value for chisquare:
In [20]:
""" arguments are: chisquare, number of constraints, number of parameters """
def prob(xsq,n,p):
return 1-chi2.cdf(xsq,n-p)
In [21]:
""" exponential decay """
def decay(x,a,tau):
return a*np.exp(-x/tau)
In [22]:
popt, pcov = curve_fit(decay,t,counts,sigma=err, absolute_sigma=True)
xsq=chisquare(counts,err,decay(t,*popt))
print('reduced chisquare={0:5.2f} probability={1:8.6f}\n'.format(
xsq/(len(counts)-len(popt)), prob(xsq,len(t),len(popt))))
reduced chisquare= 3.78 probability=0.000000
reduced chisquare is pretty high! probability < 1e-6
In [23]:
plt.errorbar(t,counts,yerr=err, marker='o', markersize=3, capsize=3, linestyle='none')
plt.plot(t,decay(t,*popt),lw=3,color='red')
plt.xlabel('time [d]')
plt.ylabel('counts / sec')
Out[23]:
Text(0, 0.5, 'counts / sec')
In [24]:
""" now with background """
def decay_w_background(x,a,tau,bg):
return a*np.exp(-x/tau)+bg
In [25]:
popt, pcov = curve_fit(decay_w_background,t,counts,sigma=err)
xsq=chisquare(counts,err,decay_w_background(t,*popt))
print('reduced chisquare={0:5.2f} probability={1:8.6f}\n'.format(
xsq/(len(counts)-len(popt)), prob(xsq,len(t),len(popt))))
reduced chisquare= 0.99 probability=0.512117
reduced chisquare is about one!
In [26]:
plt.errorbar(t,counts,yerr=err, marker='o', markersize=3, capsize=3, linestyle='none')
plt.plot(t,decay_w_background(t,*popt),lw=3,color='red')
plt.xlabel('time [d]')
plt.ylabel('counts / sec')
Out[26]:
Text(0, 0.5, 'counts / sec')
In [ ]: