12.1, 12.2
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')))
(2017, 4, 5, 17, 43, 12) 2017-04-05 17:43:12
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.stats import chi2
14.01
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)
def poly1 (x,a0,a1):
return a0*x+a1
def poly2 (x,a0,a1,a2):
return a0*x**2+a1*x+a2
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))
[<matplotlib.lines.Line2D at 0x7f7a41376400>]
14.02
def temp_seasonal(x,a0,a1,a2):
return a0*np.cos(2*np.pi*x+a1)+a2
# 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)
(2008, 2012)
print ("temp=",popt[0],"*cos(2*pi*x+",popt[1],")+",popt[2],")")
temp= -9.955183582006475 *cos(2*pi*x+ 12.313405907548399 )+ 9.040845452179305 )
The average temperature is 9.0 degrees, the minimum/maximum are -3.3/21.4 degrees.
The meaning of popt[1]:
diff=(4*np.pi-popt[1])
print (365*diff/(2*np.pi))
14.695112983605249
This means that according to this model the coldest time of year is about 15 days after New Year.
15.01
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))
quad(gaussian,-10,10,args=(0,1))
(1.0, 8.671029987439099e-10)
x=np.linspace(-10,10,100)
y=gaussian(x,0,1)
simps(y,x=x)
1.0
15.02
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
#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])
[<matplotlib.lines.Line2D at 0x7f7a412cba90>]
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])
[<matplotlib.lines.Line2D at 0x7f7a4122f780>]
# 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,))
"""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|')
Text(0.5, 1.0, '|psi* psi|')
t,counts,err=np.loadtxt('data/new_decay_data.txt',unpack=True)
def chisquare(data,sigma,model):
return ((data-model)**2/sigma**2).sum()
plt.errorbar(t,counts,yerr=err)
plt.xlabel('time [d]')
plt.ylabel('counts / sec')
Text(0, 0.5, 'counts / sec')
Define a function that calculates the probability to obtain a greater value for chisquare:
""" arguments are: chisquare, number of constraints, number of parameters """
def prob(xsq,n,p):
return 1-chi2.cdf(xsq,n-p)
""" exponential decay """
def decay(x,a,tau):
return a*np.exp(-x/tau)
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
plt.errorbar(t,counts,yerr=err)
plt.plot(t,decay(t,*popt),lw=3,color='red')
plt.xlabel('time [d]')
plt.ylabel('counts / sec')
Text(0, 0.5, 'counts / sec')
""" now with background """
def decay_w_background(x,a,tau,bg):
return a*np.exp(-x/tau)+bg
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!
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')
Text(0, 0.5, 'counts / sec')
