Curve Fitting Using Python

1. AIM:

• To perform a curve fitting on a given data set using Python software.

2. GOVERNING EQUATION:

1. For quadratic equations, `C_p=aT^2+bT+c` and so on for other degrees of equations.
2. For `R^2` method, `R^2 = \frac{SSR}{SST}where \sum(f(x(i))-mean)^2 and SST = SSR + \sum(Y(i) - f(x(i)))^2`
3. For RMSE method, RMSE=`\sqrt{\frac{SSE}{n}}`
   • where a, b, c are coefficients of some degree of
     • `f(x(i))` is the y coordinate of the fitted curve having equation y = f(x)
     • Y(i) is the y coordinate of the data set ( `C_p`) here in this case.

3. OBJECTIVES OF THE PROJECT:

1. Write Python code to fit a polynomial through the cp data
2. Plotting the data and checking error using different methods
3. To find what does popt and pcov mean?
4. To find what does np.array(temperature) do?
5. To find what does the * in *popt mean?
6. What needs to be done in order to make the curve fit perfect?

4. FLOW OF THE PROJECT:

• Collect sample data of which you want to get a curve.
• Differentiate the data and assign some variables.
• Use Python to fit a polynomial degree curve through those variables and thus obtain coefficients of the variable Temperature.
• Plot the original dataset as well as curve that is being fitted through those datasets.
• Find error so as to get an approximation of how good your curve fits the data sets.

5. PYTHON CODE:

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

#Curve fit function quadratic  
def func(t, a, b, c):
	return a*pow(t,2) + b*t + c

#Cubic function
def cubic(t,a,b,c,d):
	return a*pow(t,3) + b*pow(t,2) + c*t + d

#4th order polynomial
def fourth(t,a,b,c,d,e):
	return a*pow(t,4) + b*pow(t,3) + c*pow(t,2) + d*t + e

#8th order polynomial
def eight(t,a,b,c,d,e,f,g,h,i):
	return a*pow(t,8)+b*pow(t,7)+c*pow(t,6)+d*pow(t,5)+e*pow(t,4)+f*pow(t,3)+g*pow(t,2)+h*t+i

#Reading thermodynamic data file
def read_file():
	temperature = []
	cp = []
	for line in open('data','r'):             #for every line in the file data
		values = line.split(',')              #split the line in parts having a ',' in between
		temperature.append(float(values[0]))  #add the first column in temperature empty array
		cp.append(float(values[1]))           #add the second column in cp empty array
	return [temperature, cp]

def error(fit_cp1,cp):                        #R and rmse calculator function
	#checking error by R-square method
	#finding the SSR
	ssr=0
	sse=0
	meancp = sum(cp)/len(cp)                  #calculating mean of cp dataset using sum and len functions
	for i in range(0,len(cp)):
		b = pow((fit_cp1[i]-meancp),2);       #looping so as to generate values of b for every value of fit_cp in dataset
		ssr=ssr+b                             #adding all the values of b

	for j in range(0,len(cp)):
		c = pow((cp[j]-fit_cp1[j]),2)         #looping so as to generate value of c for every value of cp and fit_cp
		sse = sse+c                           #adding all values of c

	sst = ssr+sse
	R = ssr/sst                               #final value of R-square

	#finding RMSE
	rmse = pow((sse/len(cp)),0.5)             #final value of RMSE]
	print('R^2 =',R,'  and  RMSE =',rmse,'\n')


#Main code
temperature, cp = read_file()                    #call the function 
popt, pcov = curve_fit(func, temperature, cp)    #obtain values of polynomial coefficients
fit_cp1 = func(np.array(temperature),*popt)      #find the value of the polynomial
z1=error(fit_cp1,cp)                             #calling the function to calculate R^2 and rmse values

popt, pcov = curve_fit(cubic, temperature, cp)
fit_cp2 = cubic(np.array(temperature),*popt)
z2 = error(fit_cp2,cp)

popt, pcov = curve_fit(fourth, temperature, cp)
fit_cp3 = fourth(np.array(temperature),*popt)
z3 = error(fit_cp3,cp)

popt, pcov = curve_fit(eight, temperature, cp)
fit_cp4 = eight(np.array(temperature),*popt)
z4 = error(fit_cp4,cp)


plt.figure()                                                  #opening a new plot window
plt.plot(temperature, cp, color='blue',linewidth=3)           #plotting original data set
plt.plot(temperature, fit_cp1, color='red',linewidth=3)       #Plotting the obtained curve on the same graph
plt.legend(['Actual data','Curve fit'])
plt.xlabel('Temperature (K)')
plt.ylabel('Cp')
plt.savefig('plot 1')
plt.show()

plt.figure()
plt.plot(temperature, cp, color='blue',linewidth=3)
plt.plot(temperature, fit_cp2, color='yellow',linewidth=3)
plt.legend(['Actual data','Curve fit'])
plt.xlabel('Temperature (K)')
plt.ylabel('Cp')
plt.savefig('plot 2')
plt.show()

plt.figure()
plt.plot(temperature, cp, color='blue',linewidth=3)
plt.plot(temperature, fit_cp3, color='orange',linewidth=3)
plt.legend(['Actual data','Curve fit'])
plt.xlabel('Temperature (K)')
plt.ylabel('Cp')
plt.savefig('plot 3')
plt.show()

plt.figure()
plt.plot(temperature, cp, color='blue',linewidth=3)
plt.plot(temperature, fit_cp4, color='black',linewidth=3)
plt.legend(['Actual data','Curve fit'])
plt.xlabel('Temperature (K)')
plt.ylabel('Cp')
plt.savefig('plot 4')
plt.show()

Code explaination:

1. First we have loaded the data given in the challenge itlself, by using open() function used under user created function read_file().
2. Then we have assigned the Cp and temperature values to variables using .append() function because the data set is a CSV i.e Comma Seperated Values file having 2 columns, i.e. seperated by a comma. Thus line splitting is done on that ','
3. Then use curve_fit function whose inputs will be the polynomial function, temperature and CP
4. From the output of curve_fit() function, we obtain the coefficients a, b and c.
5. Next function error() is used for checking the maximum error of the fitting of the curve. 2 methods are used here, namely 1. R-square method and 2. Root mean square method.
6. The value of R in R-square method should lie between 0 and 1 where 0 means poor fit and 1 means . Typically a value greater than 0.95 means a good fit.
7. The RMSE value in the next method should be as less as possible to be called as a good fit.
8. Later on, using plot function, plot the graph of Cp of original dataset and fit_cp1 got from puting the values obtained from the curve_fit into the polynomial. Compare the fit of the curve visually.
9. The plotting and error calculation for different degree of polynomial is done by calling the function in the main code.

6. OBSERVATION:

For quadratic and cubic polynomials,

As the degree is increasing, the curve has more flexiblity and thus can take the shape of the initial dataset in a much better fasion than that of lower polynomial

For 4th and 8th order polynomial,

• Command window:

In increasing order of degree of polynomial: (2,3,4,8)

7. ANSWER CONCLUSIONS:

1. What does popt and pcov mean?

• popt is the array of coefficients of polynomials taken into consideration.
• pcov is a 2d array containing the estimated covariance of popt.

2. What does np.array(temperature) do?

• Helps convert float type array into an integer type array. If this command is not used, it will throw an error saying that it cant use float values

3. What does the * in *popt mean?

• means that the function is referring to the address of these values in order to retrieve them and use them in some other function.

5. As per command window and visual inspection, we can see that as the order of polynomial increases, `R^2` value starts approaching 1 and RMSE value also decreases. This is a good sign. A good fit is generally that one whose `R^2` value is above 0.95. Thus all fittings are good fits as above.

So, in order to make your curve a perfect fit, go on increasing the degree of the polynomial. Also, if at a certain area, even after increasing the order does not help, then you can break the dataset in multiple parts and use higher order polynomials only where there is a complex shape of the curve and lower order polynomials where there is a simpler shape of polynomial

Also, changing the function can also help achieve a perfect fit.

A perfect is the fit where the curve of values of cp found via curve fitting perfectly matches on top of the initial dataset curve.