Curve fitting (MATLAB)

AIM:

To write a program to fit a linear and cubic polynomial for the specific heat data set then calculate the goodness of fit using different parameters and different ways to improve the fit in MATLAB

THEORY:

Curve fitting is the process of constructing a curve or mathematical function that best fits the data points, possibly subjected to constraints. It helps us in understanding the type of relationship between the set of data and helps in predicting the unknown values.

The Goodness of fit of statistical data determines how well it fits a set of observations. Measures of Goodness of fit typically summarize the discrepancy between observed values and the values expected under the model in question. There are four methods to quantify the goodness of fit.

They are:

The Sum of Squares due to Error (SSE) – This statistic measures the total deviation of the response values from the fit to the response values. It is also called the summed square of residuals and is usually labeled as SSE.

$SSE={\sum }_{i=1}^n{w}_i{({}_i-y{}_i)}^2$

R-square – This statistic measures how successful the fit is in explaining the variation of the data. It can also be understood as the square of the correlation between the response values and the predicted response values. It is defined as the ratio of the sum of squares of the regression (SSR) and the total sum of squares (SST).

SSR is defined as  $SSR={\sum }_{i=1}^n{w}_i{(y{}_i-y̅)}^2$

SST is also called the sum of squares about the mean. $SST={\sum }_{i=1}^n{w}_i{({}_i-y̅)}^2$

R – square $R^2=\frac{SSR}{SST}$

Its value range is in between 0 and 1

Adjusted R-square – This statistic uses the R – square statistic defined above, and adjusts it based on the residual degrees of freedom. It is generally the best indicator of the fir quality when you compare two models that are nested – that is, a series of models each of which adds additional coefficients to the previous model.  

Adjusted R-square $=1-\left[\frac{(1-R^2)(n-1)}{n-k-1}\right]$  

where, n – number of data points and k – number of independent variables

Root Mean Squared Error (RMSE) – This statistic is also known as the fit standard error and the standard error of the regression. It is an estimate of the standard deviation of the random component in the data.

RMSE $=\surd \left(\frac{SSE}{n}\right)$   where,  n – total number of data points 

Using these four quantities MATLAB effectively deduces the goodness of fit, typically R – square is greater than 0.95 the fit is assumed to be good.

PROCEDURE:

FOR LINEAR, AND CUBIC FIT

• Start programming in MATLAB with these three commands, as all the previous calculations are erased. “clear” this command clears all the data stored in the workspace previously. “close all” this command closes all the previous figures and plots. “clc” this command clears the command window.
• To start with the curve fit first data is loaded in a variable by using “load” command. Then the data of first column is assigned to a variable “temperature” and the second column contains the specific heat that is assigned to the variable “cp”.
• Now curve fit for 1 degree polynomial is created, “polyfit” command is used to calculate the values of the coefficients and then with “polyval” the values of “predicted_cp1” is calculated with coefficients and the temperature. Now follow the same steps for 3 degree polynomial.
• Now for curve fitting the data of 1 degree polynomial data. Open a new figure window using “figure” command, with “plot” command plot the original data, and then plot the 1 degree polynomial data. Now for labelling the x – axis and y – axis use “xlabel” and “ylabel” respectively. Now provide legends for differentiating the data by using “legends” command, for giving a title using “title” command. Follow the same steps for 3 degree polynomial data.
• Now for calculating SSE, SSR, SST, R_sq, R_sq_adj, RMSE. First for 1 degree polynomial and then for 3 degree polynomial. Calculate the “error”, then mean “M” using “mean” command also provide with the number of independent variables “k”. Now use the formulas for calculating different parameters that helps us determine the Goodness of fit.

FOR IMPROVING THE FIT – The fit is divided into 3 parts so that it can fit properly, if needed it can be divided in more parts to get more precise result

• Start with “clear”, “close all”, “clc” commands to clear the workspace, closing all the figures, and clear the command window.
• To start with the curve fit first data is loaded in a variable by using “load” command. Then the data of first column is assigned to a variable “temperature” and the second column contains the specific heat that is assigned to the variable “cp”.
• Now divide the whole data into 4 parts that is 1 – 800, 801 – 1600, 1601 – 2400, and 2401 – 3200. Now assign the values of temperature (with suffix 1, 2, 3, 4) and cp (with suffix 1, 2, 3, 4) in the variables.
• Now for curve fit use the “polyfit” command for determining the coefficients of the polynomial and “polyval” command calculates the value of cp on the basis of coefficients and temperature. Here I have used a 5th degree polynomial and for that centering and scaling of “polyfit” command is used so that there are no warnings. Output of “polyfit” here is saved in three variables, as the degree of polynomial is therefore it requires four inputs for evaluating the new “cp”
• Now plot the data use the “plot” command, start with the original data, then using “hold on” command plot the data for the ranges 1 – 800, 801 – 1600, 1601 – 2400, and 2401 – 3200. Now for labelling the x – axis and y – axis use “xlabel” and “ylabel” respectively. Now provide legends for differentiating the data by using “legends” command, for giving a title using “title” command.

CODE:

FOR LINEAR, AND CUBIC FIT

% curve fitting
clear 
close all
clc
 
% preparing the data
cp_data = load('data'); %loads the data and then saves it in the form of array
temperature = cp_data(:,1); % stores first column as temperature
cp = cp_data(:,2); % stores second column as cp
 
% curve fit for 1 degree polynomial    
coeffs1 = polyfit(temperature,cp,1); % calculates the value of coefficients of 1 degree polynomial equation
predicted_cp1 = polyval(coeffs1,temperature); % calculates the value of cp based on the value of coefficient and temperature
 
% curve fit for 3 degree polynomial
coeffs3 = polyfit(temperature,cp,3); % calculates the value of coefficients of 3 degree polynomial equation
predicted_cp3 = polyval(coeffs3,temperature);  % calculates the value of cp based on the value of coefficient and temperature
 
% comparing my curve fit for 1 degree polynomial with original data
 
figure(1) % to open a new figure window
plot(temperature,cp,'linewidth',3) % plotting original data
hold on % to plot another set of data in the same graph
plot(temperature,predicted_cp1,'linewidth',3,'color','r')  % plotting the 1 degree polynomial data
xlabel('Temperature [K]') % labelling x-axis
ylabel('Specific Heat [KJ/Kmol-K]')  % labelling y-axis
legend('Original data','Curve fit') % providing legends to differentiate the data
title('1 Degree Polynomial') % giving a title
hold off 
 
% comparing my curve fit for 3 degree polynomial with original data
 
figure(3) % to open a new figure window
plot(temperature,cp,'linewidth',3)   % plotting original data
hold on % to plot another set of data in the same graph
plot(temperature,predicted_cp3,'linewidth',3,'color','r')  % plotting the 3 degree polynomial data
xlabel('Temperature [K]') % labelling x-axis
ylabel('Specific Heat [KJ/Kmol-K]') % labelling y-axis
legend('Original data','Curve fit') % providing legends to differentiate the data
title('3 Degree Polynomial') % giving a title
hold off
 
% calculating the SSE, SSR, SST, R_sq, R_sq_adj, RMSE for 1 dgree polynomial
 
error1 =  predicted_cp1 - cp;
 
n = length(cp) % number of data
k = 2 % number of independent variables
 
%mean
M = mean(cp)
 
% calculating SSE1
sse1 = (error1).^2;
SSE1 = sum(sse1)
 
% calculating SSR1
ssr1 = (predicted_cp1 - M).^2;
SSR1 = sum(ssr1);
 
% calculating SST1
SST1 = SSE1 + SSR1;
 
% calculating R_sq1
R_sq1 = SSR1/SST1
 
% calculating R_sq_adj1
R_sq_adj1 = 1-(((1-R_sq1)*(n-1))/(n-k-1))
 
% calculating RMSE1
RMSE1 = (SSE1/n)^0.5
 
% calculating the SSR, SST, R_sq, R_sq_adj, RMSE for 3 dgree polynomial
 
error3 =  predicted_cp3 - cp;
 
n = length(cp)
k = 3 % number of independent variables
 
%mean
M = mean(cp)
 
% calculating SSE3
sse3 = (error3).^2;
SSE3 = sum(sse3)
 
 
% calculating SSR3
ssr3 = (predicted_cp3 - M).^2;
SSR3 = sum(ssr3);
 
% calculating SST3
SST3 = SSE3 + SSR3;
 
% calculating R_sq3
R_sq3 = SSR3/SST3
 
% calculating R_sq_adj3
R_sq_adj3 = 1-(((1-R_sq3)*(n-1))/(n-k-1))
 
% calculating RMSE3
RMSE3 = (SSE3/n)^0.5

FOR IMPROVING THE FIT

% curve fitting - improving the fit by using centering and scaling
clear 
close all
clc
 
% preparing the data
cp_data = load('data');
temperature = cp_data(:,1); % first column as the temperature
cp = cp_data(:,2); % second column as the cp
 	
temperature1 = temperature(1:800); % temperature is divided in 1 to 800
predicted_cp1 = cp(1:800);  % cp is divided in 1 to 800
 
temperature2 = temperature(801:1600); % temperature is divided in 801 to 1600
predicted_cp2 = cp(801:1600);  % cp is divided in 801 to 1600
 
temperature3 = temperature(1601:2400); % temperature is divided in 1601 to 2400
predicted_cp3 = cp(1601:2400);  % cp is divided in 1601 to 2400
 
temperature4 = temperature(2401:3200); % temperature is divided in 2401 to 3200
predicted_cp4 = cp(2401:3200);  % cp is divided in 2401 to 3200
 
% curve fit for 1 to 800 using centering and scaling
[coeffs1,~,mu] = polyfit(temperature1,predicted_cp1,5); %  calculates the value of coefficients of 5 degree polynomial equation
cp1 = polyval(coeffs1,temperature1,[],mu);  % calculates the value of cp1 based on the value of coefficient1 and temperature1
 
% curve fit for 801 to 1600 using centering and scaling
[coeffs2,~,mu] = polyfit(temperature2,predicted_cp2,5);  % calculates the value of coefficients of 5 degree polynomial equation 
cp2 = polyval(coeffs2,temperature2,[],mu); % calculates the value of cp2 based on the value of coefficient2 and temperature2
 
% curve fit for 1601 to 2400 using centering and scaling
[coeffs3,~,mu] = polyfit(temperature3,predicted_cp3,5);  % calculates the value of coefficients of 5 degree polynomial equation 
cp3 = polyval(coeffs3,temperature3,[],mu); % calculates the value of cp3 based on the value of coefficient3 and temperature3
 
% curve fit for 2401 to 3200 using centering and scaling
[coeffs4,~,mu] = polyfit(temperature4,predicted_cp4,5);  % calculates the value of coefficients of 5 degree polynomial equation 
cp4 = polyval(coeffs4,temperature4,[],mu); % calculates the value of cp4 based on the value of coefficient4 and temperature4
 
 
% plotting 
plot(temperature,cp,'linewidth',3,'o') % original data
 
hold on
plot(temperature1,cp1,'linewidth',3,'color','r') % data from 1 to 800
 
plot(temperature2,cp2,'linewidth',3,'color','g') % data from 801 to 1600
 
plot(temperature3,cp3,'linewidth',3,'color','y') % data from 1601 to 2400
 
plot(temperature4,cp4,'linewidth',3,'color','m') % data from 2401 to 3200
 
axis([0 3600 900 1350]) % setting the axiz limits
hold off
 
xlabel('Temperature [K]') % labelling x-axis
ylabel('Specific Heat [KJ/Kmol-K]') % labelling y-axis
title('Curve Fitting of Specific Heat vs Temperature') % giving a title to the plot
legend('Original Data','Data for 1-800','Data for 801-1600','Data for 1601-2400','Data for 2401-3200') 
grid on

RESULTS:

FOR LINEAR, AND CUBIC FIT

USING THE CURVE FITTING TOOLBOX IN MATLAB TO VALIDATE THE FITNESS CHARACTERISTICS

FOR IMPROVING THE FIT

CONCLUSIONS:

Thus we can see that R – square parameter for the linear polynomial equation is 0.9249 whereas for the cubic polynomial equation it is 0.9967 which means the degree of the polynomial increases the curve fitting also increases. For R – square the higher the number the better the fit, the lower the number the lousier the fit, its value ranges from 0 to 1. By using centering and scaling in “polyfit” command polynomial equation of the 5^th degree can be plotted without getting any warnings.