5.6. Solutions to exercises#
Exercise 5.1:
# You can also use the option "unpack=True" to have loadtxt send back all columns separately
t2,v2 = np.loadtxt("exercise_data.dat", unpack=True)
print("Loaded", len(v2), "points")
Exercise 5.2:
x = np.arange(1,100.1,0.1)
def function(x):
return 3.05*x+0.95
y = function(x)
data = np.zeros([len(x),2])
data[:,0] = x
data[:,1] = y
np.savetxt('testfile.csv',data,delimiter=';')
Exercise 5.3:
plt.plot(t2,v2)
# ALL PLOTS MUST HAVE LABELS WITH UNITS!!!!
plt.xlabel("$t$ (s)")
plt.ylabel("$v$ (V)")
plt.show()
Exercise 5.4:
# p = 2 looks ok
p = 2
plt.plot(t2**p,v2)
# A decent label(we can use the variable value automatically, and use latex for superscripts)
plt.xlabel("$t^"+ str(p) + "$ (s$^" + str(p) +"$)")
plt.ylabel("v (V)")
plt.show()
Exercise 5.5:
# p = 2 looks ok
plt.plot(np.log(t2),np.log(v2),'o', label="Data")
# Now a straight line
# This is a bit tricky, we need the same x-range as the data
# My best estimate of the slope is 2.3
x = np.log(t2)
p = 2.3
y = p*x
plt.plot(x,y, label="Slope = " + str(p))
plt.legend()
# A decent label(we can use the variable value automatically, and use latex for superscripts)
plt.xlabel("log(t) (log(s))")
plt.ylabel("log(v) (log(V))")
plt.show()
Exercise 5.6:
a=2
p=2.34
plt.plot(t2,v2, 'o', label='Data')
plt.plot(t2,a*t2**p, label='Model')
plt.legend()
plt.xlabel("t (s)")
plt.ylabel("v (V)")
plt.show()
Exercise 5.7:
def f2(t,a,p):
return a*t**p
values, covariance = curve_fit(f2, t2, v2)
a_fit = values[0]
p_fit = values[1]
print(values)
plt.plot(t2,v2, 'o', label='Data')
plt.plot(t2,f2(t2,a_fit,p_fit), label='Model')
plt.legend()
plt.xlabel("t (s)")
plt.ylabel("v (V)")
plt.show()
Exercise 5.8
x_meas = np.arange(0,10.1,0.3)
y_meas = 0.2345*x_meas**2 + 3.045*x_meas + np.random.normal(0,1.5,len(x_meas))
#your code here
plt.figure(figsize=(3,3))
plt.plot(x_meas,y_meas,'k.')
plt.show()
#first guess
def function_fit(x,a):
return a*x**2
val, cov = curve_fit(function_fit,x_meas,y_meas)
x_test = np.linspace(min(x_meas),max(x_meas),1000)
plt.figure(figsize=(3,3))
plt.plot(x_meas,y_meas,'k.')
plt.plot(x_test,function_fit(x_test,*val))
plt.show()
plt.figure(figsize=(3,3))
plt.plot(x_meas,y_meas-function_fit(x_meas,*val),'k.')
plt.show()
#there is a pattern in the residuals, so our formula does not sufffice
#second guess
def function_fit(x,a,b,c):
return a*x**2+b*x+c
val, cov = curve_fit(function_fit,x_meas,y_meas)
x_test = np.linspace(min(x_meas),max(x_meas),1000)
plt.figure(figsize=(3,3))
plt.plot(x_meas,y_meas,'k.')
plt.plot(x_test,function_fit(x_test,*val))
plt.show()
plt.figure(figsize=(3,3))
plt.plot(x_meas,y_meas-function_fit(x_meas,*val),'k.')
plt.show()
plt.figure(figsize=(3,3))
plt.hist(y_meas-function_fit(x_meas,*val),bins=int(len(x_meas)/3))
plt.show()
#looks a bit gauss distributed
print('the formula is: a*x**2+b*x+c with parameters a = %.2f +/- %.2f, b= %.1f+/-%.1f, c=%.1f+/-%.1f' %(val[0],np.sqrt(cov[0,0]),val[1],np.sqrt(cov[1,1]),val[2],np.sqrt(cov[2,2])))
#Note that the number of decimal placed can differ as each time a new dataset is made!
Exercise 5.8*
a_err = np.sqrt(covariance[0,0])
p_err = np.sqrt(covariance[1,1])
# Uncomment to debug :)
#print(a_err)
#print(p_err)
# A bit nerdy but cool: automatically format the string output with the
# correct number of digits :)
# (you can also just read off the number of digits by hand and "hard code" it...)
a_dig = -(int(np.log10(a_err))-1)
p_dig = -(int(np.log10(p_err))-1)
# Uncomment to debug :)
#print(a_dig)
#print(p_dig)
a_fmt = "%." + str(a_dig) + "f"
p_fmt = "%." + str(p_dig) + "f"
# The units on a are a bit funny, so we'll skip the units on this one...
a_msg = "a has the value %s +/- %s" % (a_fmt, a_fmt)
a_msg = a_msg % (a_fit, a_err)
print(a_msg)
p_msg = "p has the value %s +/- %s" % (p_fmt, p_fmt)
p_msg = p_msg % (p_fit, p_err)
print(p_msg)
Exercise 5.9
y = 2.064 * np.cos(0.25*(2*np.pi)*x+0.45)+0.0075*x
#load our data
data = np.loadtxt('dataset.csv',delimiter=';')
x_data = data[:,0]
y_data = data[:,1]
x_test = np.linspace(min(x_data),max(x_data),1000)
#see what it looks like
plt.figure()
plt.plot(x_data,y_data,'k.')
plt.show()
#First guess:
def functie(x,a,b,c):
return a*np.cos(b*x+c)
val, cov = curve_fit(functie,x_data,y_data,p0=(2.05,1.57,0.5))
print(val)
#see what it yields, playing with the inital guesses is required!
plt.figure()
plt.plot(x_data,y_data,'k.')
plt.plot(x_test,functie(x_test,*val))
plt.show()
#looking at the residuals
plt.figure()
plt.plot(x_data,y_data-functie(x_data,*val))
plt.show()
#Our data clearly shows a hidden pattern!!!
#second guess
def functie(x,a,b,c,d):
return a*np.cos(b*x+c)+d*x
val, cov = curve_fit(functie,x_data,y_data,p0=(2.05,1.57,0.5,0.05)) #estimate for d from slope
print(val)
plt.figure()
plt.plot(x_data,y_data,'k.')
plt.plot(x_test,functie(x_test,*val))
plt.show()
#looking at the residuals
plt.figure()
plt.plot(x_data,y_data-functie(x_data,*val))
plt.show()
x = np.linspace(0,6*np.pi,200)
y = 2.064 * np.cos(0.25*(2*np.pi)*x+0.45)+0.0075*x
dataset = np.zeros([len(x),2])
dataset[:,0] = x
dataset[:,1] = y
np.savetxt('dataset.csv',dataset,delimiter=';')
plt.figure()
plt.plot(x,y,'k.')
plt.show()