To write a code for curve fitting and plotting Cp data and to validate it’s fitness characteristics using toolbox on MATLAB
AIM:
To write a MATLAB code to fit a linear, quadratic and cubic polynomial for the Cp data variation with temperature, to plot the data for the same and to validate the fitness characteristics using Curve fitting toolbox in MATLAB.
DESCRIPTION:
What is Curve Fitting?
Curve fitting is the process of constructing a curve, or mathematical functions, which possess the closest proximity to the real series of data. By curve fitting, we can mathematically construct the functional relationship between the observed dataset and parameter values, etc. It is highly effective in the mathematical modeling of some natural processes.
Why is curve fitting used?
Fitted curves can be used as an aid for data visualization, to infer values of a function where no data are available and to summarize the relationships among two or more variables. Extrapolation refers to the use of a fitted curve beyond the range of the observed data and is subject to a degree of uncertainty since it may reflect the method used to construct the curve as much as it reflects the observed data.
And it also helps to do additional mathematical operations which might not be possible in the initial data set.
In engineering and the sciences is often of use to fit a model to measured data, in order to predict behaviors of the modeled system, to understand how measurement errors affect the confidence in model parameters, and to quantify the confidence in the predictive capabilities of the model.
For example, let us take the population data set of a particular country, as there is no formula that can tell you what the population of the country is going to be next year. However, if you do a curve fit you will get a mathematical function on how the population is going to change with time. In other words, you can predict the population in the future time using this. Based on the accuracy of your curve fit the final result can change, but even then, you can get a rough trend, since the curve fit is a mathematical function, you can do additional operations such as derivatives and so on.
In this project I have taken a data set which has 2 columns, the 1st column is the temperature and the 2nd column is the specific heat column
The file link is for the data is –
https://drive.google.com/file/d/1vWi9Vj3AdQf0gqUK1u9jwS8zg92BHGR0/view
How to choose the perfect fit?
For the curve to fit perfectly over the raw data points, the errors(including the centering and scaling (mu)) should also be zero.
To obtain the best fit curve the errors should be minimum but in case if we want the curve to be a perfect fit, the errors should be 0/null.
How to choose the best fit?
Assume data points (x1,y1) , (x2,y2) , (x3,y3) , ... (xn,yn) and assume the best fit for these set of points will be y = f(x), the criterion under which the equation of the curve will be the best fit possible is as follows
Upon finding a curve that satisfies the data points there are 4 quantities that help us measure the goodness of fit criteria or how well the equation is representing the data points. There are,
- The sum of squares due to error (SSE)
- R-square
- Adjusted R-square
- Root mean squared error (RMSE)
1) The Sum of Squares due to Error (SSE):
Suppose each of these n points has an individual error (difference in y co-ordinate values)
The error for a general "i" point can be given by the function.
Error = |y(i) − value of f(x) at x(i) |= | y(i) − f(x(i)) |
To remove the modulus sign, it is effective if we square the error term,
SSE (sum of squared errors)= `∑_(i=1)^n`Error`(i)^2`
assume the proposed equation is an nth degree polynomial of the form,
`f(x)=A1+A2 .x+A3.x^2+…+An.x^(n-1)`
When we differentiate the SSE with each and every coefficient of the proposed equation f(x) we and equate it to zero, we can find the minima and the value of that coefficient.
Upon solving the multiple systems of equations, we can find the values of all the coefficients.
This method is known as the Least Squares Method and is the most commonly used technique as long as the fitting function is a polynomial.
2) R-square:
R-square is the square of the correlation between the response values and the predicted response values. It is also called the square of the multiple correlation coefficient and the coefficient of multiple determination.
Higher the value of R-square the better the fit, the lower it is the lousier the fit, R-square ranges from a value of 0 to 1, 0 for the worst possible fit for the data set and 1 for the best possible fit for the data set.
The value of R-square can be calculated using the following formula,
`R^2=(SSR)/(SST)`
SSR = sum of squares of the regression, which is basically the square of the difference between the value of the fitting function at a given point and the arithmetic mean of the data set.
3) Adjusted R-square (`R_(adj)^2`):
Adjusted R-square is a normalizing variation of R-square which neglects the variation of R-square due to insignificant independent variables.
for, y=f (x1, x2, x3, ..., xn) here, x1, x2, ..., xn are the independent variables and y is the dependent variable.
- If value of x(i) contributes to the output y, then x(i) is a significant dependent variable
- If value of x(i) does not contribute to the output y, then x(i) is an insignificant dependent variable
Adjusted R-square disregards these infinitesimally small increments to R-square, and only considers the variables which actually contribute to the output of the system in a drastic way.
R-square is formulated as shown below,
`R_(adj)^2=1-[((1-R^2)(n-1))/(n-k-1)]`
where,
n - the number of data points.
k - the number of independent variables
4) RMSE (Root Mean Squared Error):
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, and is defined as,
where,
n is the total number of data points available.
Using these four quantities MATLAB effectively deduces the goodness of fit, typically is R-square is greater than 0.95 the fit is assumed to be good.
What could be done to improve the cubic fit?
This centering and scaling transformation (mu) improves the numerical properties of both the polynomial and the fitting algorithm.
Centering and scaling values, specified as a two-element vector. This vector is an optional output from [p,S,mu] = polyfit(x,y,n) that is used to improve the numerical properties of fitting and evaluating the polynomial p.
OBJECTIVES:
- To write a MATLAB code to fit a linear, quadratic and cubic polynomial for the Cp data, to plot the data for the same.
- Explain the parameters used to measure the fitness characteristics for all the curves.
- To validate the fitness characteristics using the Curve fitting toolbox in MATLAB.
CODE:
%%%%%%%%%%%%%%%%% CURVE FITTING CODE %%%%%%%%%%%%%%%%%
clear all
close all
clc
%inputs
cp_data=load('data');
%extracting values
temperature_array=cp_data(:,1); %extracting 1st column of cp_data
cp_array=cp_data(:,2); %extracting 2nd column of cp_data
% for diffent polynomial curves from linear, quadratic and cubic
[coefficents1,S,mu]=polyfit(temperature_array,cp_array,1); %for linear curve, polyfit used to fit n degree polynomial in x and y data points
[coefficents2,S,mu]=polyfit(temperature_array,cp_array,2); %for quadratic curve, polyfit used to fit n degree polynomial in x and y data points
[coefficents3,S,mu]=polyfit(temperature_array,cp_array,3); %for cubic curve, polyfit used to fit n degree polynomial in x and y data points
cp_predicted1= polyval(coefficents1,temperature_array,S,mu); %for linear curve, polyval used to determine the output value of curve when coeffiecnts and the input variable are given
cp_predicted2= polyval(coefficents2,temperature_array,S,mu); %for quadratic curve, polyval used to determine the output value of curve when coeffiecnts and the input variable are given
cp_predicted3= polyval(coefficents3,temperature_array,S,mu); %for cubic curve, polyval used to determine the output value of curve when coeffiecnts and the input variable are given
%%%%%%%%%%%%%%% DETERMINE THE FITNESS CHARACTERISTICS %%%%%%%%%%%%%%%%%%
%%%%%%%%%% Determine SSE (SUM SQUARED ERROR) and R-SQUARE %%%%%%%%%%%%%%%%
%R-SQUARE METHOD - HIGHGHER THE VALUE OF R-SQUARE (NEAR TO 1) BETTER THE FIT
%THE LOWER THE VALUE OF R-SQUARE (NEAR TO 0) THE LOUSIER THE FIT
%R-SQAURE VALUE RANGES FROM 0 TO 1
% 0 IS THE WORST FIT
% 1 IS THE BEST POSSIBLE FIT FOR THE DATA SET
SSE1=0;
SSE2=0;
SSE3=0;
SSR1=0;
SSR2=0;
SSR3=0;
for i=1:length(temperature_array)
SSE1=SSE1+(cp_array(i)-cp_predicted1(i))^2; %SSE for linear curve
SSE2=SSE2+(cp_array(i)-cp_predicted2(i))^2; %SSE for quadratic curve
SSE3=SSE3+(cp_array(i)-cp_predicted3(i))^2; %SSE for cubic curve
mean_cp=mean(cp_array);
SSR1= SSR1+(cp_predicted1(i)-mean_cp)^2; %(SSR)SUM OF SQUARES OF REGRESSION for linear curve
SSR2= SSR2+(cp_predicted2(i)-mean_cp)^2; %(SSR)SUM OF SQUARES OF REGRESSION for quadratic curve
SSR3= SSR3+(cp_predicted3(i)-mean_cp)^2; %(SSR)SUM OF SQUARES OF REGRESSION for cubic curve
end
SSE1
SSE2
SSE3
SSR1
SSR2
SSR3
%calculation of SST
SST1= SSR1+SSE1;
SST2= SSR2+SSE2;
SST3= SSR3+SSE3;
SST1
SST2
SST3
%calculation of R-SQUARE
R_SQAURE1=(SSR1/SST1)
R_SQAURE2=(SSR2/SST2)
R_SQAURE3=(SSR3/SST3)
% calculation of ADJUSTED R-SQUARE
R_SQAURE_ADJ1= 1-((1-R_SQAURE1)*(length(temperature_array)-1))/(length(temperature_array)-1-1)
R_SQAURE_ADJ2= 1-((1-R_SQAURE2)*(length(temperature_array)-1))/(length(temperature_array)-1-1)
R_SQAURE_ADJ3= 1-((1-R_SQAURE3)*(length(temperature_array)-1))/(length(temperature_array)-1-1)
% USING RMSE (ROOT MEAN SQAURE ERROR)- STANDARD ERROR OF REGRESSION
% it is the estimate of standard deviation of the random compoent in data
RMSE1= sqrt(SSE1/length(temperature_array))
RMSE2= sqrt(SSE2/length(temperature_array))
RMSE3= sqrt(SSE3/length(temperature_array))
%%%%%%%%%%%%%%%%%% PLOTTING %%%%%%%%%%%%%%%%%%%%%%%%%
%plotting for linear curve fitting
figure(1); %to have each plot in each individual figure
plot(temperature_array,cp_array,'linewidth',3); %plotting the actual cp curve for linear curve
hold on %to display all figures not just the recent ones
plot(temperature_array,cp_predicted1,'linewidth',3,'color','r'); %plotting the predicted cp curve for linear curve
xlabel('Temperature [k]'); %to provide label in x-axis of plot
ylabel('Specific Heat [KJ/Kmol-K]'); %to provide label in y-axis of plot
title('Linear curve fitting'); %to provide a title to plot
legend('cp actual curve','cp predicted linear curve'); %to distiguish between actual and predicted curve
% num2str used to converts a numeric array into a character array that represents the numbers and strcat concatenates strings horizontally.
string_1= strcat('SSE =',num2str(SSE1),' R-square =',num2str(R_SQAURE1));
string_2= strcat('R-SqaureAdj =',num2str(R_SQAURE_ADJ1),' RMSE =',num2str(RMSE1));
text(800,950,string_1,'color','b','fontsize',9); %is used to write the individual strings
text(800,980,string_2,'color','b','fontsize',9);
%plotting for Quadratic curve fitting
figure(2); %%to have each plot in each individual figure
plot(temperature_array,cp_array,'linewidth',3); %plotting the actual cp curve for quadratic curve
hold on %to display all figures not just the recent ones
plot(temperature_array,cp_predicted2,'linewidth',3,'color','r'); %plotting the predicted cp curve for quadratic curve
xlabel('Temperature [k]'); %to provide label in x-axis of plot
ylabel('Specific Heat [KJ/Kmol-K]'); %to provide label in y-axis of plot
title('Quadratic curve fitting'); %to provide a title to plot
legend('cp actual curve','cp predicted quadratic curve'); %to distiguish between actual and predicted curve
% num2str used to converts a numeric array into a character array that represents the numbers and strcat concatenates strings horizontally.
string_1= strcat('SSE =',num2str(SSE2),' R-square =',num2str(R_SQAURE2));
string_2= strcat('R-SqaureAdj =',num2str(R_SQAURE_ADJ2),' RMSE =',num2str(RMSE2));
text(800,950,string_1,'color','b','fontsize',9); %is used to write the individual strings
text(800,980,string_2,'color','b','fontsize',9);
%plotting for Cubic curve fitting
figure(3); %to have each plot in each individual figure
plot(temperature_array,cp_array,'linewidth',3); %plotting the actual cp curve for cubic curve
hold on %to display all figures not just the recent ones
plot(temperature_array,cp_predicted3,'linewidth',3,'color','r'); %plotting the predicted cp curve for cubic curve
xlabel('Temperature [k]'); %to provide label in x-axis of plot
ylabel('Specific Heat [KJ/Kmol-K]'); %to provide label in y-axis of plot
title('Cubic curve fitting'); %to provide a title to plot
legend('cp actual curve','cp predicted cubic curve'); %to distiguish between actual and predicted curve
% num2str used to converts a numeric array into a character array that represents the numbers and strcat concatenates strings horizontally.
string_1= strcat('SSE =',num2str(SSE3),' R-square =',num2str(R_SQAURE3));
string_2= strcat('R-SqaureAdj =',num2str(R_SQAURE_ADJ3),' RMSE =',num2str(RMSE3));
text(800,950,string_1,'color','b','fontsize',9); %is used to write the individual strings
text(800,980,string_2,'color','b','fontsize',9);
RESULTS AND PLOTS AND VALIDATION OF RESULTS USING CURVE FITTING TOOLBOX:
1.1) LINEAR CURVE FIT RESULT FROM CODE:
1.2) LINEAR CURVE FIT OBTAINED FROM TOOLBOX:
2.1) QUADRATIC CURVE FIT RESULT FROM CODE:
2.2) QUADRATIC CURVE FIT OBTAINED FROM TOOLBOX:
3.1) CUBIC CURVE FIT RESULT FROM CODE:
3.2) CUBIC CURVE FIT OBTAINED FROM TOOLBOX:
EXPLANATION FOR CODE:
- clc, clear all, and close all is used to start the code as it clears the command window, clears the variables from the workspace, and closes all the plot respectively.
- cp_data=load('data'), is used to load the data file directly into MATLAB using load function
- Temperature array and specific heat array is formed by extracting the 1st and 2nd column of the data array loaded in Matlab.
- For determining linear, quadratic and cubic curve we use polyfit and polyval function.
- p = polyfit(x,y,n) finds the coefficients of a polynomial p(x) of degree n that fits the data, p(x(i)) to y(i), in a least-squares sense. The result p is a row vector of length n+1 containing the polynomial coefficients in descending powers.
- y = polyval(p,x) returns the value of the polynomial p evaluated at x . Polynomial p is a vector whose elements are the coefficients of a polynomial in descending powers. x can be a matrix or a vector. In either case, polyval evaluates p at each element of x.
- The fitness characteristics of the data is done by 4 methods using the formulas given in the description, the for loop is used to run the formulas for all the 3 curves for length of temperature array which is 3200 times.
- Plotting is done for all the 3 curve fits using, figure function to have each plot in each figure, hold on is used since without that the MATLAB just plots most recent plot only and not the previous plots, The plot command is used to create a plot between temperature and specific heat for all the processes, Text function is used to write the individual strings which display the values of SSE, R-square, R-square adj, and RMSE values with the first 2 index represents its position in the plot. Xlabel and ylabel functions is used to label the x and y-axis. Title function is used to give title for the plot, Legend function is used to create legends to differentiate between the actual curve and the predicted curve fit, num2str function is used to converts a numeric array into a character array that represents the numbers and strcat concatenates strings horizontally.
- The validation of the curve fits is done using the curve fit toolbox which can be used by going to, apps > curve fitting >enter the fit name and the type of equation required e.g. polynomial equation and provide the degree too, provide the x-data as the temperature_array and the y data by cp_array and hit enter to validate different curves using this.
ERRORS IN THE CODE AND STEPS TO OVERCOME THEM:
This error was fixed using centering and scaling transformation (mu) improves the numerical properties of both the polynomial and the fitting algorithm.
Using, [coefficents3,S,mu]=polyfit(temperature_array,cp_array,3);
Instead of using [coefficents3]=polyfit(temperature_array,cp_array,3);
And using cp_predicted3= polyval(coefficents3,temperature_array,S,mu);
Instead of cp_predicted3= polyval(coefficents3,temperature_array);
CONCLUSION:
1) In this project, I have learned to solve the code to fit a linear, quadratic, and cubic polynomial for the Cp data, to plot the data for the same.
2) From the results we can see that the linear fit has the R-square to be less than 0.95 and the R-square for quadratic and the cubic fit is greater than 0.95, hence quadratic and cubic are good curve fit.
3) From the curve fit toolbox we can validate that the results from the code is similar to the results of the toolbox hence this verifies the code written is correct and since all the results are consistent.
REFERENCES:
https://skill-lync.com/projects/Curve-Fitting-and-criterion-to-choose-best-fit-41860
https://en.wikipedia.org/wiki/Curve_fitting