To fit a linear and a cubic polynomial to the given Cp data and determine the goodness of fit.
Aim:
- To fit a linear and cubic polynomial to the given c_p data
- To plot the linear and cubic fit curves along with linear data points
- To make the curve fit perfectly and to get the best fit
- To improve the cubic fit
Curve fitting
It is the process of constructing a curve or a mathematical function that fits to a series of data points subjected to constraints, as accurately as possible.
Best curve fit
A curve is said to be the best fit when it passes through a scatter plot of data while best expressing the relationships between the points. It is usually determined by finding the root mean square error pf the data points from the curve and then minimizing it. This involves a process known as linear regression.
Perfect curve fit
A curve is said to be the perfect fit when it fits the original data without showing any errors such as centring and scaling errors for that particular degree of polynomial. A perfect fit can always be a best fit but the vice versa is not true.
Program
Given below is the program to fit a linear and a cubic curve on the given data. The load command is used to get the data, which is stored at the same location as the script, and assign it to an array. To generate a curve fit, the coefficients of the equation of the curve are obtained using the polyfit command.
For example, in this case, the x value of the equation is represented by temperature (T) and the y value is represented by `C_p` . For a cubic polynomial, we can write the equation of the curve as:
`C_p` = `p_1*T^3+p_2*T^2+p_3*T+p_4`
Hence, the value of coefficients are predicted using the polyfit command and a predicted value of `C_p` is obtained using polyval command. This predicted `C_p` is then plotted against temperature to get a curve fit.
% program to perform curve fit
clear all
close all
clc
%fetching data
cp_data = load('data.txt');
temperature = cp_data(:,1);
cp = cp_data(:,2);
% Curve fit
%linear polynomial
%cp = a*T + b
coeffs = polyfit(temperature, cp, 1);
predicted_cp1 = polyval(coeffs, temperature);
%compare curve fit with original data
figure(1)
plot(temperature,cp,'linewidth',3,'color','b')
hold on
plot(temperature, predicted_cp1,'linewidth',3,'color','r')
xlabel('temperature [K]')
ylabel('specific heat [KJ/Kmol-K]')
legend('original data','curve fit')
title('Linear curve fit')
% cubic polynomial
%cp = a*T^3 + b*T^2 + c*T + d
coeffs = polyfit(temperature, cp, 3);
predicted_cp2 = polyval(coeffs, temperature);
%compare curve fit with original data
figure(2)
plot(temperature,cp,'linewidth',3,'color','b')
hold on
plot(temperature, predicted_cp2,'linewidth',3,'color','r')
xlabel('temperature [K]')
ylabel('specific heat [KJ/Kmol-K]')
legend('original data','curve fit')
title('Cubic curve fit')
The linear curve fit works fine. However, when we try to perform the cubic curve fit, the curve is generated but a warning, as shown below, is displayed.
In order to rectify the warning, scaling and centring option is used, as described in the code given below.
In the give program, s is the structure obtained using polyfit command whhich is then passed to the polyval command to obtain error estimates. mu is the 2-element vector with centring and scaling values, where the first element, mu(1) represents the mean of x value and the second element, mu(2) represents the standard deviation of the x value. Using these values, polyfit centres x value at 0 and scales it to have unit standard deviation.
Running this code eleminates the warning displayed above.
clear all
close all
clc
%fetching data
cp_data = load('data.txt');
temperature = cp_data(:,1);
cp = cp_data(:,2);
% Improving the cubic curve fit
[CP_coeffs,s,mu] = polyfit(temperature,cp,3);
predicted_CP = polyval(CP_coeffs,temperature,s,mu);
%compare curve fit with original data
figure(1)
plot(temperature,cp,'linewidth',3,'color','b')
hold on
plot(temperature, predicted_CP,'linewidth',3,'color','r')
xlabel('temperature [K]')
ylabel('specific heat [KJ/Kmol-K]')
legend('original data','curve fit')Spliitwise polynomial could be defined in order to obtained a better curve fit. A split-wise polynomial refers to different polynomials for different range of values. Each polynomial would then individually be operated to determine the values of their coefficients according the range for which they are defined. This would allow a better fit as compared to a single polynomial across the entire range of values.
Obtaining a perfect fit
A perfect fit could be obtained using the split-wise polynomial, while considering centring and scaling. The code given below shows the perfect fit of a cubic polynomial using split-wise polynomial at 4 different intervals, considering centring and scaling.
clear all
close all
clc
%fetching data
cp_data = load('data.txt');
temperature = cp_data(:,1);
cp = cp_data(:,2);
%splitwise polynomial for curve fit
%for temperature between 1-800K
temp_1 = cp_data(1:800,1);
cp_1 = cp_data(1:800,2);
%for temperature between 801-1200K
temp_2 = cp_data(801:1600,1);
cp_2 = cp_data(801:1600,2);
%for temperature between 801-2400K
temp_3 = cp_data(1601:2400,1);
cp_3 = cp_data(1601:2400,2);
%for temperature between 2401-3200K
temp_4 = cp_data(2401:3200,1);
cp_4 = cp_data(2401:3200,2);
[coeffs_cp1,s1,mu1] = polyfit(temp_1,cp_1,3);
predicted_cp1 = polyval(coeffs_cp1,temp_1,s1,mu1);
[coeffs_cp2,s2,mu2] = polyfit(temp_2,cp_2,3);
predicted_cp2 = polyval(coeffs_cp2,temp_2,s2,mu2);
[coeffs_cp3,s3,mu3] = polyfit(temp_3,cp_3,3);
predicted_cp3 = polyval(coeffs_cp3,temp_3,s3,mu3);
[coeffs_cp4,s4,mu4] = polyfit(temp_4,cp_4,3);
predicted_cp4 = polyval(coeffs_cp4,temp_4,s4,mu4);
%compare curve fit with original data
figure(2)
plot(temperature,cp,'linewidth',3,'color','b')
hold on
plot(temp_1, cp_1,'linewidth',3,'color','r')
plot(temp_2, cp_2,'linewidth',3,'color','r')
plot(temp_3, cp_3,'linewidth',3,'color','r')
plot(temp_4, cp_4,'linewidth',3,'color','r')
xlabel('temperature [K]')
ylabel('specific heat [KJ/Kmol-K]')
legend('original data','curve fit')
title('splitwise cubic curve fit')
Goodness of fit
The goodness of describes how well the curve fits the given data. Its measurement summarizes the discrepancy between the actual values and the expected values under the model in question. The parameters used to determine the fitness characteristics of the curves are:
1) Sum of squares due to error (SSE): It measures the total deviation of the raw data points to the response values from the fit. A value closer to 0 indicates a better fit.
2) R-squared (`R^2`): It is the square of correlation between the actual and the predicted response values. It is known as Coefficient of determination (R-squared) and it indicates the proportionate amount of variation in the response variable y (`C_p`) explained by the independent variables X (temperature) in the linear regression model.
`R^2` = `(SSR)/( SST) =1− (SSE)/( SST)` where,
SSR = Sum of squared regression
SST = Total sum of squares = SSR + SSE
The larger the R-squared is, the more variability is explained by the linear regression model. It can take on any value between 0 and 1, with a value closer to 1 indicating a better fit. For example, R2 value of 0.9967 means that the fit explains 99.67% of the total variation in the data about the average.
3) Degrees of freedom of adjusted `R^2` : It uses `R^2` and adjusts it based upon the residual degrees of freedom which is defined as (number of response values) - (number of fitted coefficients) estimated from the response values.
4) Adjusted `R^2` : Obtained in (3), it is the best indicator of the fit quality when additional coefficients are added to the model. It only takes the values less than or equal to 1, with value closesr to 1 indicating a better fit.
`(R_(adj)^2) = 1-((n-1)/(n-p))*(SSE)/(SST)` where,
n = number of observations
p = number of regression coefficients. It includes the intercept value, `therefore` p=2 for linear fit.
5) Root Mean Squared Error (RMSE) : It is also known as fit standard error and a value of RMSE closer to 0 indicates a better fit.
% program to perform curve fit
clear all
close all
clc
%fetching data
cp_data = load('data.txt');
temperature = cp_data(:,1);
cp = cp_data(:,2);
%defining goodness of fit (gof)
[cp1,gof1] = fit(temperature,cp,'poly1','Normalize','on') %linear fit
[cp3,gof3] = fit(temperature,cp,'poly3','Normalize','on') %cubic fit
%plotting curves
figure(1)
hold on
plot(cp1,temperature,cp)
plot(cp3,'g')
legend('temp vs cp','linear','cubic')
title('curve fit')
xlabel('temperature')
ylabel('cp')
hold off
%plotting residuals
figure(2)
hold on
plot(cp1,temperature,cp,'residuals')
plot(cp3,temperature,cp,'g','residuals')
hold off
legend('linear','zero line','cubic')
title('residuals')In the above figure, the goodness of fit for both linear and the cubic fit is shown. Based on what we have learned about the parameters of the fitness characterstics and assuming that a fit with the `R^2` > 0.95 can be considered as a good fit, we can say that:
i) Linear fit is not a good fit because `R^2`<0.95, SSE is higher and RMSE is not close to 0.
ii) Cubic fit is a good fit because `R^2`> 0.95, SSE is lesser and RMSE is also close to 0.
Conclusion
It is not always necessary to get the perfect fit because sometimes, it may take a lot of time to get the polynomial for the exact fit while approximate fit could provide satisfactory results. Increasing the number of coefficients in the model, in other words, the degree of the polynomial does not necessarily provide the best fit.
References