Fitness Characteristics of Linear & Cubic polynomial Curve Fitting for Cp data using Python & MATLAB
OBJECTIVE:
To perform and plot the curve fitting for a linear & cubic polynomial for cp data and also explaining the parameters to measure the fitness for both the curves using Python & MATLAB.
THEORY:
Curve fitting is a technique of finding an algebraic relationship that "best" fits a given set of data. It can be mathematically constructed to determine the relationship between observed dataset and its parameter values.
Regression analysis is a way of mathematically sorting out which of those variables does indeed have an impact. It answers the questions: Which factors matter most? Which can we ignore? How do those factors interact with each other? And, perhaps most importantly, how certain are we about all of these factors?
The above pic represents an example of linear Regression line for 50 random points around the line y=1.5x+2.
There are many types of Regression analysis method, here we find the goodness of fit using polynomial regression method. The general formula for function of type y = f(x).
Linear polynomial :
Quadratic polynomial :
Cubic polynomial :
Below figure represents figures in linear,quadratic & cubic fits.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~MATLAB~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
For polynomial curve fit of any order n, the problem is to find its coefficients from the given data of length(n+1) or more. MATLAB provides an easy way through with built in functions polyfit and polyval.
polyfit : Given two vectors x and y, the command a = polyfit(x,y,n) fits a polynomial of order n through the data points(xi,yi) and returns (n+1) coefficients of powers of x in the row vector a.
polyval : Given a data vector x and the coefficients of a polynomial in a row vector a, the command y = polyval(a,x) evaluates the polynomial at the data points xi and generates the values of yi.
Apart finding a curve which satisfies the data points the GOODNESS OF FIT can be evaluated by the stastical parameters given below
- Sum of Squares due to error (SSE)
- R-square
- Adjusted R-square
- Root mean squared error (RMSE)
MAIN PROGRAM
clear all
close all
clc
%To find & load the directory of the file 'data'
%'data' contains temperature & specific heat
cp_data = load('data');
%Assigning both columns
temperature = cp_data(:,1);
cp = cp_data(:,2);
%Plotting
plot(temperature,cp,'linewidth',3,'color','r');
xlabel('Temperature (K)');
ylabel('Specific Heat (KJ/Kmol - K)');
title('Temperature vs Specific Heat');
grid on
hold on
%Linear Curve fit
coefficient_linear = polyfit(temperature,cp,1);
predicted_cp_linear = polyval(coefficient_linear,temperature);
plot(temperature,predicted_cp_linear,'linewidth',2,'color','g')
%Cubic Curve fit
coefficient_cubic = polyfit(temperature,cp,3);
predicted_cp_cubic = polyval(coefficient_cubic,temperature);
plot(temperature,predicted_cp_cubic,'linewidth',2,'color','b')
legend('Original Dataset','Linear','Cubic','location','east')
%To check the goodness of the fit
n = 3200;
%For Linear fit
error_linear = cp - predicted_cp_linear;
SSE_linear = sum(error_linear.^2);
SSR_linear = sum((predicted_cp_linear - mean(cp)).^2);
SST_linear = SSR_linear + SSE_linear;
Rsquare_linear = SSR_linear/SST_linear;
RMSE_linear = sqrt(SSE_linear/n);
%For Cubic fit
error_cubic = cp - predicted_cp_cubic;
SSE_cubic = sum(error_cubic.^2);
SSR_cubic = sum((predicted_cp_cubic - mean(cp)).^2);
SST_cubic = SSR_cubic + SSE_cubic;
Rsquare_cubic = SSR_cubic/SST_cubic;
RMSE_cubic = sqrt(SSE_cubic/n);
Explanation:
- The code is started with standard procedure command such as clear all,close all & clc.
- The file 'data' which contains temperature & specific heat has been loaded into the Matlab program and assigned to variable cp_data.
- The x & y values has been allocated with temperature & cp respectively.
- The linear curve fit has been determined by the commands polyfit and polyval (syntax is explained in theory section).
- Vairable is been assigned for both and is now plotted.
- Now the same method is done to execute cubic curve fit and then plotted.
- Goodness of the fit is now done by finding the error and using the known formulas for SSE, SSR, SST, R^2 & RMSE.
OUTPUT:
Here we can see that the 'cubic' fit closely represents the initial dataset.
TO IMPROVE THE CUBIC POLYNOMIAL FIT
Eventhough the previous program closely represents the original dataset it is now however not perfect, hence we are now using the method of splitting the plot into eight equal portion and then combine them together by using Splitwise method.
PROGRAM for determining Splitwise method
clear all
close all
clc
%PROGRAM DONE BY SPLITWISE METHOD
%To find & load the directory of the file 'data'
%'data' contains temperature & specific heat
cp_data = load('data');
%Assigning both columns
temperature = cp_data(:,1);
cp = cp_data(:,2);
%Splitwise method
temperature_1 = cp_data(1:400,1);
cp_1 = cp_data(1:400,2);
coefficient_1 = polyfit(temperature_1,cp_1,3);
predicted_cp_1 = polyval(coefficient_1,temperature_1);
temperature_2 = cp_data(401:800,1);
cp_2 = cp_data(401:800,2);
coefficient_2 = polyfit(temperature_2,cp_2,3);
predicted_cp_2 =polyval(coefficient_2,temperature_2);
temperature_3 = cp_data(801:1200,1);
cp_3 = cp_data(801:1200,2);
coefficient_3 = polyfit(temperature_3,cp_3,3);
predicted_cp_3 = polyval(coefficient_3,temperature_3);
temperature_4 = cp_data(1201:1600,1);
cp_4 = cp_data(1201:1600,2);
coefficient_4 = polyfit(temperature_4,cp_4,3);
predicted_cp_4 = polyval(coefficient_4,temperature_4);
temperature_5 = cp_data(1601:2000,1);
cp_5 = cp_data(1601:2000,2);
coefficient_5 = polyfit(temperature_5,cp_5,3);
predicted_cp_5 = polyval(coefficient_5,temperature_5);
temperature_6 = cp_data(2001:2400,1);
cp_6 = cp_data(2001:2400,2);
coefficient_6 = polyfit(temperature_6,cp_6,3);
predicted_cp_6 = polyval(coefficient_6,temperature_6);
temperature_7 = cp_data(2401:2800,1);
cp_7 = cp_data(2401:2800,2);
coefficient_7 = polyfit(temperature_7,cp_7,3);
predicted_cp_7 = polyval(coefficient_7,temperature_7);
temperature_8 = cp_data(2801:3200,1);
cp_8 = cp_data(2801:3200,2);
coefficient_8 = polyfit(temperature_8,cp_8,3);
predicted_cp_8 = polyval(coefficient_8,temperature_8);
%Plotting
%original dataset
plot(temperature,cp,'linewidth',3,'color','r');
hold on
%Splitwise method dataset
predicted_cp = [predicted_cp_1; predicted_cp_2; predicted_cp_3; predicted_cp_4; predicted_cp_5; predicted_cp_6; predicted_cp_7; predicted_cp_8]
plot(temperature,predicted_cp,'linewidth',1.5,'color','k')
xlabel('Temperature (K)');
ylabel('Specific Heat (KJ/Kmol - K)');
title('Temperature vs Specific Heat');
legend('Original dataset','Splitwise dataset','location','east')
grid on
n = 3200;
for i = 1:length(cp);
error_splitwise = cp(i) - predicted_cp(i);
SSE_cubic_splitwise = sum(error_splitwise).^2;
SSR_cubic_splitwise = sum(predicted_cp(i) - mean(cp).^2);
end
%Goodness of split
SST_cubic_splitwise = SSR_cubic_splitwise + SSE_cubic_splitwise;
Rsquare_cubic_splitwise = (SSR_cubic_splitwise / SSR_cubic_splitwise);
RMSE_cubic_splitwise = sqrt(SSE_cubic_splitwise / n);
Explanation
- In this program we are executing the same command but with a difference of splitting the cp into 8 equal parts.
- These 8 equal parts have been individually assigned to the variable temperature_1, cp_1, coefficient_1, predicted_cp_1 and so on till ....._8, respectively for all the eight data points.
- Now these 8 points are joined together and plotted against temperature and compared with the original dataset.
- Since it contains splitted data points, for loop has been assigned to find the error, SSE & SSR and SST, R^2 & RMSE has been found with respective to that.
- Goodness of fit has now been calculated.
OUTPUT:
We can see that splitwise dataset coincides exactly with the original dataset.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~PYTHON~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
MAIN PROGRAM
For fitting a Linear polynomial:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
#Deriving functions:
#Linear fit
def func(t,a,b):
return a*t + b
#Thermodynamic data
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 Polynomial
temperature,cp = read_file()
popt,pcov = curve_fit(func,temperature,cp)
fit_cp = func(np.array(temperature),*popt)
#Plotting
plt.plot(temperature,cp,linewidth = 2, color = \"blue\")
plt.plot(temperature,fit_cp,linewidth = 2, color = \"red\")
plt.title(\'Linear Ploynomial Fit\')
plt.legend([\'Original Data\',\'Curve Fit - Linear\'])
plt.xlabel(\'Temperature\')
plt.ylabel(\'Cp\')
plt.grid()
plt.show()Explanation:
- import function is used to import all the required functions such as math, matplotlib.py, numpy & scipy.
numpy is used for assigning same values of same common difference between two values.
matplotlib.py is used for plotting 2d graphs.
scipy is used from which curve_fit is imported.
- Linear polynomial is found by defining the formula given above and returning the values.
- The given data file is now read by for loop command where data represents the filename and r represents the read command.
- Values present inside the data file is splitted by two with a \',\' in between.
- The index O is given to the first value temperature and index 1 to cp.
- In order to use these values in the calculation we have to cast them to the float command.
- Now they are append to the empty arrays of temperature & cp and returned.
- Part of the main program, we know have temperature & cp for the given data file.
- curve_fit function in the scipy module is used for the parameters func, temperature, cp. func is the defined as the formula for linear fit.
- These are now assigned to the output popt, pcov (detailed explanation given in inquest section).
- Now the defined func is passed with temperature array along with coefficient popt passing with a numpy array.(i.e popt is passed instead of a & b in func)
- Values are now been plotted with title, legend, xlabel, ylabel, grid for linear fit and original data.
OUTPUT:
For fitting a Cubic polynomial:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
#Deriving functions:
#Cubical fit
def func(t,a,b,c,d):
return a*pow(t,3) + b*pow(t,2) + c*t + d
#Thermodynamic data
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 Cubic Polynomial
temperature,cp = read_file()
popt,pcov = curve_fit(func,temperature,cp)
fit_cp = func(np.array(temperature),*popt)
#Plotting
plt.plot(temperature,cp,linewidth = 2, color = \"blue\")
plt.plot(temperature,fit_cp,linewidth = 2, color = \"red\")
plt.title(\'Cubical Ploynomial Fit\')
plt.legend([\'Original Data\',\'Curve Fit - Cubical\'])
plt.xlabel(\'Temperature\')
plt.ylabel(\'Cp\')
plt.grid()
plt.show()Explanation:
Procedure is identical as the previous program except for one change where we have to def func for cubic polynomial having four coefficients a,b,c,d.
OUTPUT:
INQUEST:
1)How to make a curve fit perfectly?
Perfect fit means the curve should fit the original curve without showing any errors without any errors. In this case we have used splitwise method which exactly matches with the original data. Apart from this we can also find the perfect fit by increasing the order of the polynomial (i.e) Cubic polynomial degree 3 almost coincides with the original data.
2)How to get the best fit?
In order for a curve to be best fit it should have least sum of square of numbers. There are different methods to reduce that error.
If the R^2 > 0.95, Curve is said to be best fit.
3)What could be done to improve the cubic fit?
Cubic fit can be improved by two methods
- by increasing the polynomials by determining the prediction intervals.
- by splitting them using splitwise method.
4)What does popt and pcov mean?
- popt is an array which returns the optimal values of the co-efficients which is present in the polynomial(linear, quadratic, cubic)
- pcov is an 2d array which returns estimated covariance of the respective plot.
5)What does np.array(temperature) do?
The function of np.array is to input all the values of temperature as an array and not as a single value. This way entire temperature array is given as a argument to the function instead of creating a loop.
6)What does the * in *popt mean?
* in *popt is a pointer operator which indicates the address of the coefficients (a,b,c,d...) which is stored as an array and passed as an argument.
RESULTS:
Hence the plot for curve fitting has been performed for a linear & cubic polynomial for cp data and also explaining the parameters to measure the fitness for both the curves using Python & MATLAB.
REFERENCES:
BOOK : Getting Started with MATLAB - Rudra Pratap
WEBSITE: