Python Program: Curve fitting

Aim: Write a Python code to fit a curve for given data set.

Objective:

• Write a code to obtain a linear and cubic polynomial for given data set. (here, Cp vs Temperature values (specific heat at constant pressure))
• Plotting the obtained curves.
• Studying fitness characteristics for the same.

Theory:

Curve fitting:

Several equations of different types can be obtained to express the given data approximately. But the problem is to find the equation of the curve of best fit which may be most suitable for predicting the unknown values. The process of finding such an equation by best fit is known as curve-fitting. There are four methods to do so:

1. Method of Least Squares
2. Method of Moments
3. Method of Group
4. Graphical Methods

Among these Method of Least Square, is probably, the best fit to a unique curve to a given data.

Method of Least Squares:

The method of least squares is a standard approach in regression nalysis to approximate solution of overdetermined systems by minimizing the sum of the residuals made in the results of every single equation.

`->` To fit a straight line `y = a +bx`

1. Substitute the observed set of n values in this equation.
2. Form normal equations for each constant `sum y = na + bsumx; sumxy = asumx + bsumx^2.`
3. Solve these normal equations as simultaneous equations for `a` and `b`.
4. Substitute the values of `a` & `b` in `y = a + bx`, which is the required best fit.

`->` To fit a parabola `y = a + bx + cx^2`

1. Solve the derived normal equations for this polynomial.
2. Substitute the values of `a, b, c.`
3. We have `y = a + bx + cx^2`, the required polynomial of best fit.

Measuring Fitness Characteristics:

After performing curve fitting, we need to check its fitness, i.e. how well our derived function is able to represent the data set. The following are most common methods to check fitness characteristics:

• The sum of squares due to error (SSE)
• R-square method
• Adjusted R-square method
• Root mean square error (RMSE)

1. Sum of squares due to error (SSE)

It is the sum of squares of results or deviations predicted from actual values of data. It is measure of the discrepancy between the data and an estimation model. A small SSE indicates a tight fit of the model to the data.

`SSE = sum_(i=1)^n`Error`(i)^2`

2. R-square is a statistical measure of how close the data are to the fitted regression line. R-square is always between 0 and 1. 0 indicates that the model explains none of the variability of the response data around its mean; 1 indicates that the model explains all the variability of the response data around its mean. In general, the higher the R-square, the better the model fits the real data.

`R^2 = (SSR)/(SST)`

where,

`SSR` = Sum of squares of regression

`SST` = `SSR` + `SST`

3. Adjusted R-square is a normalizing variation of R-square due to insignificant independant variables. Adjusted R-square disregards infinitesimally small increments to R-square and only considers the variable which actually contribute to the output of the system in a drastic way.

`R^2(adj) = 1- [((1-R^2)(n-1))/(n-k-1)]`

where,

`n` = number of data points

`k` = the number of independant

4. Root Mean Square Error (RMSE) is an estimate of the standard deviation of the random component in the data, and is defined as,

`RMSE = sqrt((SSE)/n)`

Program:

#Curve Fitting

#Importing required libraries
from numpy import *
from scipy.optimize import curve_fit
from matplotlib.pyplot import *

#Curve fit function
def linear(t,a,b):
	return(a*t + b)

def quadratic(t,a,b,c):
	return(a*pow(t,2) + b*t + c)

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

def quartic(t,a,b,c,d,e):
	return(a*pow(t,4) + b*pow(t,3) + c*pow(t,2) + d*t + e)

#Reading Thermodynamic Data File
def read_file():
	temperature = []
	cp = []
	for line in open('data','r'):
		values = line.split(',')
		temperature.append(float(values[0]))
		cp.append(float(values[1]))
	return [temperature, cp]

#Main Program for linear fit
temperature, cp = read_file()
popt, pcov = curve_fit(linear,temperature,cp)
fit_cp1 = linear(array(temperature), *popt)

#Main Program for quadratic fit
temperature, cp = read_file()
popt, pcov = curve_fit(quadratic,temperature,cp)
fit_cp2 = quadratic(array(temperature), *popt)

#Main Program for cubic fit
temperature, cp = read_file()
popt, pcov = curve_fit(cubic,temperature,cp)
fit_cp3 = cubic(array(temperature), *popt)

#Main Program for quartic fit
temperature, cp = read_file()
popt, pcov = curve_fit(quartic,temperature,cp)
fit_cp4 = quartic(array(temperature), *popt)

#Plotting
figure(1)
plot(temperature,cp,label="Actual Data")
plot(temperature,fit_cp1,label="Linear Curve Fit")
legend()
title("Linear Curve Fit")
xlabel("Temperature (in K)")
ylabel("Specific Heat (in kJ/kg-K)")
savefig("linear-curve-fit.jpg")
show()

figure(2)
plot(temperature,cp,label="Actual Data")
plot(temperature,fit_cp2,label="Quadratic Curve Fit")
legend()
title("Quadratic Curve Fit")
xlabel("Temperature (in K)")
ylabel("Specific Heat (in kJ/kg-K)")
savefig("quadratic-curve-fit.jpg")
show()

figure(3)
plot(temperature,cp,label="Actual Data")
plot(temperature,fit_cp3,label="Cubic Curve Fit")
legend()
title("Cubic Curve Fit")
xlabel("Temperature (in K)")
ylabel("Specific Heat (in kJ/kg-K)")
savefig("cubic-curve-fit.jpg")
show()

figure(4)
plot(temperature,cp,label="Actual Data")
plot(temperature,fit_cp4,label="Quartic Curve Fit")
legend()
title("Quartic Curve Fit")
xlabel("Temperature (in K)")
ylabel("Specific Heat (in kJ/kg-K)")
savefig("quartic-curve-fit.jpg")
show()

#Calculation of Fitness Characteristics (R^2 method)
mean = sum(cp)/len(cp)

for i in range(len(cp)):
    SSR1 = sum(pow(fit_cp1[i]-mean,2))
    SSE1 = sum(pow(cp[i]-fit_cp1[i],2))
    
    SSR2 = sum(pow(fit_cp2[i]-mean,2))
    SSE2 = sum(pow(cp[i]-fit_cp2[i],2))
    
    SSR3 = sum(pow(fit_cp3[i]-mean,2))
    SSE3 = sum(pow(cp[i]-fit_cp3[i],2))
    
    SSR4 = sum(pow(fit_cp4[i]-mean,2))
    SSE4 = sum(pow(cp[i]-fit_cp4[i],2))
    
print("\n- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -")
print("| R-SQUARE FITNESS CHARACTERISTICS                        |")
print("- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -")
print("| R-square fitness for linear fit    :", SSR1/(SSR1+SSE1),"|")
print("| R-square fitness for quadratic fit :", SSR2/(SSR2+SSE2),"|")
print("| R-square fitness for cubic fit     :", SSR3/(SSR3+SSE3),"|")
print("| R-square fitness for quartic fit   :", SSR4/(SSR4+SSE4),"|")
print("- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -")

Working of Code:

1. The code begins with importing of required libraries. In this program we'll use numpy and matlpotlib.pyplot for calculations and plotting respectively. A new module is introduced, i.e. scipy.optimize and we'll import curve_fit from this module which will generate coeffficients of polynomials.

2. Now, for each type of fit (linear, cubic,....) we'll define a function which would take in the variable (here, t `->` temperature) and the coefficients of respective polynomial. The function would then return the deduced polynomial.

3. We need to extract the file which has the temperature and corresponding cp values. Again, we'll define a function for this purpose named read_file. Here, an empty array of temperature and cp is created. We'll read the file line by line and as the required data is separated by a comma ',' we'll separate the either sides of ',' using split command. The result is stored in an array called values which now has two columns.

4. The first column contains temperature which is appeneded to temperature array using append command. And the second column has cp values which are appended in cp array. The function returns combination of this array.

5. Now, we do the curve fitting. We retrieve temperature and cp values and send it to curve_fit along with the polynomial. This gives us popt and pcov which has coefficients and variance, covariance values. Later on, we pass temperauture and coefficients to the polynomial function and derive cp values, which are stored in fit_cp.

6. Now, we plot the temperature vs cp data (the one read from the file and the curve fitted one). The plotted curves are given suitable lables alog with a title and the plots are saved using savefig command.

7. Although we have fitted the curve but we need to evaluate the fitness parameters. We have used R-square method for this. The R-square for each curve was displayed.

Outputs:

The various fits obtained are as follows:

R-square chracteristics:

Result & Conclusion:

• The R-square fitness was found to increase with increase in degree of polynomial.
• Best fits have higher values of R-square (i.e. difference between actual and predicted values are minimized) whereas perfect fit has R-square value 1 (the traced curve has exact values with zero deviation)
• For complex relations, we can split the data into several parts and evaluate fit for each part to fit the curve.

Questions:

1.What does popt and pcov mean?

`->` The popt arguement stores the best fit coefficients for a given polynomial. The pcov contains the covariance matrix, which indicates the uncertainities and corelations between parameters.

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

`->`It creates a numpy array of temperature so that all the values of temperature are passed sequentially to evaluate its respective cp (by fitting).

3. What does the * in *popt mean?

`->` The * denotes all. The popt arguements stores coefficients of polynomial which needs to be passed into the polynomial's function. And depending on the degree of equation, the number of coeffiients changes. Hence, to ensure that all the required number of coefficients are passed sequentially '*' is provided.

4. Write a code to fit a linear amd cubic polynomial for the Cp data. Explain if your results are good or bad.

`->` The code has been included above. The linear fit barely represents the data accurately whereas cubic fit had some accuracy. Although it was clearly visible from the plots but we can get insights from R-square. The R-square value improved with increase in degree of polynomial.

5. What needs to be done in order to make the curve fit perfect?

`->` For obtaining a perfect curve fit, the data needs to be splitted in many parts and the individual parts needs to be evaluated.

References:

• Python for Mechanical Engineers - Course provided by Skill Lync
• Anaconda - Anaconda Inc.
• JupyterLab - Project Jupyter
• Supporting Material: Higher Engineering Mathematics - BS Grewal
• Supporting Material: Scipy Documentation