MATLAB Program: Curve Fitting

Aim: Write a MATLAB program to fit a curve for given data set.

Objective:

• Write a code to obtain a linear and cubic polynomial for given data set. (here, Cp vs Temperature values (specific heat at constant pressure))
• Plotting the obtained curves.
• Studying fitness characteristics for the same.

Theory:

Curve fitting:

Several equations of different types can be obtained to express the given data approximately. But the problem is to find the equation of the curve of best fit which may be most suitable for predicting the unknown values. The process of finding such an equation by best fit is known as curve-fitting. There are four methods to do so:

1. Method of Least Squares
2. Method of Moments
3. Method of Group
4. Graphical Methods

Among these Method of Least Square, is probably, the best fit to a unique curve to a given data.

Method of Least Squares:

The method of least squares is a standard approach in regression nalysis to approximate solution of overdetermined systems by minimizing the sum of the residuals made in the results of every single equation.

`->` To fit a straight line `y = a +bx`

1. Substitute the observed set of n values in this equation.
2. Form normal equations for each constant `sum y = na + bsumx; sumxy = asumx + bsumx^2.`
3. Solve these normal equations as simultaneous equations for `a` and `b`.
4. Substitute the values of `a` & `b` in `y = a + bx`, which is the required best fit.

`->` To fit a parabola `y = a + bx + cx^2`

1. Solve the derived normal equations for this polynomial.
2. Substitute the values of `a, b, c.`
3. We have `y = a + bx + cx^2`, the required polynomial of best fit.

Measuring Fitness Characteristics:

After performing curve fitting, we need to check its fitness, i.e. how well our derived function is able to represent the data set. The following are most common methods to check fitness characteristics:

• The sum of squares due to error (SSE)
• R-square method
• Adjusted R-square method
• Root mean square error (RMSE)

1. Sum of squares due to error (SSE)

It is the sum of squares of results or deviations predicted from actual values of data. It is measure of the discrepancy between the data and an estimation model. A small SSE indicates a tight fit of the model to the data.

`SSE = sum_(i=1)^n`Error`(i)^2`

2. R-square is a statistical measure of how close the data are to the fitted regression line. R-square is always between 0 and 1. 0 indicates that the model explains none of the variability of the response data around its mean; 1 indicates that the model explains all the variability of the response data around its mean. In general, the higher the R-square, the better the model fits the real data.

`R^2 = (SSR)/(SST)`

where,

`SSR` = Sum of squares of regression

`SST` = `SSR` + `SST`

3. Adjusted R-square is a normalizing variation of R-square due to insignificant independant variables. Adjusted R-square disregards infinitesimally small increments to R-square and only considers the variable which actually contribute to the output of the system in a drastic way.

`R^2(adj) = 1- [((1-R^2)(n-1))/(n-k-1)]`

where,

`n` = number of data points

`k` = the number of independant

4. Root Mean Square Error (RMSE) is an estimate of the standard deviation of the random component in the data, and is defined as,

`RMSE = sqrt((SSE)/n)`

Program:

The program required to clear the above objectives is as shown. The program was written in MATLAB programming language.

%Curve Fitting
clear all
close all
clc

%data Preparation
cp_data = load('data');
temperature  = cp_data(:,1);
cp = cp_data(:,2);

%Curve fit

%Linear Curve Fitting 
%cp = a*T + b
[co_eff1,S,mu] = polyfit(temperature,cp,1);
predicted_cp1 = polyval(co_eff1,temperature,[],mu);

%Quadratic Curve Fitting
%cp = a*T^2 + b*T + c
[co_eff2,S,mu] = polyfit(temperature,cp,2);
predicted_cp2 = polyval(co_eff2,temperature,[],mu);

%Cubic Curve Fitting
%cp = a*T^3 + b*T^2 + c*T + d
[co_eff3,S,mu] = polyfit(temperature,cp,3);
predicted_cp3 = polyval(co_eff3,temperature,[],mu);

%Quartic Curve Fitting
%cp = a*T^4 + b*T^3 + c*T^2 + d*T + e
[co_eff4,S,mu] = polyfit(temperature,cp,4);
predicted_cp4 = polyval(co_eff4,temperature,[],mu);

%Compare curve fit with the original data and plot
figure(1)
plot(temperature,cp,'linewidth',5,'color','r')
hold on
plot(temperature,predicted_cp1,'linewidth',3,'color','b')
axis([0 4000 900 1400])
xlabel('Temperature {K}')
ylabel('Specific Heat {KJ/Kmol-K}')
legend('Original Data','Linear Curvefit','Location','Southeast')
title('Linear Curve Fit')
grid on
pause(4)

figure(2)
plot(temperature,cp,'linewidth',5,'color','r')
hold on
plot(temperature,predicted_cp2,'linewidth',3,'color','b')
axis([0 4000 900 1400])
xlabel('Temperature {K}')
ylabel('Specific Heat {KJ/Kmol-K}')
legend('Original Data','Quadratic Curvefit','Location','Southeast')
title('Quadratic Curve Fit')
grid on
pause(4)

figure(3)
plot(temperature,cp,'linewidth',5,'color','r')
hold on
plot(temperature,predicted_cp3,'linewidth',3,'color','b')
axis([0 4000 900 1400])
xlabel('Temperature {K}')
ylabel('Specific Heat {KJ/Kmol-K}')
legend('Original Data','Cubic Curvefit','Location','Southeast')
title('Cubic Curve Fit')
grid on
pause(4)

figure(4)
plot(temperature,cp,'linewidth',5,'color','r')
hold on
plot(temperature,predicted_cp4,'linewidth',3,'color','b')
axis([0 4000 900 1400])
xlabel('Temperature {K}')
ylabel('Specific Heat {KJ/Kmol-K}')
legend('Original Data','Quartic Curvefit','Location','Southeast')
title('Quartic Curve Fit')
grid on
pause(4)

%R^2 method to find the Fitness Characteristics of the curve
M = mean(cp);

for i = 1:length(cp)
    SSR1 = sum((predicted_cp1(i)-M)^2);
    SSE1 = sum((cp(i)-predicted_cp1(i))^2);
end
SST1 = SSR1 + SSE1;
R1 = SSR1/SST1;

for i = 1:length(cp)
    SSR2 = sum((predicted_cp2(i)-M)^2);
    SSE2 = sum((cp(i)-predicted_cp2(i))^2);
end
SST2 = SSR2 + SSE2;
R2 = SSR2/SST2;

for i = 1:length(cp)
    SSR3 = sum((predicted_cp3(i)-M)^2);
    SSE3 = sum((cp(i)-predicted_cp3(i))^2);
end
SST3 = SSR3 + SSE3;
R3 = SSR3/SST3;

for i = 1:length(cp)
    SSR4 = sum((predicted_cp4(i)-M)^2);
    SSE4 = sum((cp(i)-predicted_cp4(i))^2);
end
SST4 = SSR4 + SSE4;
R4 = SSR4/SST4;

printf("The R-square fitness for linear fit is %d \n",R1)
printf("The R-square fitness for quadratic fit is %d \n",R2)
printf("The R-square fitness for cubic fit is %d \n",R3)
printf("The R-square fitness for quartic fit is %d \n",R4)

Working of Code:

Luckily, we don't have to derive or code the normal equations and let's see how we evaded it. First, we imported the data file using load command (it contained data in tabular form). Then, segragated the Cp values and Temperature values and stored them in the required variables. For curve fitting we used the following commands:

1. polyfit: Returns the coefficients of a polynomial `p(x)` of degree `n` that minimizes the least square error of the fit to the points `[x,y]`. The polynomial coefficinets are returned in a row vector.

2. polyval: Evaluates the polynomial `p` at the specified values of `x`. If `x` is a vector or matrix, the polynomial is evaluated for each of the elements of `x`.

To determine fitness characteristics, we used R-square method. To evaluate this we used for loop (it helped find squared errors and simultaneously summing it). Hence, plots were generated and R-square values of each fit was displayed.

Outputs:

1. Linear Fit

2. Quadratic Fit

3. Cubic Fit

4. Quartic Fit

5. R-Squares

Result and Conclusion:

• The R-square fitness was found to increase with increase in degree of polynomial.
• Best fits have higher vaules of R-sqaure (i.e. difference between actual and predicted values are minimized) whereas perfect fit has R-square value 1 (the traced curve has exact values with zero deviation)
• For complex relations, we can split the data into several parts and evaluate fit for each part to fit the curve.

References:

• MATLAB for Mechanical Engineers - Course provided by Skill Lync
• GNU Octave
• Supporting Material: Higher Engineering Mathematics - BS Grewal
• Supporting Material: Evaluation of Goodness of Fit - Project by Thanuj @ Skill Lync Projects

Material: Files can be accessed here

---