Exercise solutions Day 04 (April 11th)¶
Python course summer semester 2024 University of Heidelberg
12.1, 12.2
In [15]:
def date2str(li):
[year,month,day,hour,minute,second] = li
date="{0:04d}-{1:02d}-{2:02d}".format(year,month,day)
time="{0:02d}:{1:02d}:{2:02d}".format(hour,minute,second)
s="{0} {1}".format(date,time)
return s
def str2date(s):
date, time =s.split()
year,month,day = date.split('-')
hour,minute,second = time.split(':')
return [int(year),int(month),int(day),int(hour),int(minute),int(second)]
# foward and backward
print (str2date(date2str([2017,4,5,17,43,12])))
print (date2str(str2date('2017-04-05 17:43:12')))
In [16]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.stats import chi2
14.01
In [17]:
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 [18]:
def poly1 (x,a0,a1):
return a0*x+a1
def poly2 (x,a0,a1,a2):
return a0*x**2+a1*x+a2
In [19]:
xp=np.linspace(0,10,100)
plt.errorbar(x,y,yerr=e,fmt="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[19]:
14.02
In [20]:
def temp_seasonal(x,a0,a1,a2):
return a0*np.cos(2*np.pi*x+a1)+a2
In [21]:
# 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.plot(date,temperature)
plt.plot(date,temp_seasonal(date,*popt),color='red',lw=3)
plt.xlim(2008,2012)
Out[21]:
In [22]:
print ("temp=",popt[0],"*cos(2*pi*x+",popt[1],")+",popt[2],")")
The average temperature is 9.0 degrees, the minimum/maximum are -3.3/21.4 degrees.
The meaning of popt[1]:
In [23]:
diff=(4*np.pi-popt[1])
print (365*diff/(2*np.pi))
This means that according to this model the coldest time of year is about 15 days after New Year.
15.01
In [24]:
from scipy.integrate import quad
from scipy.integrate import simps
def gaussian(x,mu,sig):
return (1/(np.sqrt(2*np.pi)*sig))*np.exp(-(x-mu)**2/(2*sig**2))
In [25]:
quad(gaussian,-10,10,args=(0,1))
Out[25]:
In [26]:
x=np.linspace(-10,10,100)
y=gaussian(x,0,1)
simps(y,x=x)
Out[26]:
15.02
In [27]:
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [28]:
#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[28]:
In [29]:
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[29]:
In [30]:
# 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 [31]:
"""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[31]:
Radioactive Decay¶
In [32]:
t,counts,err=np.loadtxt('data/new_decay_data.txt',unpack=True)
In [33]:
def chisquare(data,sigma,model):
return ((data-model)**2/sigma**2).sum()
In [34]:
plt.errorbar(t,counts,yerr=err)
plt.xlabel('time [d]')
plt.ylabel('counts / sec')
Out[34]:
Define a function that calculates the probability to obtain a greater value for chisquare:
In [35]:
""" arguments are: chisquare, number of constraints, number of parameters """
def prob(xsq,n,p):
return 1-chi2.cdf(xsq,n-p)
In [36]:
""" exponential decay """
def decay(x,a,tau):
return a*np.exp(-x/tau)
In [37]:
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 is pretty high! probability < 1e-6
In [38]:
plt.errorbar(t,counts,yerr=err)
plt.plot(t,decay(t,*popt),lw=3,color='red')
plt.xlabel('time [d]')
plt.ylabel('counts / sec')
Out[38]:
In [39]:
""" now with background """
def decay_w_background(x,a,tau,bg):
return a*np.exp(-x/tau)+bg
In [40]:
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 is about one!
In [41]:
plt.errorbar(t,counts,yerr=err)
plt.plot(t,decay_w_background(t,*popt),lw=3,color='red')
plt.xlabel('time [d]')
plt.ylabel('counts / sec')
Out[41]:
