Specific heat Vs Temperature curve fitting by MATLAB
There are the 4 quantities to measure the goodness of the fit criteria. These are
(i) SSE ( sum of squares due to error)
(ii) R-square
(iii) Adjusted R-square
(iv) Root Mean square error (RMSE)
R-square : `(SSR)/(SST)`
`SSR = sum_(i=1)^N *(f(x(i)-Mean))^2`
SSR = sum of square regression
Where f(x(i)) represent the value of the fitting function at the given point
`SSE = sum_(i=1)^N *(Y(i)-f(x(i)))^2`
SSE = sum of squares due to error
Where Y(i) represents the value of the given data point at i.
SST = SSE + SSR
R-squared measures the proportion of the variation in your dependent variable (Y) explained by your independent variables (X) for a linear regression model.
R-squared ranges from a value of 0 to 1
0 - for the worst possible fit for the data
1 - for the best possible fit for the data
The image taken from "http://www.biostathandbook.com/linearregression.html".
Adjusted R-squared adjusts the statistic based on the number of independent variables in the model.
`R^2 adj = 1 - [((1-R^2)(n-1))/(n-k-1)]`
n = number of data points
K = number of independent variables
There are two possibilites
(i) The newly added independent variables is insignificant
(ii) The newly added independent variables is significant
Take (i) possiblilites in the consideration
`(R^2 adj)darr = 1 - [((1-R^2uarr)(n-1))/(n-kuarr-1)]`
the increase value of k makes slight changes in the value of R-square due to which R-square adjusted value decreases.
Take (ii) possiblilites in the consideration
`(R^2 adj)uarr = 1 - [((1-R^2uarruarr)(n-1))/(n-kuarr-1)]`
the R-square increases greatly with k value increases in this possibilities which makes R-square adjusted value increases.
Root mean square error (RMSE):
`RMSE = sqrt ((sum_(i=1)^N*(Y(i) - f(x(i)))^2)/n)`
RMSE = `sqrt((SSE)/n`
n = total number of data points
The BEST FIT is achieved when the geometric distances from the given points to the fitting curve are minimized, in the least squares sense.
The matlab program for fitting the curve with the data points.
% Preparing the data
CP_data = load('data')
load command use to bring the data from the saved file
300 940.365951202235
301 940.539121845390
302 940.713818755900
303 940.890032724206
304 941.067754561463
305 941.246975099535
306 941.427685191000
307 941.609875709146
308 941.793537547974
309 941.978661622197
310 942.165238867237
311 942.353260239231
312 942.542716715026
313 942.733599292180
314 942.925898988965
315 943.119606844363
316 943.314713918068
317 943.511211290485
318 943.709090062733
319 943.908341356639
320 944.108956314746
321 944.310926100305
322 944.514241897281
323 944.718894910351
324 944.924876364900
325 945.132177507030
326 945.340789603550
327 945.550703941984
328 945.761911830566
329 945.974404598243
330 946.188173594671
temperature = CP_data(:,1); % temperature data
(:,1) all row and 1st colum of the array.
CP = CP_data(:,2); % specific heat data
(:,2) all row of the the 2nd colunm of the array
% Curve fit
% CP = a*T^3+b*T^3+c*T^2+d*T+e
coeff_cp_3 = polyfit(temperature,CP,3)
polyfit is used to fit a a polynomial curve between temperature and CP and the degree of the polynomial is 3.
predict_cp_3 = polyval(coeff_cp_3,temperature);
To obtain the value of the output of the polynomial fit curve "polyval" command is used.
% Curvefit with the linear fit data
% CP = a*T+b
coeff_cp_1 = polyfit(temperature,CP,1);
predict_cp_1 = polyval(coeff_cp_1,temperature);
% Plot the linear curve fit
plot(temperature,CP,'linewidth',2);
plot is used to plot the curve between temperature in xaxis and y-axis having linewidth 2.
hold on
plot(temperature,predict_cp_3,'linewidth',2,'color','g')
plot(temperature,predict_cp_1,'linewidth',2,'color','r')
xlabel('temperature')
ylabel('specific heat at const. pressure')
legend('original CP','cubic CP','linear CP')
legend is used to name the different plot used in the same figure to distinct it.
% calculation of error
Error_1 = CP - predict_cp_1 % Error calculation b/w actual and linear specific heat
sqr_error_1 = Error_1.^2 % squaring the error value
sum_error_1 = sum(sqr_error_1); % sum of the square error value
n = length(CP) % total no of data point of CP
RMSE_1 = (sum_error_1/n).^0.5 % calculating the root mean square error value
Error_3 = CP - predict_cp_3 % Error calculation b/w actual and cubic specific heat
sqr_error_3 = Error_3.^2 % squaring the error value
sum_error_3 = sum(sqr_error_3); % sum of the square error value
n = length(CP) % total no of data point of CP
RMSE_3 = (sum_error_3/n).^0.5 % calculating the root mean square error value
# comparing the value of RMSE betweeen different degree of polynomial.
RMSE_1 = 25.9991
RMSE_3 = 5.4277
% find the R-square of the fit curve
% R^2 = SSR/SST,
% R^2 ranges from 0 to 1
% R^2 0 values for the worst possible fit for the data
% R^2 1 values for the best possible fit for the data
%R^2 calculation for the linear fit curve
mean_cp = sum(CP)/n; % calculating the mean vaule of the specific heat of the data
error_cp_1 = predict_cp_1 - mean_cp;
sqr_error_cp_1 = (error_cp_1).^2;
SSR_sum_error_cp_1 = sum(sqr_error_cp_1);
SST_error_cp_1 = (SSR_sum_error_cp_1) + sum_error_1
R_square_cp_1 = (SSR_sum_error_cp_1)/(SST_error_cp_1);
%R^2 calculation for the cubic fit curve
mean_cp = sum(CP)/n; % calculating the mean vaule of the specific heat of the data
error_cp_3 = predict_cp_3 - mean_cp;
sqr_error_cp_3 = (error_cp_1).^2;
SSR_sum_error_cp_3 = sum(sqr_error_cp_1);
SST_error_cp_3 = (SSR_sum_error_cp_1) + sum_error_3
R_square_cp_3 = (SSR_sum_error_cp_3)/(SST_error_cp_3);
# The value of R-square to characterisze the fiteness of the curve fit .
R_square_cp_1 = 0.9249
R_square_cp_3 = 0.9965
the R-square value for third degree polynomial is close to 1 so it may be used for the best fit curve for the above data point.
clear all
close all
clc
% Preparing the data
CP_data = load('data')
temperature = CP_data(:,1); % temperature data
CP = CP_data(:,2); % specific heat data
% Curve fit
% CP = a*T^3+b*T^3+c*T^2+d*T+e
coeff_cp_3 = polyfit(temperature,CP,3)
predict_cp_3 = polyval(coeff_cp_3,temperature);
% Curvefit with the linear fit data
% CP = a*T+b
coeff_cp_1 = polyfit(temperature,CP,1);
predict_cp_1 = polyval(coeff_cp_1,temperature);
% Plot the linear curve fit
plot(temperature,CP,'linewidth',2);
hold on
plot(temperature,predict_cp_3,'linewidth',2,'color','g')
plot(temperature,predict_cp_1,'linewidth',2,'color','r')
xlabel('temperature')
ylabel('specific heat at const. pressure')
legend('original CP','cubic CP','linear CP')
% calculation of error
Error_1 = CP - predict_cp_1 % Error calculation b/w actual and linear specific heat
sqr_error_1 = Error_1.^2 % squaring the error value
sum_error_1 = sum(sqr_error_1); % sum of the square error value
n = length(CP) % total no of data point of CP
RMSE_1 = (sum_error_1/n).^0.5 % calculating the root mean square error value
Error_3 = CP - predict_cp_3 % Error calculation b/w actual and cubic specific heat
sqr_error_3 = Error_3.^2 % squaring the error value
sum_error_3 = sum(sqr_error_3); % sum of the square error value
n = length(CP) % total no of data point of CP
RMSE_3 = (sum_error_3/n).^0.5 % calculating the root mean square error value
% find the R-square of the fit curve
% R^2 = SSR/SST,
% R^2 ranges from 0 to 1
% R^2 0 values for the worst possible fit for the data
% R^2 1 values for the best possible fit for the data
%R^2 calculation for the linear fit curve
mean_cp = sum(CP)/n; % calculating the mean vaule of the specific heat of the data
error_cp_1 = predict_cp_1 - mean_cp;
sqr_error_cp_1 = (error_cp_1).^2;
SSR_sum_error_cp_1 = sum(sqr_error_cp_1);
SST_error_cp_1 = (SSR_sum_error_cp_1) + sum_error_1
R_square_cp_1 = (SSR_sum_error_cp_1)/(SST_error_cp_1);
%R^2 calculation for the cubic fit curve
mean_cp = sum(CP)/n; % calculating the mean vaule of the specific heat of the data
error_cp_3 = predict_cp_3 - mean_cp;
sqr_error_cp_3 = (error_cp_1).^2;
SSR_sum_error_cp_3 = sum(sqr_error_cp_1);
SST_error_cp_3 = (SSR_sum_error_cp_1) + sum_error_3
R_square_cp_3 = (SSR_sum_error_cp_3)/(SST_error_cp_3);
To improve the cubic polynomial fit we have to use the split-fit curve. The data point is plot by splitting in two and more set of data point to obtain the best fit curve.
clear all
close all
clc
% split wise polynomial fitting
cp = load('data');
temperature = cp(:,1);
cp_data = cp(:,2);
temp_1 = cp(1:820,1);
cp_data_1 = cp(1:820,2);
temp_2 = cp(821:1650,1);
cp_data_2 = cp(821:1650,2);
temp_3 = cp(1651:2200,1);
cp_data_3 = cp(1651:2200,2);
temp_4 = cp(2201:3200,1);
cp_data_4 = cp(2201:3200,2);
% split curve_fit
coeff_split_1 = polyfit(temp_1,cp_data_1,3);
pred_cp_1 = polyval(coeff_split_1,temp_1);
coeff_split_2 = polyfit(temp_2,cp_data_2,3);
pred_cp_2 = polyval(coeff_split_2,temp_2);
coeff_split_3 = polyfit(temp_3,cp_data_3,3);
pred_cp_3 = polyval(coeff_split_3,temp_3);
coeff_split_4 = polyfit(temp_4,cp_data_4,3);
pred_cp_4 = polyval(coeff_split_4,temp_4);
plot(temperature,cp_data,'linewidth',4,'color','b')
hold on
plot(temp_1,cp_data_1,'linewidth',2,'color','r')
plot(temp_2,cp_data_2,'linewidth',2,'color','r')
plot(temp_3,cp_data_3,'linewidth',2,'color','r')
plot(temp_4,cp_data_4,'linewidth',2,'color','r')
xlabel('temperature[K]');
ylabel('specific heat at const. pressure [KJ/kg-mol]');
legend('orignal curve','split fit curve')polyfit Fit polynomial to data.
P = polyfit(X,Y,N) finds the coefficients of a polynomial P(X) of degree N that fits the data Y best in a least-squares sense. P is a row vector of length N+1 containing the polynomial coefficients in
descending powers, P(1)*X^N + P(2)*X^(N-1) +...+ P(N)*X + P(N+1).
[P,S] = polyfit(X,Y,N) returns the polynomial coefficients P and a structure S for use with POLYVAL to obtain error estimates for predictions. S contains fields for the triangular factor (R) from a QR
decomposition of the Vandermonde matrix of X, the degrees of freedom (df), and the norm of the residuals (normr). If the data Y are random, an estimate of the covariance matrix of P is (Rinv*Rinv')*normr^2/df, where Rinv is the inverse of R.
[P,S,MU] = polyfit(X,Y,N) finds the coefficients of a polynomial in
XHAT = (X-MU(1))/MU(2) where MU(1) = MEAN(X) and MU(2) = STD(X). This centering and scaling transformation improves the numerical properties of both the polynomial and the fitting algorithm.
clear all
close all
clc
% split wise polynomial fitting
cp = load('data');
temperature = cp(:,1);
cp_data = cp(:,2);
temp_1 = cp(1:820,1);
cp_data_1 = cp(1:820,2);
temp_2 = cp(821:1650,1);
cp_data_2 = cp(821:1650,2);
temp_3 = cp(1651:2200,1);
cp_data_3 = cp(1651:2200,2);
temp_4 = cp(2201:3200,1);
cp_data_4 = cp(2201:3200,2);
% split curve_fit
[coeff_split_1,s,mu] = polyfit(temp_1,cp_data_1,3);
pred_cp_1 = polyval(coeff_split_1,temp_1,s,mu);
[coeff_split_2,s,mu] = polyfit(temp_2,cp_data_2,3);
pred_cp_2 = polyval(coeff_split_2,temp_2,s,mu);
[coeff_split_3,s,mu] = polyfit(temp_3,cp_data_3,3);
pred_cp_3 = polyval(coeff_split_3,temp_3,s,mu);
[coeff_split_4,s,mu] = polyfit(temp_4,cp_data_4,3);
pred_cp_4 = polyval(coeff_split_4,temp_4,s,mu);
plot(temperature,cp_data,'linewidth',4,'color','b')
hold on
plot(temp_1,cp_data_1,'linewidth',2,'color','r')
plot(temp_2,cp_data_2,'linewidth',2,'color','r')
plot(temp_3,cp_data_3,'linewidth',2,'color','r')
plot(temp_4,cp_data_4,'linewidth',2,'color','r')
xlabel('temperature[K]');
ylabel('specific heat at const. pressure [KJ/kg-mol]');
legend('orignal curve','split fit curve')