Curve fitting
1. What does popt and pcov mean?
The value POPT in the curve fit function returns the optimised coeffiencients for the chosen polynomial that is fitted to the data. See the screenshot below.
The curve fit function returns a list of the coefficient values.
PCOV, this is a matrix of variances and covariances. The diagonal of the matrix is the variance of the fit function to the data it has been fit to and the no diagonal values are the covariances of the fitted function to the data. This matrix can give us some insight into how the fitted function behaves with respect to the data.
2. What does np.array(temperature) do?
The numpy array function converts a list of values into an array of the same size and with the same type of values.
3. What does the * in *popt mean?
The star operator allows you to send a list of values into the function using a single container. The * lets the program know that it needs to read into the container for the values inside it and to send those in the order that it reads it.
4. Write code to fit a linear and cubic polynomial for the Cp data. Explain if your results are good or bad.
The linear fit of the data can be seen below.
It is very clear that the line does not correctly approximate the data it was given to fit. A straight line will not be suitable.
A cubic fit to the data can be seen below. The cubic function manages to follow the trend of the data closer that the straight line. It does, however,
fail to accurately approximate the data. The function over shoots and under shoots the data serveral times.
5. What needs to be done in order to make the curve fit perfect?
In order to remedy the inaccuracy of the of the plot, the function used to fit to the data can be changed to one of a higher order. Also the data can be broken up into smaller intervals and the cubic function, the more accurate of the two, can be fit to those smaller intervals. The results of which can be seen below.
The data was broken up into \"boxes\" of intervals. The function has been broken up into intervals of similar gradients. The similarity is decided by a tolerance value that is assumed to be acceptable. The image above on the left is the new fit with the data broken up into intervals, it is seen to closely match the data that is plotted on the right.
The plot above is the data plot superimposed on the functions fitted to the smaller intervals. This proves that the method of using smaller intervals with the cubic function to fit is more appropriate.
The code used to acheive the results above can be seen below
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
#function receives an array of values in t with the function coefficients a,b,c and evaluates the function for each value in t
def poly_Cube(t,a,b,c,d):
return a*pow(t,3) + b*pow(t,2) + c*pow(t,1) + d
def poly_Straight(t,a,b):
return a*pow(t,1) + b*pow(t,0)
#function read the data from a text file and returns two arrays containing the values from each column in the text file.
def read():
temperature = [] # empty array to hold temperature values
cp = [] # empty array to hold specific heat values
#for loop used to open and read file
for line in open(r\'C:\\Users\\Nicole\\Documents\\Chaz Fenner\\Skill Lync\\Python\\Challenges\\Curve fitting\\data\',\'r\'):
values = line.split(\',\') #split the data in each line by a comma
temperature.append(float(values[0])) # append value in loop to array
cp.append(float(values[1])) # append value in loop to array
return(temperature,cp) # reutn to arrays and count
\"\"\"
#this function is used to split a set of values into invervals separated by inflection points
def inflection(temp,cp, count):
grad = [] #empty array to store gradient values
hold = [] #Hold variable
count_b = 1 #initalising a count value
counts = [] # count array to store positions of inflections points
# this for loop is used to determine the gradient of each point
for i in range(0,count):
grad.append((cp[i]- cp[i-1])/(temp[i]- temp[i-1]))
# here the gradients are compared to determine where there is a straight section or a curved section of the curve
for i in range(0,count):
if grad[i] < grad[i-1]: # if the next gradient is less than the previous
hold.append(0)
elif grad[i] > grad[i-1]:
hold.append(1)
while hold[count_b] == hold[count_b - 1]:
count_b = count_b + 1
if hold[count_b] != hold[count_b - 1]:
counts.append(count_b)
#print(count_b)
numb = len(counts)
T = []
TF =[\'NONE\']*(numb+1)
tempA = []
for i in range(numb+1):
if i == 0:
tempA = []
for j in range(counts[i]):
tempA.append(temp[j])
elif i > 0 and numb > 1:
tempA= []
for j in range(counts[i], counts[i + 1]):
tempA.append(temp[j])
elif i == numb:
tempA = []
for j in range(counts[i-1], count):
tempA.append(temp[j-1])
TF[i] = tempA
return(counts,TF,numb,grad,i)
\"\"\"
# this function is used to break up the arrays into intervals of similar changes in gradient
def break_int(temp,cp):
count_b = len(temp) # assing the length of the temp array to count_b
count_ParaT = [] #empty array to store similar intervals of values
count_StrT = [] #empty array to store similar intervals of values
count_ParaCP = [] #empty array to store similar intervals of values
count_StrCP = [] #empty array to store similar intervals of values
grad = [] #array to store
tol = 0.0000001 # tolerance value that will be used to determine similar interval gradients
k = 0 #counter used in for loop
# this loop will calculate the gradient and compare them, the comparisions will be used to separate similar
# dissimilar intervals of gradients
for i in temp:
grad.append((cp[k]- cp[k-1])/(temp[k]- temp[k-1])) # determining the gradient
#finding interval to apply parabola shape
# only initiate code when more than one gradient value has been determined
if k > 0:
if (grad[k] - grad[k-1]) > tol: # if the gradient difference is greater than tolerance then its parabolic
count_ParaT.append(i) # assign temperature value at loop position
count_ParaCP.append(cp[k]) # assign corresponding cp
else : # else the section difference is acceptable and is assumed to be a straight line
count_StrT.append(i) # assign temperature value at loop position
count_StrCP.append(cp[k ]) # assign corresponding cp
k = k + 1
return(count_StrT,count_ParaT,count_StrCP,count_ParaCP) # return parabolic and straight intervals to main function
# here, we split the intervals into separate arrays so that they can be sent separately into
# we simply use a jump in temperature values as an indicator to split
def split_array(count_StrT,count_ParT,count_StrCP,count_ParCP):
count_a = len(count_StrT) # determine the length of the array
temps =[] # create empty array
cps = [] # create empty array
int_countStr = 0 # count value to move to new position in dictionary
temp_dictStr = {\'box_{}\'.format(q): 0 for q in range(int_countStr)} # initialise empty dictionary
cp_dictStr = {\'box_{}\'.format(q): 0 for q in range(int_countStr)} # initialise empty dictionary
# loop through for as many time as length of array
for i in range(1,count_a):
# if the difference between the two positions is one, the interval is similar, save it
if (count_StrT[i] - count_StrT[i-1]) == 1:
temps.append(count_StrT[i]) # add the similar value to the array
cps.append(count_StrCP[i]) # add the corresponding cp value to the array
temp_dictStr[\"box_%d\" % (int_countStr)] = temps # dictionary at position in for loop must equal array
cp_dictStr[\"box_%d\" % (int_countStr)] = cps # dictionary at position in for loop must equal array
else: # else, this is a new interval and we need a new position to save the next inverval
temps = [] # empty the array for the new set of values
cps = [] # empty the array for the new set of values
int_countStr = int_countStr + 1 # move to the next dictionary position
temps.append(count_StrT[i]) # append first value to the array that broke the if statement
cps.append(count_StrCP[i]) # append first value to the array that broke the if statement
temp_dictStr[\"box_%d\" % (int_countStr)] = temps # dictionary at position in for loop must equal array
cp_dictStr[\"box_%d\" % (int_countStr)] = cps # dictionary at position in for loop must equal array
#now for parabola repeat the same process
temps =[] # create empty array
cps = [] # create empty array
count_a = len(count_ParT) # determine the length of the array
int_countPar = 0 # count value to move to new position in dictionary
temp_dictPar = {\'box_{}\'.format(q): 0 for q in range(int_countPar)} # initialise empty dictionary
cp_dictPar = {\'box_{}\'.format(q): 0 for q in range(int_countPar)} # initialise empty dictionary
# loop through for as many time as length of array
for i in range(1,count_a):
# if the difference between the two positions is one, the interval is similar, save it
if (count_ParT[i] - count_ParT[i-1]) == 1:
temps.append(count_ParT[i]) # add the similar value to the array
cps.append(count_ParCP[i]) # add the corresponding cp value to the array
temp_dictPar[\"box_%d\" % (int_countPar)] = temps # dictionary at position in for loop must equal array
cp_dictPar[\"box_%d\" % (int_countPar)] = cps # dictionary at position in for loop must equal array
else:
temps = [] # create empty array
cps = [] # create empty array
int_countPar = int_countPar + 1
temps.append(count_ParT[i]) # add the similar value to the array
cps.append(count_ParCP[i]) # add the corresponding cp value to the array
temp_dictPar[\"box_%d\" % (int_countPar)] = temps # dictionary at position in for loop must equal array
cp_dictPar[\"box_%d\" % (int_countPar)] = cps # dictionary at position in for loop must equal array
return(temp_dictStr,cp_dictStr,temp_dictPar,cp_dictPar,int_countStr) # return parabolic and straight intervals to main function
# call read function to get arrays of values
temp,cp= read()
# send value that were read into the funtion to be split
count_StrT,count_ParaT,count_StrCP,count_ParaCP = break_int(temp,cp)
# send arrays into this function to be splintered into separate arrays in a dictionary
temp_dictStr,cp_dictStr,temp_dictPar,cp_dictPar,int_countStr = split_array(count_StrT,count_ParaT,count_StrCP,count_ParaCP)
#ploting a straight function to the data
plt.figure(1)
var, pcov = curve_fit(poly_Straight,temp,cp)
# fit funcion using coeffiencets determined by curve fit
fit_cp = poly_Straight(np.array(temp), *var)
# plot fitted function
plt.plot(temp,fit_cp, label= \' Cp fit of Straight line\' )
#plot data
plt.plot(temp,cp, label=\'Data\')
#position legend in the best location
plt.legend(loc= \'best\')
#assing a lable to the x axis
plt.xlabel(\'Temperature (K) \')
#assign a lable to the x axis
plt.ylabel(\'Cp (J / kgK) \')
# fit function using cubic function
plt.figure(2)
var, pcov = curve_fit(poly_Cube,temp,cp)
# fit funcion using coeffiencets determined by curve fit
fit_cp = poly_Cube(np.array(temp), *var)
# plot fitted function
plt.plot(temp,fit_cp, label= \' Cp fit of Cubic line\' )
#plot data
plt.plot(temp,cp, label=\'Data\')
#position legend in the best location
plt.legend(loc= \'best\')
#assing a lable to the x axis
plt.xlabel(\'Temperature (K) \')
#assign a lable to the x axis
plt.ylabel(\'Cp (J / kgK) \')
plt.figure(3)
#loop through dictionary
for i in temp_dictStr:
# do curve fitting for intervals in dictionary
var, pcov = curve_fit(poly_Cube,temp_dictStr[i],cp_dictStr[i])
# fit funcion using coeffiencets determined by curve fit
fit_cp = poly_Cube(np.array(temp_dictStr[i]), *var)
# plot fitted function
plt.plot(temp_dictStr[i],fit_cp, label= \'Cp fit in %s\' % (i))
#loop through dictionary
for i in temp_dictPar:
# do curve fitting for intervals in dictionary
var, pcov = curve_fit(poly_Cube,temp_dictPar[i],cp_dictPar[i])
# fit funcion using coeffiencets determined by curve fit
fit_cp = poly_Cube(np.array(temp_dictPar[i]), *var)
# plot fitted function
plt.plot(temp_dictPar[i],fit_cp, label= \'Cp fit in %s\' % (i))
#plot data
plt.plot(temp,cp, label=\'Data\')
#position legend in the best location
plt.legend(loc= \'best\')
#assing a lable to the x axis
plt.xlabel(\'Temperature (K) \')
#assign a lable to the x axis
plt.ylabel(\'Cp (J / kgK) \')
#show plot
plt.show()
The inflection function in the code above should be ignored please.