Aim: Write a python program to perform curve fitting and find goodness of a function.
Objective:
1)To write codes to fit a linear and cubic polynomial for the Cp data
2)To plot the linear and cubic fit curves along with the raw data points.
3)To measure the fitness characteristics for both the curves.
Procedure:
1)First, we import inbuilt libraries necessary for the code such as math, numpy, mathplotlib.pyplot and from scipy.optimize we import curve_fit.
2)We define a function for a linear as well as cubic curve fit then with the return command we get the linear and cubic equation.
3)Then we load the data and two empty arrays are created for both temperature and cp.
4)For loop is used to read each line and with the help of open command we open the file and r is used to read it.
5)The temperature and cp in the data file are separated by a (‘,’) operator.
6)Unlike the text file we have to float the data in the input file and append them in to the empty array’s. Which returns the values of temperature and cp as a column matrix.
7)Then we come to the main program where 'popt' and 'pcov' are used as outputs for storing parameters obtained by 'curve_fit' which takes input as the required function(linear or cubic) , temperature, and Cp.
8)Then 'fit_cp' stores the new values of Cp that are calculated with the help of a temperature array.
9)Plots for both the data i.e predicted as well as the raw data are compared to find the fitness characteristics of a linear and cubic polynomial.
10)To measure the fitness characteristics of both the curves various parameters such as SSE,SST,SSR,R^2,RMSE are calculated.
11)Value that is close to one has been considered as a good fit.
Python code:
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# curve fit function
# linear function
def function_0(t,a,b):
return a*t + b
# cubic function
def function_2(t,a,b,c,d):
return a*pow(t,3)+b*pow(t,2)+c*t+d
# 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
temperature , cp = read_file()
popt,pcov=curve_fit(function_0,temperature,cp)
fit_cp=function_0(np.array(temperature),*popt)
plt.plot(temperature,cp,color='red',linewidth=3)
plt.plot(temperature,fit_cp,color='blue',linewidth=3)
plt.legend(['Actual Data','curve fit'])
plt.xlabel('Temperature[K]')
plt.ylabel('Cp')
plt.title('Linear Curve Fitting')
plt.show()
# for cubic polynomial
popt,pcov=curve_fit(function_2,temperature,cp)
fit_cp_c=function_2(np.array(temperature),*popt)
plt.plot(temperature,cp,color='red',linewidth=3)
plt.plot(temperature,fit_cp,color='blue',linewidth=3)
plt.legend(['Actual Data','curve fit'])
plt.xlabel('Temperature[K]')
plt.ylabel('Cp')
plt.title('Cubic Curve Fitting')
plt.show()
# Measuring the fitness characetristics of the curves
# Methods to calculate goodness of fit
# SSE (sum of error squared)
# SSR (sum of squares of the regression)
# SST = SSE + SSR
# R^2 = SSR/SST
# RMSE = (SSE/l)^0.5
# Finding mean of Cp data
# mean = ( sum of all elements)/(total number of elements)
s = np.sum(cp)
l = np.size(cp)
m = s/l
print('Mean of all cp :',m)
print('') # This will print an empty line
# For Linear curve fit
linear_sse = 0
linear_ssr = 0
for i in range(l):
linear_error = abs((np.sum((cp[i] - fit_cp[i]))))
linear_sse = linear_sse+pow(linear_error,2)
linear_ssr = linear_ssr+ np.sum(pow((fit_cp[i] - m),2))
linear_sst = linear_sse + linear_ssr
print('linear_sst :',linear_sst)
linear_R2 = linear_ssr/linear_sst
print('linear_R2 :',linear_R2)
linear_rmse = pow((linear_sse/l),0.5)
print('linear_RMSE :',linear_rmse)
print('')
# For Cubic curve fit
cubic_sse = 0
cubic_ssr = 0
for j in range(l):
cubic_error = abs((np.sum((cp[i] - fit_cp_c[i]))))
cubic_sse = cubic_sse+pow(cubic_error,2)
cubic_ssr = cubic_ssr+ np.sum(pow((fit_cp_c[i] - m),2))
cubic_sst = cubic_sse + cubic_ssr
print('cubic_sst :',cubic_sst)
cubic_R2 = cubic_ssr/cubic_sst
print('cubic_R2 :',cubic_R2)
cubic_rmse = pow((cubic_sse/l),0.5)
print('cubic_RMSE :',cubic_rmse)
#program for split-wise method
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# curve fit function
# linear function
def function_0(t,a,b):
return a*t + b
# cubic function
def function_2(t,a,b,c,d):
return a*pow(t,3)+b*pow(t,2)+c*t+d
# 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
temperature , cp = read_file()
X=0
Y=799
fit_cp=[]
for i in range(0,4):
popt,pcov=curve_fit(function_0,temperature[X:Y+1],cp[X:Y+1])
fit=function_0(np.array(temperature[X:Y+1]),*popt)
X=Y + 1
Y = Y + 800
for z in fit:
fit_cp.append(z)
print(fit_cp)
#plt.plot(temperature,cp,color='red',linewidth=3)
plt.plot(temperature,fit_cp,color='blue',linewidth=3)
plt.legend(['Actual Data','curve fit'])
plt.xlabel('Temperature[K]')
plt.ylabel('Cp')
plt.title('Linear Curve Fitting')
plt.show()
# for cubic polynomial
for i in range(0,4):
popt,pcov=curve_fit(function_2,temperature[X:Y+1],cp[X:Y+1])
fit=function_2(np.array(temperature[X:Y+1]),*popt)
X=Y + 1
Y = Y + 800
for z in fit:
fit_cp.append(z)
plt.plot(temperature,cp,color='red',linewidth=3)
plt.plot(temperature,fit_cp,color='blue',linewidth=3)
plt.legend(['Actual Data','curve fit'])
plt.xlabel('Temperature[K]')
plt.ylabel('Cp')
plt.title('Cubic Splitwise Curve Fitting')
plt.show()
Results:
1)Curve fit for both the equations are plotted and are shown below.
2)Goodness of a fit for linear and cubic curves are as follows
Mean of all cp : 1182.785685525818
linear_sst : 28793673.249584787
linear_R2 : 0.9248776260524184
linear_RMSE : 25.999093068195325
cubic_sst : 62243756.575895965
cubic_R2 : 0.9907133759605076
cubic_RMSE : 13.440079591073248
1. What does popt and pcov mean?
Ans) Popt represents the matrix that stores the optimal values which are nothing but the coefficients of the function.
Pcov is a covariance matrix which stores the values of the coefficients of the function which is nothing but the values of popt in a square matrix.
2. What does np.array(temperature) do?
Ans) This command converts each value of temperature in the dta file in to an array using the NumPy module .
3. What does the * in *popt mean?
Ans)*popt indicates and returns each coefficient that are stored in the popt array.
4. What needs to be done in order to make the curve fit perfect?
Ans)If we increase the order of the polynomial we can make the perfect curve.
Conclusion:
In the end we can conclude that if the value of R^2 is close to 1 we can say that the fit is good. Here the calculated R^2 value for cubic curve fit is 0.99 which is close to 1 compared to 0.92 of a linear fit. With this we can say that with the increase in degree of a polynomial we can achieve a good and a perfect fit.