Analysis on curve fitting of variation of specific heat with change in temperature

Aim of the project:- To analyse the case study of variation of specific heat with increasing temperature for creating a curve fit and evaluating the goodness of the curve fit in matlab.

Theory:- Curve fit is used to derieve the most accurate relation between the given set of data to further use it for processing in multiple applications, to predict the future models that are extension of the original data set etc. Curve fitting, also called regression analysis, is a process of fitting a function to a set of data points. The function can then be used as a mathematical model of the data. Since there are many types of functions (linear, polynomial, power, exponential,etc.), curve fitting can be a complicated process. Many times one has some idea of the type of function that might fit the given data and will need only to determine the coefficients of the function. In other situations, where nothing is known about the data, it is possible to make different types of plots that provide information about possible forms of functions that might fit the data well.

There are two methods: Interpolation and curve fitting,

If data is reliable, it can be plotted by connect the dots. This is piece-wise, linear interpolation. This has limited use as a general function `f(x)`. Since its really a group of small `f(x)`s, connecting one point to the next it doesn’t work very well for data that has built in random error (scatter). A general representation of the interpolation can be seen as,

Curve fitting is capturing the trend in the data by assigning a single function across the entire range. A general graphical difference between interpolation and curve fitting for a simple linear equation `f(x) = ax + b` can be seen as,

Polynomials are functions that have the form: `f(x) = a_nx^n + a_(n-1)x^(n-1) + . . . . +a_1x +a_0`, The co-efficients `a_n,a_(n-1),. . . . ,a_1,a_0` are real numbers and `n` which is a non-negative integer is the degree or order of the polynomial.

In matlab, polynomials are represented by a row vector in which the elements are the coefficients `a_n,a_(n-1),. . . . ,a_1,a_0`. The first element is the coefficient of the x with the highest power. The vector has to include all the coefficients, including the ones that are equal to 0.

The value of a polynomial at a point x can be calculated with the function polyval which has the form:

x can also be a vector or a matrix. In such a case the polynomial is calculated for each element (element-by-element), and the answer is a vector, or a matrix, with the corresponding values of the polynomial.

Polynomials can be used to fit data points in two ways. In one the polynomial passes through all the data points, and in the other the polynomial does not necessarily pass through any of the points, but overall gives a good approximation of the data. The two options are described below.

Polynomials that pass through all the points:
When `n` points `(x_i,y_i)` are given, it is possible to write a polynomial of degree `n-1` that passes through all the points. For example, if two points are given it is possible to write a linear equation in the form `y = mx + b` of that passes through the points. With three points the equation has the form of `y = ax^2+bx+c`. With `n` points the polynomial has the form `a_n,a_(n-1),. . . . ,a_1,a_0`. The coefficients of the polynomial are determined by substituting each point in the polynomial and then solving the `n` equations for the coefficients. Polynomials of high degree might give a large error if they are used to estimate values between data points.

Polynomials that do not necessarily pass through any of the points:
When `n` points are given, it is possible to write a polynomial of degree less than `n-1` that does not necessarily pass through any of the points, but overall approximates the data. The most common method of finding the best fit to data points is the method of least squares. In this method the coefficients of the polynomial are determined by minimizing the sum of the squares of the residuals at all the data points. The residual at each point is defined as the difference between the value of the polynomial and the value of the data.

For example, consider the case of finding the equation of a straight line that best fits four data points as shown in figure below. The points are `(x_1,y_1),(x_2,y_2),(x_3,y_3) and (x_4,y_4)`, and the polynomial of the first degree can be written as `f(x) = a_1x +a_0`. The residual, `R_i`, at each point is the difference between the value of the function at `x_i and y_i, R_i = f(x_i)-y_i`. An equation for the sum of the squares of the residuals `R_i` of all the points is given by

`R = [f(x_1)-y_1]^2 + [f(x_2)-y_2]^2 + [f(x_3)-y_3]^2 + [f(x_4)-y_4]^2`, or, after substituting the equation of the polynomial at each point, by: `R = [a_1x_1+a_0-y_1]^2 + [a_1x_2+a_0-y_2]^2 + [a_1x_3+a_0-y_3]^2 + [a_1x_4+a_0-y_4]^2`.

At this stage R is a function of `a_1 and a_0` . The minimum of `R` can be determined by taking the partial derivative of `R` with respect to `a_1 and a_0` (two equations) and equating them to zero.

`(delR)/(dela_1) = 0 and (delR)/(dela_0) = 0`

This results in a system of two equations with two unknowns `a_1 and a_0`. The solution of these equations gives the values of the coefficients of the polynomial that best fits the data. The same procedure can be followed with more points and higher-order polynomials.

Curve fitting with polynomials is done in MATLAB with the polyfit function, which uses the least squares method. The basic form of the polyfit function is:

For the same set of `m` points, the polyfit function can be used to fit polynomials of any order up to `m-1`. If `n=1` the polynomial is a straight line, if `n=2` the polynomial is a parabola, and so on. The polynomial passes through all the points if `n = m-1`(the order of the polynomial is one less than the number of points). It should be pointed out here that a polynomial that passes through all the points, or polynomials with higher order, do not necessarily give a better fit overall. High-order polynomials can deviate significantly between the data points.

An example for the different types of fit can be seen in the images below considering some random set of seven points: `(0.9,0.9),(3,2.5),(4,5.1),(6,4.5),(8,4.9) and (9.5,6.3)`, The points are fitted using the polyfit function with polynomials of degrees 1 through 6. Each plot in figures shows the same data points, marked with circles, and a curve-fitted line that corresponds to a polynomial of the specified degree.

It can be seen that the polynomial with `n=1` is a straight line, and with `n=2` is a slightly curved line. As the degree of the polynomial increases, the line develops more bends such that it passes closer to more points. When `n=6`, which is one less than the number of points, the line passes through all the points. However, between some of the points, the line deviates significantly from the trend of the data.

The curve fit is a approximation method for finding the best fit for the data points and due to which there are inherent errors in the process, these errors can be computed and treated as handles for controlling the parameters of the matlab program to obtain the best fit. This is done by evaluating the goodness of the fitted curve.

After fitting data with one or more models, next step is to evaluate the goodness of fit. A visual examination of the fitted curve can be performed in the Curve Fitting Toolbox. Understanding the working of the toolbox that provides these goodness of fit measures for both linear and nonlinear parametric fits and carrying out the evaluation methods manually is essential. The goodness of the fit can be evaluated by the following measures:

1. Residuals

2. Goodness of fit statistics

3. Confidence and prediction bounds

The residuals and prediction bounds are graphical measures, while the goodness of fit statistics and confidence bounds are numerical measures. Graphical measures are more beneficial than numerical measures because they allow the entire data set to view at once, and they can easily display a wide range of relationships between the model and the data. The numerical measures are more narrowly focused on a particular aspect of the data and often try to compress that information into a single number. In practice, depending on the data points and analysis requirements, both types of measures can be used to determine the best fit.

It is possible that none of the fits performed can be considered the best one. In this case, it might be that there is a need to select a different model. Conversely, it is also possible that all the goodness of fit measures indicate that a particular fit is the best one.

Residuals:

The residuals from a fitted model are defined as the differences between the response data and the fit to the response data at each predictor value, `residual (r) = data - fit``

Mathematically, the residual for a specific predictor value is the difference between the response value `y` and the predicted response value `hat y`, `r = y - hat y`. If the model that is fit to the data is correct, the residuals approximate the random errors. Therefore, if the residuals appear to behave randomly, it suggests that the model fits the data well. However, if the residuals display a systematic pattern, it is a clear sign that the model fits the data poorly.

Goodness of Fit Statistics:

Apart from the graphical methods to evaluate the goodness of fit, the goodness of fit statistics can be examined out for optimising the curve fit. The Curve Fitting Toolbox supports different measures of goodness of fit statistics for parametric models, they are:

1. Sum of squares due to error (SSE)

2. R-square

3. Adjusted R-square

4. Root mean squared error (RMSE)

Sum of Squares Due to Error:

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(y_i - haty_i)^2`, a value closer to 0 indicates a better fit.

R-Square:

This statistic measures how successful the fit is in explaining the variation of the data. Put another way, 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. R-square 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=0)^n w_i(haty_i - bary)^2`, where `bar y` is the mean of the response value `y`.

SST is also called the sum of squares about the mean, and is defined as: `SST = sum_(i=0)^n w_i(y_i - bary)^2`, where SST = SSR + SSE.

R-square is expressed as: R-square = `(SSR)/(SST) = 1 - (SSE)/(SST)`, R-square can take on any value between 0 and 1, with a value closer to 1 indicating a better fit. For example, an R-square value of 0.8234 means that the fit explains 82.34% of the total variation in the data about the average.

If there is increase in the number of fitted coefficients in the selected model, R-square might increase although the fit may not improve. To avoid this situation, there is a need to use the degrees of freedom adjusted R-square statistic.

Degrees of Freedom Adjusted R-Square:

This statistic uses the R-square statistic defined above, and adjusts it based on the residual degrees of freedom. The residual degrees of freedom is defined as the number of response values `n` minus the number of fitted coefficients `m` estimated from the response values.

`v = n - m`, where `v` indicates the number of independent pieces of information involving the `n`
data points that are required to calculate the sum of squares. The adjusted R-square statistic is generally the best indicator of the fit quality when additional coefficients are included in the model.

Adjusted R-square = `1 - (SSE(n-1))/(SST(v-1))`

The adjusted R-square statistic can take on any value less than or equal to 1, with a value closer to 1 indicating a better fit.

Root Mean Squared Error:

This statistic is also known as the fit standard error and the standard error of the regression,

`RMSE = s = sqrt(MSE)`, where `MSE` is the mean square error or the residual mean square.

`MSE = (SSE)/v`, A RMSE value closer to 0 indicates a better fit.

Confidence and Prediction Bounds:

Confidence bounds are calculated for the fitted coefficients, and prediction bounds are for new observations or for the fitted function. The confidence bounds are numerical, while the prediction bounds are displayed graphically.

Confidence and prediction bounds define the lower and upper values of the associated interval, and define the width of the interval. The width of the interval indicates how uncertain you are about the fitted coefficients, the predicted observation, or the predicted fit. For example, a very wide interval for the fitted coefficients can indicate that you should use more data when fitting before you can say anything very definite about the coefficients and vice-versa. The bounds are defined with a level of certainty that you specify. The level of certainty is often 95%, but it can be any value such as 90%, 99%, 99.9%, and so on. For example, you might want to take a 5% chance of being incorrect about predicting a new observation. Therefore, you would calculate a 95% prediction interval. This interval indicates that you have a 95% chance that the new observation is actually contained within the lower and upper prediction bounds.

Now for understanding the maths behind the curve fitting process, a general polynomial equation of `j th` order can be used. This case will specifically be treated as polynomial curve fitting. The polynomial curve fitting can be solved by gaussian elimination method.

The general form of the polynomial of order `j` can be written as:

The curve that gives minimum error between data `y` and the fit `f(x)` is ‘best’.

By least error squares approach method,

The general expression for any error using the least squares approach is

,

where we want to minimize this error. Now substituting the general polynomial form in the equation above,

into the genral least squares error equation,

, where `n` is the number of data points given, `i` is the current data point being summed and `j` is the polynomial order, re-writing the equation,

, the best line will have minimum error(squared distance) between the line and the data points. The co-efficients `a_k, a_0` are found out first, to minimize the error equation above, we take derivative with respect to each co-efficient `a_0,...,a_k` set to zero, where `k = 1,....,j` ,

re-writing these `j+1` equations, and puttting this in matrix form,

, where all summations are over `i = 1,....,n`. Data points `(x_i,y_i) for i = 1,....,n` is known. Co-efficients `a_0,...,a_k`, `k = 1,....,j` can be found out by using gaussian elimination method as,

where the summations above are over `i = 1,...., n` data points. The order of the polynomial can be varrying as per models selected but we get the linear equations at the end which can be solved for obtaining the solution vector `X` that contains the co-efficients of the polynomial,

`AX = B`, which can be re-written as `X = A^(-1)*B`, this can be written in matlab as `X =` inv`(A)*B`

Objectives of the project:-

1. To analyse the case study problem statement and to create a program that creates different curve fit models for the original data set.

2. Plotting the curve fit models along with the original data set curve to compare for the better fit.

3. Perform evaluation of the goodness of the curve fit for the selected fitted models.

4. Perform methods to find the highest order of polynomial that would fit the original data set with lowest possible error.

5. Validate the outcomes of the curve fit models created in this matlab program with the standard results in the curve fitting toolbox in matlab.

Tools Used:- Matlab

Case Study:- The following matlab program is focussed on the variation of the specific heat (`C_p`) with the variation in temperature and for this a set of 3200 original data points are predefined in the program to compare and create different curve fit models that would best fit the original data curve with lowest possible error.

Before jumping straight to using the matlab built-in library functions to carry out the curve fitting process, the same process can be carried out manually by mathemetical computation as discussed earlier in the theory section.

The matlab program explaining the mathematical computation for the 1st, 2nd and 3rd order polynomial curve fits i.e. linear, quadratic and cubic curve fit respectively are explained below in piecewise with the code snippet and their output plots,

The code snippet of the header lines of the matlab program is written as,

%Aim of the program is to carry out the curve fitting using mathemetical computation

clear all
close all
clc

%DOwnloading the data in this program
cp_data = load('data');
temperature = cp_data(:,1);
cp = cp_data(:,2);

%Finding out the correct co-efficients/constants of the relation of the parameters of
%the original data set using the gaussian elemination method
%Basis of gaussian elimination method: X = inv(A)*B
%Considering temperature as my x axis data and cp as my y axis data

The code snippet for linear fit is written in matlab as,

%%For a first order polynomial equation i.e linear polynomial equation
disp('For a first order polynomial equation i.e linear polynomial equation')
%Elements of matrix A
disp('Elements of matrix A')
format long g
f_n = length(cp);  %To know the number of data sets
f_x = sum(temperature)
f_x_2 = sum(temperature.^2)

%Elements for matrix B
disp('Elements for matrix B')
f_y = sum(cp)
f_x1_y1 = sum(temperature.*cp)

%Matrix A
disp('Matrix A')
format short e
f_A = [f_n, f_x;f_x, f_x_2]

%Matrix B
disp('Matrix B')
f_B = [f_y; f_x1_y1]

%Elements for solution matrix
%f_sol = [f_sol(1), f_sol(2)]

%The constants of the first order polynomial are
disp('The constants of the first order polynomial are')
format long g
f_sol = inv(f_A)*f_B

disp('Co-efficient of x^1')
f_a = f_sol(2)

disp('Co-efficient of x^0')
f_b = f_sol(1)

%For curve fitting
f_t_c = linspace(temperature(1),temperature(end),length(temperature)); 
%f_t_c is same as t, just for uniformity with original and curve fit data
%set, f_t_c is choosen
f_cp_c = f_a.*f_t_c + f_b;

%compare the curve fit with the original data 
figure(1)
plot(temperature,cp,'linewidth',3)  %original plot between temperature and cp values
hold on
plot(f_t_c, f_cp_c, 'linewidth',3,'color','r') %curve fit relation using method 2
xlabel('temperature (K)')
ylabel('specific heat (KJ/Kmol-K)')
title('Curve fit for original data set for 1st order polynomial')
legend('original data curve','linear curve fit')

The output of the code snippet in command window can be seen as,

The output plot of the code snippet can be seen as below,

Inference - The linear curve fit is a simplest curve fit and hence it does not comply with the original original data set well.

The code snippet for quadratic fit is written in matlab as,

%%For a second order polynomial equation i.e quadratic polynomial eqaution
disp(' ')
disp('For a second order polynomial equation i.e quadratic polynomial eqaution')
%Elements of matrix A
disp('Elements of matrix A')
format long g
n = length(cp);  %To know the number of data sets
x = sum(temperature)
x_2 = sum(temperature.^2)
x_3 = sum(temperature.^3)
x_4 = sum(temperature.^4)

%Elements for matrix B
disp('Elements for matrix B')
y = sum(cp)
x1_y1 = sum(temperature.*cp)
x1_2_y1 = sum((temperature.^2).*cp)

%Matrix A
disp('Matrix A')
format short e
A = [n, x, x_2;x, x_2, x_3;x_2, x_3, x_4]

%Matrix B
disp('Matrix B')
B = [y; x1_y1; x1_2_y1]

%Elements for solution matrix
%sol = [sol(1), sol(2), sol(3)]

%The constants of the second order polynomial are
disp('The constants of the second order polynomial are')
format long g
sol = inv(A)*B

disp('Co-efficient of x^2')
a = sol(3)

disp('Co-efficient of x^1')
b = sol(2)

disp('Co-efficient of x^0')
c = sol(1)

%For curve fitting
t_c = linspace(temperature(1),temperature(end),length(temperature)); 
%t_c is same as t, just for uniformity with original and curve fit data
%set, t_c is choosen
cp_c = a.*t_c.^2 + b.*t_c + c;

%compare the curve fit with the original data 
figure(2)
plot(temperature,cp,'linewidth',3)  %original plot between temperature and cp values
hold on
plot(t_c, cp_c, 'linewidth',3,'color','r') %curve fit relation using method 2
xlabel('temperature (K)')
ylabel('specific heat (KJ/Kmol-K)')
title('Curve fit for original data set for 2nd order polynomial')
legend('original data curve','quadratic curve fit')

The output of the code snippet in command window can be seen as,

The output plot of the code snippet can be seen as below,

Inference - The quadratic curve fit does a better job than the linear curve fit by including more and more original data points into consideration.

The code snippet for cubic fit is written in matlab as,

%%For a third order polynomial equation i.e cubic polynomial eqaution
disp(' ')
disp('For a third order polynomial equation i.e cubic polynomial eqaution')
%Elements of matrix A
disp('Elements of matrix A')
format long g
a_n = length(cp);  %To know the number of data sets
a_x = sum(temperature)
a_x_2 = sum(temperature.^2)
a_x_3 = sum(temperature.^3)
a_x_4 = sum(temperature.^4)
a_x_5 = sum(temperature.^5)
a_x_6 = sum(temperature.^6)

%Elements for matrix B
disp('Elements for matrix B')
b_y = sum(cp)
b_x1_y1 = sum(temperature.*cp)
b_x1_2_y1 = sum((temperature.^2).*cp)
b_x1_3_y1 = sum((temperature.^3).*cp)

%Matrix A
disp('Matrix A')
format short e
A1 = [a_n, a_x, a_x_2, a_x_3;a_x, a_x_2, a_x_3, a_x_4;a_x_2, a_x_3, a_x_4, a_x_5;a_x_3, a_x_4, a_x_5, a_x_6]

%Matrix B
disp('Matrix B')
B1 = [b_y; b_x1_y1; b_x1_2_y1; b_x1_3_y1]

%Elements for solution matrix
%sol1 = [sol1(1), sol1(2), sol1(3), sol1(4)]

%The constants of the second order polynomial are
disp('The constants of the second order polynomial are')
format long g
sol1 = inv(A1)*B1

disp('Co-efficient of x^3')
a1 = sol1(4)

disp('Co-efficient of x^2')
b1 = sol1(3)

disp('Co-efficient of x^1')
c1 = sol1(2)

disp('Co-efficient of x^0')
d1 = sol1(1)

%For curve fitting
t_c1 = linspace(temperature(1),temperature(end),length(temperature)); 
%t_c1 is same as t, just for uniformity with original and curve fit data
%set, t_c1 is choosen
cp_c1 = a1.*t_c1.^3 + b1.*t_c1.^2 + c1.*t_c1 + d1;

%compare the curve fit with the original data
figure(3)
plot(temperature,cp,'linewidth',3)  %original plot between temperature and cp values
hold on
plot(t_c1, cp_c1, 'linewidth',3,'color','r') %curve fit relation using method 2
xlabel('temperature (K)')
ylabel('specific heat (KJ/Kmol-K)')
title('Curve fit for original data set for 3rd order polynomial')
legend('original data curve','cubic curve fit')

The output of the code snippet in command window can be seen as,

The output plot of the code snippet can be seen as below,

Inference - The cubic curve fit does a further better job than the quadratic curve fit by including more and more original data points into consideration.

Now, the following code snippet is for the same matlab program written above but while using matlab built - in library functions such as polyfit, polyval etc,

Note- Two methods has been used for better understanding of the matlab library functions, method 1 is short form of the use of library function and method 2 is the expanded form of the library function for understanding and the output result is also compared on the same basis. The purpose of comparing method 1 and method 2 is for checking the exactness of the output results. The following code snippet contains the code for evaluating the goodness of fit characteristic measure as described earlier in the theory section. The following code also contains a code section explaining the accuracy of the selected co-efficients to original co-efficients for the data set, closer the values of the co-efficient correction matrix of the fitted model with the original data set, the more accurate are the selected co-efficients for the curve fit.

The following code snippet is a generic code i.e. the code can be used for all orders of polynomial by simply changing the input value of d in the code (i.e. degree of the polynomial), for iilustration purposes the degree of polynomial for curve fit assumed is 2 i.e. quadratic curve fit for comparing method 1 and method 2.

%Aim of the program is to perform curvefit on the given set of data 

clear all
close all
clc

%DOwnloading the data in this program
cp_data = load('data');
temperature = cp_data(:,1);
cp = cp_data(:,2);

%performing curve fit
%considering a polynomial function for the given problem i.e. variation of
%cp value with temperature
disp('Degree of the polynomial: ')
d = 2   %degree of the polynomial
co_effs = polyfit(temperature, cp, d);

%method 1 - short form using polyval
predicted_cp = polyval(co_effs,temperature);
disp(' ')
disp('The co_efficients of the the polynomial are: ')
disp(co_effs')
err = predicted_cp - cp;
disp(' ')
disp('Mean of error between predicted cp and original cp data set using method 1: ')
err_mean = mean(err)

%method 2 - expanded form as the basis of polyval function
disp(' ')
disp('Co-efficients calculated by method 2: ')
a = co_effs(1)
b = co_effs(2)
c = co_effs(3)
t_c = linspace(temperature(1),temperature(end),length(temperature)); 
%t_c is same as t, just for uniformity with original and curve fit data
%set, t_c is choosen
cp_c = a.*t_c.^2 + b.*t_c + c;

%sum of squares due to error
disp(' ')
disp('Sum of squares due to error by method 1:')
SSE = sum(err.^2) 
MSE = mean(err.^2)  %mean of the squared errors
RMSE = sqrt(MSE)    %root of the mean of squared errors

%sum of squares due to error
disp(' ')
disp('Sum of squares due to error by method 2:')
err1 = (cp_c)' - cp;
disp(' ')
disp('Mean of error between predicted cp and original cp data set using method 2: ')
err1_mean = mean(err1)
SSE1 = sum(err1.^2)        
MSE1 = mean(err1.^2)  %mean of the squared errors
RMSE1 = sqrt(MSE1)    %root of the mean of squared errors

%sum of squares of the regression
disp(' ')
disp('Sum of squares of the regression: ')
cp_mean = mean(cp);
s = predicted_cp - cp_mean;
SSR = sum(s.^2)     

%total sum of squares
disp(' ')
disp('Total sum of squares: ')
SST = SSR + SSE     

%R_square
disp(' ')
disp('R square is calculated as: ')
R_square = SSR/SST

%Degrees of Freedom Adjusted R-Square
disp(' ')
disp('Degrees of Freedom Adjusted R-Square: ')
n = length(cp)      %number of response values
m = length(co_effs) %number of co-efficients
v = n-m             %degrees of freedom
adjusted_R_square = 1 - ((SSE*(n-1))/(SST*(v-1)))


%Comparing the accuracy of selected co-efficients to original co-efficients
%of the relation of the data set 
%original data set
disp(' ')
disp('co-efficient correction for original data set: ')
original_relationship_co_efficients = corrcoef(temperature,cp)

%Method 1 co-efficient relationship
disp(' ')
disp('co-efficient correction for curve fit co-efficients using method 1: ')
curvefit_method1_relationship_co_efficients = corrcoef(temperature,predicted_cp) 

%Method 1 co-efficient relationship
disp(' ')
disp('co-efficient correction for curve fit co-efficients using method : ')
curvefit_method2_relationship_co_efficients = corrcoef(temperature,cp_c)


%compare the curve fit with the original data using method 1
figure(1)
plot(temperature,cp,'linewidth',3)  %original plot between temperature and cp values
hold on
plot(temperature, predicted_cp, 'linewidth',3,'color','r')  %curve fit relation using method 1
xlabel('temperature (K)')
ylabel('specific heat (KJ/Kmol-K)')
title(['Curve fit for original data set using method 1 for order - ',num2str(d),' polynomial'])
legend('original data curve','curve fit')

%compare the curve fit with the original data using method 2
figure(2)
plot(temperature,cp,'linewidth',3)  %original plot between temperature and cp values
hold on
plot(t_c, cp_c, 'linewidth',3,'color','r') %curve fit relation using method 2
xlabel('temperature (K)')
ylabel('specific heat (KJ/Kmol-K)')
title(['Curve fit for original data set using method 2 for order - ',num2str(d),' polynomial'])
legend('original data curve','curve fit')

The output of the following code snippet can be seen in the command window as,

The output plots can be seen as below for the above code snippet,

The same matlab program can be run for higher order polynomial by providing input to 'd' as 3,4,5 and so on..., some of the higher order curve fits can seen as below,

From these we get an important inference that as we increase the order of the polynomial, we get a better curve fit for the original data points, which brings the point of finding out the optimum and highest order of polynomial that could be the best fit for the data points. As discussed earlier in theory section, the goodness of fit measures indicate that a particular fit is the best one. Hence instead of using this evaluation of fit measure to calculate the errorrs between the predicted values and original values, conversely this error calculation can be used to determine the highest and optimum order of polynomial that could provide a curve fit which is having 0.1% variation, or its acccuracy is 99.9%. The following matlab code snippet describes the calculation of such order of the polynomial as below,

%Aim of the program is to find the most accurate degree of polynomial for  99.9% polynomial curve fit to the original data set

clear all
close all
clc

%DOwnloading the data in this program
cp_data = load('data');
temperature = cp_data(:,1);
cp = cp_data(:,2);

flag = 0;
for i = 1:inf
    d = i;   %degree of the polynomial
    co_effs = polyfit(temperature, cp, d);
    predicted_cp = polyval(co_effs,temperature);
    err = predicted_cp - cp;
    %sum of squares due to error
    SSE = sum(err.^2); 
    MSE = mean(err.^2);  %mean of the squared errors
    RMSE = sqrt(MSE);    %root of the mean of squared errors 
    if RMSE < 0.1     %Error in the selection of co-efficients is less than 0.1% 
        flag = 1;
    end
    if flag == 1
            RMSE
            degree_of_polynomial = d
            break
    end
end

%compare the curve fit with the original data 
figure(1)
plot(temperature,cp,'linewidth',3)  %original plot between temperature and cp values
hold on
plot(temperature, predicted_cp, 'linewidth',3,'color','r')  %curve fit relation 
xlabel('temperature (K)')
ylabel('specific heat (KJ/Kmol-K)')
title(['Curve fit for original data set by using order - ',num2str(d),' polynomial'])
legend('original data curve','curve fit')

The output value of the RMSE that is less than 0.1 and the order of the polynomial that could possibly best fit the curve can be seen in the command window as below,

And the output curve fit plot can be seen as below,

The curve fits are also used for estimating or predicting the behaviour of the case study at higher or extended data points. A simple observation of the same can be performed in matlab by writing the following matlab code snippet to predict the values of `C_p` with variation of data points of temperature extended from 3200 K to 7000 K.

%AIm of the program is to predict the behaviour of the curve fit for higher extended values

clear all
close all
clc

%DOwnloading the data in this program
cp_data = load('data');
temperature = cp_data(:,1);
cp = cp_data(:,2);

%estimating higher values
temp = temperature';
temp_1 = 3201:1:7000;
temp_2 = [temp temp_1];
temp_extend = temp_2';
co_effs1 = polyfit(temperature, cp, 3);
predicted_cp_extend = polyval(co_effs1,temp_extend);

plot(temperature,cp,'linewidth',3)  %original plot between temperature and cp values
hold on
plot(temp_extend, predicted_cp_extend, 'linewidth',3,'color','r')  %curve fit relation 
xlabel('temperature (K)')
ylabel('specific heat (KJ/Kmol-K)')
legend('original data curve','curve fit')

The output curve can be seen as below,

Note - The above code is written for 3rd order polynomial as for higher order i.e. 4, 5 and so on, the curve becomes unstable that can be seen as below which is not much of use for prediction purposes. Hence for predicting the behaviour of the curve for extended values, 3rd order polynomial serves useful one.

As discussed in theory before, the prediction intervals for the curve fits can be found out as well by writing the following code snippet,

%Aim of the program is to perform curvefit on the given set of data with prediction interval

clear all
close all
clc

%DOwnloading the data in this program
cp_data = load('data');
temperature = cp_data(:,1);
cp = cp_data(:,2);

%performing curve fit
%considering a polynomial function for the given problem i.e. variation of
%cp value with temperature
d = 2;  %degree of the polynomial
[co_effs, S] = polyfit(temperature, cp, d);
[predicted_cp, delta] = polyval(co_effs, temperature, S);
disp(co_effs')
err = predicted_cp - cp;
MSE = mean(err.^2)
RMSE = sqrt(MSE)

%compare the curve fit with the original data
plot(temperature,cp,'linewidth',3)  %original plot between temperature and cp values
hold on
plot(temperature, predicted_cp, 'linewidth',3,'color','r')
plot(temperature,predicted_cp+2*delta,'m--',temperature,predicted_cp-2*delta,'m--')
xlabel('temperature (K)')
ylabel('specific heat (KJ/Kmol-K)')
legend('original data curve','curve fit','95% Prediction Interval')
title(['Prediction interval for curve fits for order - ',num2str(d),' polynomial'])

The sample output curve fit plots with prediction interval can be seen as below,

As discussed in the theory section, residuals are a part of the evaluation of the goodness of the curve fit, for better visualisation on how the data points are scattered around the curve fit line or the residuals of the curve fit can be seen as below by the following code snippet,

%Aim of the program is evaluate and plot the residuals of the curve fit

clear all
close all
clc

%DOwnloading the data in this program
cp_data = load('data');
temperature = cp_data(:,1);
cp = cp_data(:,2);

%Sorting the array as in general before plotting for any case study
[temperature, sortOrder] = sort(temperature, 'ascend'); 
cp = cp(sortOrder); 

disp('Degree of the polynomial: ')
d = 12   %degree of the polynomial
co_effs = polyfit(temperature, cp, d);
predicted_cp = polyval(co_effs,temperature);

plot(temperature, cp, 'linewidth',10,'color','r');
hold on;
plot(temperature, predicted_cp, 'LineWidth', 6,'color','g');
grid on
xlabel('temperature (K)')
ylabel('specific heat (KJ/Kmol-K)')

% Plot residuals as lines from actual data to fitted line.
for k = 1 : length(temperature)
  yActual = cp(k);
  yFit = predicted_cp(k);
  x = temperature(k);
  line([x, x], [yFit, yActual], 'Color', 'm');
end

legend('original data curve','curve fit','residuals')
title(['Residuals of the curve fit for order - ',num2str(d),' polynomial'])

The sample output plot for the residual of curve fits of order - 1, 2, 3 and 12 can be seen as below,

The results of the above calculations for generating various curve fit models and curve fit evaluation measures are validated with the standard results by comparing with the results obtained in the curve fitting toolbox. The same program is run in curve fitting toolbox for order of the polynomial - 1(linear fit), 2(quadratic fit) and 3(cubic fit) and the results can be seen below,

Hence it can be seen that the operations performed in obtaining the curve fit are in correct accordance as the observed results are same as the standard results.

Conclusion:- The objectives stated at the start of the project are extensively carried out and the results obtained are positive. The best curve fit is modeled and various evaluation measures are carried out, validation of the curve fit models are carried out as well. Curve fitting is a very powerful tool that helps in analyzing the behaviour or the trend or the nature of variation of the data points, it can also help in forecasting methods and much more.

Refernce:-

1. Skill lync basic curve fit tutorial video.

2. Multiple documents on curve fitting.

3. Matlab Documentation.