Linear & Cubic Curve fitting for Specific Heat & Temperature data

OBJECTIVE:

To write a MatLab program(s) to fit a linear and cubic polynomial for the given Specific Heat data.

Compare the linear and cubic fit curves along with the raw data points by plotting them together.

OVERVIEW:

Curve fitting is a process of constructing curve or mathematical functions, that is close to the original dataset. The data may be scattered. By curve fitting, we can match the trend of the observed dataset. It is beneficial and highly effective in deducing a mathematical model of naturally happening processes.

Here, we will be using the regression analysis method. It is a process of estimation of the relationship between the dependent variable and one or more independent variables. Linear regression and Polynomial regression are the most common forms of regression analysis. In Linear regression, a line (or a more complex linear function) is found which most closely fits the dataset. Polynomial regression is used to determine the relationship between the independent variable and the dependent variable at an nth degree polynomial.

Fitting Lines and Polynomial functions to data points:

1. The first-degree polynomial equation/function, `f(x)=y=ax+b`, is a line with a slope 'a'. A line will connect any two points, so a first-degree polynomial equation is an exact fit through any two points with distinct 'x' coordinates.
2. The second-degree polynomial function, `f(x)=y=ax^2+bx+c`, will exactly fit a simple curve to three points, with distinct 'x' coordinates.
3. Similarly third-degree polynomial function, `f(x)=y=ax^3+bx^2+cx+d`, will exactly fit four points. A more general statement would be to say it will exactly fit four constraints. Each constraint can be a point, angle, or curvature. Angle and curvature constraints are most often added to the ends of a curve, and in such cases are called end conditions.
   
   

In the figure above, the red line is a first-degree polynomial, the green line is second degree, the orange line is the third-degree and blue is a fourth-degree polynomial.

The first-degree polynomial equation could also be an exact fit for a single point and an angle while the third-degree polynomial equation could also be an exact fit for two points, an angle constraint, and a curvature constraint.

Several reasons are given to get an approximate fit when it is possible to simply increase the degree of the polynomial equation and get an exact match:

1. High-order polynomials can be highly oscillatory. If a curve runs through two points A and B, it would be expected that the curve would run somewhere near the midpoint of A and B, as well. This may not happen with high-order polynomial curves; they may even have values that are very large in positive or negative magnitude.
2. With low-order polynomials, the curve is more likely to fall near the midpoint (it's even guaranteed to exactly run through the midpoint on a first degree polynomial). Low-order polynomials tend to be smooth and high order polynomial curves tend to be "lumpy". To define this more precisely, the maximum number of inflection points possible in a polynomial curve is n-2, where n is the order of the polynomial equation. An inflection point is a location on the curve where it switches from a positive radius to negative.
   
   The degree of the polynomial curve being higher than needed for an exact fit is undesirable for all the reasons listed previously for high order polynomials.
   
   

Goodness of a fit:

The goodness of fit of a mathematical model describes how well it fits a set of observations. Measures of goodness of fit typically summarize the discrepancy/deviation between observed values and the values expected under the model in question.

In regression analysis, the coefficient of determination (`R^2`- r squared) is the proportion of the variance/deviations in the dependent variable, i.e., the predictable from the independent variables.

A data set has n values marked y₁,...,y_n (collectively known as y_i), each associated with a fitted (or modeled, or predicted) value f₁,...,f_n (known as f_i, sometimes as ŷ_i).

If `bary` is the mean of the observed data:

then the variability of the data set can be measured using three sums of squares formulas:

1. The total sum of squares (TSS) measures how much variation there is in the observed data, and to the residual sum of squares.
   `TSS=RSS+ESS`
   
   
2. The residual sum of squares (RSS) is a measure of the discrepancy between the data and an estimation model. A small RSS indicates a tight fit of the model to the data.
   
   
   
3. The explained sum of squares (ESS), also called the regression sum of squares, measures how much variation there is in the modeled values.
   
   
   `therefore` By definition, the Coefficient of determination is given by the formula:
   
   `R^2=(RSS)/(TSS)=(RSS)/(ESS+RSS)`
   
   If the value of `R^2` is more than 0.95, the Goodness of the fit is best for the given dataset.

MATLAB PROGRAM:

Progam-1 & explanation:

% Linear & Cubic CURVE FIT 

close all 
clear all
clc

% Preparing data
temp_cp_data = load('data');
t = temp_cp_data(:,1);
cp = temp_cp_data(:,2);


% Linear Polynomial
co_eff1 = polyfit(t,cp,1);
cp1 = polyval(co_eff1,t);


%Cubic Polynomial Co-efficients
[co_eff2,s,mu] = polyfit(t,cp,3);
cp2 = polyval(co_eff2,t,s,mu);

% Plot comparison

figure(1)
plot(t,cp,'color','r','linewidth',2);
hold on 
plot(t,cp1,'color',[0, 0.75, 0.75],'linewidth',2)
axis([200 3600 900 1500])
plot(t,cp2,'color','b','linewidth',2)
axis([200 3600 900 1500])
xlabel('Temperature [K]');
ylabel('Specific Heat [JK^(-1) (or) Kmol-K]');
legend('Original dataset','Linear Curvefit','Cubic Curvefit');


% Error Calculation
mean = sum(cp)/length(cp);


rss_linear = sum((cp1 - mean).^2);
ess_linear = sum((cp - cp1).^2);
Linear_rsquared = rss_linear/(rss_linear+ess_linear)


rss_cubic = sum((cp2 - mean).^2);
ess_cubic = sum((cp - cp2).^2);
Cubic_rsquared = rss_cubic/(rss_cubic+ess_cubic)

• At first the file containing the Temperature & corresponding specific heat data is loaded into the program using the load('data') command and it is stored in a variable 'temp_cp_data'. Then the temperature data is extracted from 'temp_cp_data' and stored in a variable using the command, t=temp_cp_data(:,1). Similarly, the specific heat data is extracted from 'temp_cp_data' and stored in a variable using the command cp=temp_cp_data(:,2).
  
• As explained before, a linear polynomial has two coefficients. So the best linear coefficients for the dataset is found and is stored in a variable using the command co_eff1=polyfit(t,cp,1). Using the 'co_eff1 & t' values, an array of predicted specific heat values are generated using and stored in a variable using the command cp1=polyval(co_eff1,t).
• Similarly, a cubic polynomial (third-degree) has four coefficients. The best cubic coefficients for the dataset is found and is stored in a variable using the command [co_eff2,s,mu]=polyfit(t,cp,3). Using the 'co_eff2 & t' values, an array of predicted specific heat values are generated using and stored in a variable using the command cp2=polyval(co_eff2,t,s,mu). Here using 's & mu' the polyfit() command returns Structure and a two-element vector with centering and scaling values. This centering and scaling transformation improves the numerical properties of both the polynomial and the fitting algorithm.
• Using the original dataset (t,cp), the predicted values (cp1) for linear curve-fit and the predicted values (cp2) for the cubic curve-fit, three plots are generated in a single figure window using plot() and hold on command. Different colors are defined for the three varied plots. X & Y axes are labeled using xlabel & ylabel. Axis range is defined using axis() command. The three curves are labeled using the legend() command.
• The `R^2` value is calculated by defining the formulae: mean, RSS, ESS & `R^2`.

Program-1 Output:

From the above output snippet, it can be seen that the `R^2` value for Linear Curvefit is lower than 0.95, which is not a good fit for the original dataset.

The linear and cubic curve-fits in the above graph is the actual representation of the corresponding `R^2` values. The cubic curve-fit is very close to the original dataset, whereas the linear curve-fit only intercepts the original curve at 2 points just like the values suggest.

However, the goodness of neither of the curve-fits is the best-fit.

Program-2 & explanation:

Since the previous method did not serve well in obtaining the best-fit even with the cubic polynomial, I will be using the splitwise curve-fit method with a quadratic polynomial in this part to construct the best-fit curves for the given dataset.

% Linear & Quadratic CURVEFIT

close all 
clear all
clc

% Preparing data
temp_cp_data = load('data');
t = temp_cp_data(:,1);
cp = temp_cp_data(:,2);

% Plotting Original T & CP data
plot(t,cp,'color','red','linewidth',2)
hold on

a = 1;
b = 400;

for i = 1:8
    coeff(:,i) = polyfit(t(a:b),cp(a:b),2);
    cp_predicted(:,i) = polyval(coeff(:,i),t(a:b));
    
    % Compare plot b/w Original data and curvefit
    plot(t(a:b),cp_predicted(:,i),'color','b','linewidth',2);
    axis([200 3600 900 1400])
    xlabel('Temperature [K]')
    ylabel('Specific Heat [J/K (or) Kmol-K]')
    legend('Original Dataset','Curvefit')
    
    CP(a:b,1) = cp_predicted(:,i);
    
    a = a+400;
    b = b+400;     
end


% Error Calculation
mean = sum(cp)/length(cp);

rss = sum((CP - mean).^2);
ess = sum((cp - CP).^2);
r_squared = rss/(rss+ess)

• At first the file containing the Temperature & corresponding specific heat data is loaded into the program using the load('data') command and it is stored in a variable 'temp_cp_data'. Then the temperature data is extracted from 'temp_cp_data' and stored in a variable using the command, t=temp_cp_data(:,1). Similarly, the specific heat data is extracted from 'temp_cp_data' and stored in a variable using the command cp=temp_cp_data(:,2).
• A plot is generated for the original dataset and hold on command is used to keep the plot.
• A for loop is opened and a variable is used for iterations. This number of iterations is the number of splits of the predicted curve. In the above program, the number of iterations is 8 and is defined by using the range operator ':'. Right before the for loop, two variables a=1; b=400 are defined.
• coeff(:,i)=polyfit(t(a:b),cp(a:b),2) Since this is a quadratic polynomial function, this command essentially computes 3 best coefficients for the original dataset for eight iterations.
• cp_predicted(:,i) = polyval(coeff(:,i),t(a:b)) the polyval() command in this line generates 400 new predicted cp values for eight iterations.
• For each iteration, using the cp_predicted values, a plot is generated in the same figure window as the original plot, but with a different color to differentiate itself from the original data curve.
• Axis range is defined using axis() command. X & Y axes are labeled using xlabel & ylabel commands. The original curve and the predicted curve are labeled using the legend() command.
• For each iteration, the cp_predicted is stored in 400 rows and in a new column. Before the a & b variables are incremented, the cp_predicted values are stored as a column matrix and stored in a variable using the command, CP(a:b,1)=cp_predicted(:,i).
• 'a & b' variables are incremented by 400 before the end of each iteration. And the above process is repeated eight times.
• The `R^2` value is calculated by defining the following formulae: mean, RSS, ESS & `R^2`.
  

Program-2 Output:

The `R^2` value printed in the command window suggests that the generated splitwise curve-fit is a perfect-fit curve for the original dataset. It can be inferred that by using splitwise curve-fit, using a higher degree polynomials are not obligatory to construct a best-fit curve for any given dataset.

The quadratic curve-fit in the above graph is the absolute reflection of the computed `R^2` value printed in the command window.

References:

• Understanding of the basics of the curve fitting - YouTube video by a Professor from the University of Washington.
• Explanation of error calculation - ouTube video by a Professor from the University of Washington and Wikipedia.
• Few online publications.