Curve fitting using MATLAB

clear all
close all
clc

% data
load('data')
cp_data=load('data')
Temp=cp_data(:,1);             % T represents temperature
cp=cp_data(:,2);            % cp as specific heat


co_eff1= polyfit(Temp,cp,1)
predicted_cp1 = polyval(co_eff1, Temp)

co_eff3= polyfit(Temp,cp,3)
predicted_cp3 = polyval(co_eff3, Temp)

co_eff6 = polyfit(Temp,cp,6)
predicted_cp6 = polyval(co_eff6, Temp)


% curve fit
% cp = a*T^2 + b*T + c
%a,b,c=3 variables
% using linear fit
figure(1)
plot(Temp,cp,'linewidth',3)
hold on
plot(Temp,predicted_cp1,'linewidth',3)
xlabel(' [K]')
ylabel('Specific Heat [KJ/Kmol-K]')
legend('original dataset','linear fit')

% using cubic fit
figure(2)
plot(Temp,cp,'linewidth',3)
hold on
plot(Temp,predicted_cp3,'linewidth',3)
xlabel(' [K]')
ylabel('Specific Heat [KJ/Kmol-K]')
legend('original dataset','cubic fit')

% using hexa fit
figure(3)
plot(Temp,cp,'linewidth',3)
hold on
plot(Temp,predicted_cp6,'linewidth',3)
xlabel(' [K]')
ylabel('Specific Heat [KJ/Kmol-K]')
legend('original dataset','hexa fit')

% Error calculation for linear polynomial
Error1 = cp-predicted_cp1 % just the errors
error_sq1= Error1.^2; % r^2 element wise operator squares each element in the matrix

SSE1=sum(error_sq1) % sum of the squared errors or sse
n = length(cp)
k = 2
RMSE1= (SSE1/n)^0.5

% Error calculation for 8 degree polynomial
Error8 = cp-predicted_cp6
error_sq6 = Error8.^2;
SSE6= sum(error_sq6)
RMSE6 = (SSE6/n)^0.5

The above program is written and executed in MATLAB to create a curve fit for a given set of data.

EXPLANATION:

• The load command is used to load the data file. The data file contains the specific heat data for oxygen as a function for temperature.

• The first column is for temperature T=cp_data(:,1) % T represents temperature in Kelvin
  and the second column is for specific heat cp=cp_data(:,2) % cp as specific heat in KJ/mol-Kelvin.

• Polyfit command returns coefficient for the polynomial of degree n that would be the best fit fot the data in y.
  co_eff1= polyfit(T,cp,1) returns the values of coefficients that corresponds to the fitting curve of degree 1.The equation will be of form ax +b(a,b being the coefficients).

• Polyval command is used to plug in the values of temperature into the coefficients.
  For example , predicted_cp3 = polyval(co_eff3, Temp).Doing this creates the predicted_cp3 curve as per equation- `ax^3 +bx^2 +cx +d`, where a,b,c,d are coeffecients obtained by polfit command using degree 3.
  
  
• Using the plot command to plot the graph for temperature and specific heat of the original data set.
• Using the hold on command to hold the graph for the above plot.
• Again using the plot command for the curve fit plot(T,predicted_cp1,'linewidth',3).
  
  Similary the curve fit has been done for cubic polynomial & hexa polynomial.
• Here the curve fit is carried out for 3 cases- for 1st degree polynomial,3rd degree and 6th degree.

RESULTS:

Curve fit for 1 degree polynomial (figure 1)

Curve fit for cubic polynomial

Curve fit for hexa polynomial

The errors in curve fit are calculated to measure the accuracy of the fit.This is carried out by witing the following piece of program.

% Error calculation for linear polynomial
Error1 = cp-predicted_cp1 % just the errors
error_sq1= Error1.^2; % r^2 element wise operator squares each element in the matrix

SSE1=sum(error_sq1) % sum of the squared errors or sse
n = length(cp)
k = 2
RMSE1= (SSE1/n)^0.5

% Error calculation for 8 degree polynomial
Error8 = cp-predicted_cp6
error_sq6 = Error8.^2;
SSE8= sum(error_sq6)
RMSE8 = (SSE8/n)^0.5

At first,error is calculated by subtracting the fit value by the actual value.

Then,the difference is squared.This is done for an array of values,thus elementwise operator was used. The square of differences are added individually using sum operator.

Values in workspace

Now RSME is calculated for each SSE output.

RMSE = sqrt((sum_(i=1)^n(Y(i)-f(x(i)))^2)/n) = sqrt((SSE)/n)

The less the errors(RSMEs), the best will be our curve fit.

As,RMSE1-for 1st degree polynomial =25.991

RMSE6-for 6th degree polynomial=1.6136

Thus it is proven that the error is lesser for a 6th degree polynomial than a 1st degree polynomial curve fit.

PERFECT FIT :The perfect fit is achieved when the geometric distances from the given points to the fitting curve are zero . This means that the curve fits perfectly on the points . To get a perfect fit we can increase to a higher order of the polynomials , this method increases the precession and reduces errors .

BEST FIT :The best fit is achieved when the geometric distances from the given points to the fitting curve are minimized, in the least squares sense. Finding the best fit reduces to the minimization of the objective function where denotes the fitting curve (line, circle, ellipse, etc.).Upon finding a curve which satisfies the data points there are error quantities which help us measure the goodness of fit criteria or how well the equation is a representing the datapoints.

1)The sum of squares due to error (SSE)

2)R-square

3)Adjusted R-square

4)Root mean squared error (RMSE)

To improve the curve fit by a 3rd degree polynomial,and to further reduce the Maximum error Percentage, The entire Data can be strategically split at 4 points, i.e., at T=410, 1300, 2200, 3200

Program to fit a split-wise Polynomial for the Cp Data, in order to make the curve fit perfect.

clear all
close all
clc

% Preparing the Data
cp_data = load('data');
temperature = cp_data(:,1);
cp = cp_data(:,2);

% Splitting the curve with 4 polynomials

% Temperature and cp values Range[1-400]
cp1=cp_data(1:410,2);
temperature1=cp_data(1:410,1);
% Temperature and cp values Range[401-1200]
cp2=cp_data(411:1300,2);
temperature2=cp_data(401:1300,1);
% Temperature and cp values Range[1301-2200]
cp3=cp_data(1301:2200,2);
temperature3=cp_data(1201:2200,1);
% Temperature and cp values Range[2201-3200]
cp4=cp_data(2201:3200,2);
temperature4=cp_data(2201:3200,1);

% Curve fitting with a splitwise polynomial 
[coeffs1,S1,mu1]=polyfit(temperature1,cp1,3); % Polynomial of degree3
[coeffs2,S2,mu2]=polyfit(temperature2,cp2,3); % Polynomial of degree3
coeffs3=polyfit(temperature3,cp3,2); % Polynomial of degree2
coeffs4=polyfit(temperature4,cp4,2); % Polynomial of degree2

predicted_cp1=polyval(coeffs1,temperature1,[],mu1);
predicted_cp2=polyval(coeffs2,temperature2,[],mu2);
predicted_cp3=polyval(coeffs3,temperature3);
predicted_cp4=polyval(coeffs4,temperature4);


% Creating predicted_cp vector for error calculations
predicted_cp=[predicted_cp1;predicted_cp2;predicted_cp3;predicted_cp4];

% Compare my curve fit with the original data

%Plot
% original data curve
plot(temperature,cp,'linewidth',2,'color','b')
hold on

% predicted data curve or curve_fit
plot(temperature1,predicted_cp1,'linewidth',2,'color','r')
plot(temperature2,predicted_cp2,'linewidth',2,'color','r')
plot(temperature3,predicted_cp3,'linewidth',2,'color','r')
plot(temperature4,predicted_cp4,'linewidth',2,'color','r')

ylabel('Specific Heat [KJ/Kmol-K]')
xlabel('Temperature [K]')
legend('original data set','curvefit')
%error calculations
error=(cp-predicted_cp).*100./cp;
fprintf('Maximum error percentage =');
disp(max(error));

Splitwise Curve Fit Data-Cubic Curve Fit (3rd order):