PERFORMING CURVE FITTING USING MATLAB

PROJECT NAME – PERFORMING CURVE FITTING USING MATLAB

AIM-

1. Learn curve fitting using MATLAB and can observe how data changes mathematically.
2. Evaluation of goodness of fit and error criteria.

PROCEDURE-

In this project we are going to learn curve fitting using MATLAB. Curve fits are really useful to describe how your data changes mathematically & once you have curve fit you can actually do additional mathematical operations that might not be possible in the initial data set.

As an e.g. – Take population data in a particular country there is no known formula that can actually tell you population of particular country going to be in the next year. However, if you actually do a curve fit you actually get a mathematical function which says how the population is going to change as a function of time, in other words you put in future time & predict what the population is going to be. Now based upon curve fit accuracy actual values may vary and that variation depends upon your accuracy of curve fit and you will get the rough trends & since curve fit is the mathematical function you can make additional operations like taking derivatives and so on.

Now to run this programme you should have file named data in the same directory where your MATLAB programmes are. This data file contains two columns first column is temperature and second column is of specific heat. And to make sure that data file is present at proper location type pwd command in MATLAB in command window and it will show you the path of the data file.

Now, first thing you will do is to load your data set and the best thing about MATLAB is that it has an inbuilt load command to load files and file we are going to read is data file so we will type that in editor of MATLAB. After plotting temperature and sp. Heat plot will show variation of C_p as a function of temperature and by seeing the graph it is clear that C_p varies quite a bit as a function of temperature.

Since you have plotted curve you also have to fit the curve. As you have the data set but you don’t know the mathematical nature of the data set so that’s why we will use curve fit and the idea behind using curve fit is quite simple you can assume the function either quadratic equation/ cubic equation/ n^th order polynomial it doesn’t matter so we have typed the second order quadratic equation. And your job now is to fix up what is the best value for a & b. so that I can use the particular format to use the data and to determine a & b you need to equations so you can see in your curve you will take two points for e.g. if the temperature is 300K you will know the value of C_p & when the temperature is 1000K you know the value for C_p .you can use these 2 datasets and figure out the value for a & b and this is going to be a very simple curve fit in order to make more robust curve fit one should have to pick up the values for a & b in such a way that error that will you get between your original dataset and fitted dataset is minimum and in MATLAB you can do it by using your polyfit command.

Ployfit(temperature, Cp,2)…… where n =2;

n is the order of polynomials so if you type 2 then it will be a quadratic equation and this will get the set of coefficients a, b & c. as we will assign it to a variable called co_effs and when you will run these commands your coefficients and these are coefficients of polynomials and by using these we can predict the value for C_p.

Now when we observe the new plot after comparison of the curves well, it is trying to match our dataset but it doesn’t do a good job so we can increase the order of polynomial after changing it will try almost to do much better job but in command window you can see the warning regarding polynomial is badly conditioned so again we have to figure out a better way to fit the data and for that-

1. you can try out different polynomial orders or
2. you can use split wise polynomials i.e. one polynomial for one temperature range and also try to come up with the way to measure error in the solution and try to fit in that particular condition.

clear all
close all
clc
 
% Load data
cp_data = load('data');
 
% Preparing the data
temperature = cp_data(:,1);
cp = cp_data(:,2);
% Plot
plot(temperature,cp,'linewidth', 3, 'color','b')

Fig – Nature of specific heat varying with temperature

clear all
close all
clc
 
% Load data
cp_data = load('data');
 
% Preparing the data
temperature = cp_data(:,1);
cp = cp_data(:,2);
 
% Curve fit
% cp = a*T^2 + b*T + c

co_effs = polyfit(temperature,cp,3);

predicted_cp = polyval(co_effs, temperature);

% compare my curve fit with the orginal data
plot(temperature,predicted_cp,'linewidth' , 3, 'color', 'r')
hold on

plot(temperature,cp,'linewidth', 3, 'color','b')
xlabel('temperature [k]')
ylabel('specific heat[KJ/Kmol-k]')
legend('original dataset','curve fit')

Fig – Comparison of curve fit with original data

Code to fit a linear and cubic polynomial for C_p data-

clear all
close all
clc

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

% Curve fit for first degree polynomial
%cp = a*t+b
co_effs1 = polyfit(temperature,cp,1);
predicted_cp1 = polyval(co_effs1,temperature);

% Curve fit for second degree polynomial
%cp = a*t^2 + b*t + c
co_effs2 = polyfit(temperature,cp,2);
predicted_cp2 = polyval(co_effs2,temperature);

% Curve fit for third degree polynomial
%cp = a*t^3+b*t^2+ct + d
co_effs3 = polyfit(temperature,cp,3);
predicted_cp3 = polyval(co_effs3,temperature);
 
% compare my curve fit with first degree polynomial with the original data
figure(1)
plot(temperature,cp,'linewidth',3,'color','b')
hold on
plot(temperature,predicted_cp1,'linewidth',3,'color','r')
xlabel('Temperature[k]')
ylabel('Specific Heat [KJ/Kmol-k]')
legend('original dataset','Temperature vs. predicted cp1')
title('First degree polynomial')
 
% compare my curve fit with second degree polynomial with the original data
figure(2)
plot(temperature,cp,'linewidth',3,'color','b')
hold on
plot(temperature,predicted_cp2,'linewidth',3,'color','r')
xlabel('Temperature[k]')
ylabel('Specific Heat [KJ/Kmol-k]')
legend('original dataset','Temperature vs. predicted_cp2')
title('second degree polynomial')

% compare my curve fit with third degree polynomial with the original data
figure(3)
plot(temperature,cp,'linewidth',3,'color','b')
hold on
plot(temperature,predicted_cp3,'linewidth',3,'color','r')
xlabel('Temperature[k]')
ylabel('Specific Heat [KJ/Kmol-k]')
legend('original dataset','Temperature vs. predicted cp3')
 
%Finding  the SSR, SST, R_sq, R_sq_adj for linear polynomial
Error1 = cp - predicted_cp1;
% Calculating SSE1
sse = 0;
sum1 = 0;
for i = 1:length(cp)
    SSE1 = sse+Error1(i)^2;
    SUM1 = sum1+cp(i);
end

n= length(cp)
k= 2
mean = SUM1/n
 
% calculating SSR1
ssr = 0;
for j = 1:length(predicted_cp1);
    SSR1 = ssr + (predicted_cp1(i)-mean)^2;
end

% calculating SST1
SST1 = SSR1 + SSE1
% Calculating R^2 1
R_sq1 = SSR1/SST1

% Calculating Radj^2 1
R_sq_adj1 = 1 - ((1-R_sq1)*(n-1)/(n-k-1))
RMSE1 = (SSE1/n)^0.5

%Finding  the SSR, SST, R_sq, R_sq_adj for linear polynomial
Error3 = cp - predicted_cp3;


% Calculating SSR3
sse = 0;
sum3 = 0;
for i = 1:length(cp)
SSE3 = sse+Error3(i)^2;
    SUM3 = sum3+cp(i);
end

n= length(cp)
k= 2
mean = SUM3/n
 
% calculating SSR3
ssr = 0;
for j = 1:length(predicted_cp1);
    SSR3 = ssr + (predicted_cp3(i)-mean)^2;
end

% calculating SST3
SST3 = SSR3 + SSE3
% Calculating R^2 3
R_sq3 = SSR3/SST3

% Calculating Radj^2 3
R_sq_adj3 = 1 - ((1-R_sq3)*(n-1)/(n-k-1))
RMSE1 = (SSE3/n)^0.5

Step 1 – As usual, to clear all the previous plots, data files start the process in editor window with commands clear all, close all , clc.

Step 2 – as we discussed above to load the data file in MATLAB in which the first column is temperature and second column is specific heat.

Step 3- here we are using the process of different polynomial orders so in first stage we are going to use the first degree polynomial and same as above process as we discussed at first here also we will use the polyfit and polyval command to fit the data(by least square polynomial method) as well as fitting required polynomial values in the curve and with the help of these two functions we have found out predicted_cp1.

Step 4- now we will use second degree polynomial and with the help of polyfit and polyval command we will find out predicted_cp2.

Step 5 – now we will find out cubic polynomial or 3 degree polynomial here also we will use the polyfit and polyval command to fit the data(by least square polynomial method) as well as fitting required polynomial values in the curve and with the help of these two functions we have found out predicted_cp3.

Step 6 - compare my curve fit with first, second, third polynomial with the original first, second, and third data

figure(1)

plot(temperature,cp,'linewidth',3,'color','b')

hold on

plot(temperature,predicted_cp1,'linewidth',3,'color','r')

xlabel('Temperature[k]')

ylabel('Specific Heat [KJ/Kmol-k]')

legend('original dataset','Temperature vs. predicted cp1')

title('First degree polynomial')

………….

figure(2)

plot(temperature,cp,'linewidth',3,'color','b')

hold on

plot(temperature,predicted_cp2,'linewidth',3,'color','r')

xlabel('Temperature[k]')

ylabel('Specific Heat [KJ/Kmol-k]')

legend('original dataset','Temperature vs. predicted_cp2')

title('second degree polynomial')

% compare my curve fit with third degree polynomial with the original data

……

figure(3)

plot(temperature,cp,'linewidth',3,'color','b')

hold on

plot(temperature,predicted_cp3,'linewidth',3,'color','r')

xlabel('Temperature[k]')

ylabel('Specific Heat [KJ/Kmol-k]')

legend('original dataset','Temperature vs. predicted cp3')

step 7 – now we will find following terms for linear and cubic polynomials

1. The sum of squares due to error.
2. R – Square.
3. Adjusted R – Square.
4. Root mean square error. (RMSE)

In short finding SSE, SSR, SST, R_sq, R_sq_adj, RMSE for linear polynomial.

Error1 = cp - predicted_cp1;

% Calculating SSE1

sse = 0;

sum1 = 0;

for i = 1:length(cp)

SSE1 = sse+Error1(i)^2;

SUM1 = sum1+cp(i);

end

n= length(cp)

k= 2

mean = SUM1/n

% calculating SSR1

ssr = 0;

for j = 1:length(predicted_cp1);

SSR1 = ssr + (predicted_cp1(i)-mean)^2;

end

% calculating SST1

SST1 = SSR1 + SSE1

% Calculating R^2 1

R_sq1 = SSR1/SST1

% Calculating Radj^2 1

R_sq_adj1 = 1 - ((1-R_sq1)*(n-1)/(n-k-1))

RMSE1 = (SSE1/n)^0.5

%Finding the SSR, SST, R_sq, R_sq_adj for cubic polynomial

Error3 = cp - predicted_cp3;

% Calculating SSR3

sse = 0;

sum3 = 0;

for i = 1:length(cp)

SSE3 = sse+Error3(i)^2;

SUM3 = sum3+cp(i);

end

n= length(cp)

k= 2

mean = SUM3/n

% calculating SSR3

ssr = 0;

for j = 1:length(predicted_cp1);

SSR3 = ssr + (predicted_cp3(i)-mean)^2;

end

% calculating SST3

SST3 = SSR3 + SSE3

% Calculating R^2 3

R_sq3 = SSR3/SST3

% Calculating Radj^2 3

R_sq_adj3 = 1 - ((1-R_sq3)*(n-1)/(n-k-1))

RMSE1 = (SSE3/n)^0.5

Use of curve fitting tool box to find out SSE,SST, R², SSR,RMSE.

SPLITWISE CURVE FITTING -

% Programme to make split wise  curve fit
clear all
close all
clc

% Preparing the data
data = load('data');
temp = data(:,1);
cp = data(:,2);

% Counting the frame
a =1;
% Third order curve fit
% Divide the loop into 32 different parts
for i = 1:200:3200 
    [p,~,mu] = polyfit(temp(i:i +199),cp(i:i+199),3);
    % Getting cp values foreach temperature of third order curve fit
    cp3(i:i + 199) = polyval(p,temp(i:i+199),[],mu);
end

% Loop to animate original curve
for i = 1:100:3200                 
    plot(temp(i:i + 99),cp(i:i + 99),'r','linewidth',2);
    axis([0 3500 900 1350]);
    xlabel('temperature[k]');
    ylabel('specific heat [KJ/Kmol-k]');
    title('curve fitting');
    hold on
    Frame(a) = getframe(gcf)
    a = a + 1;
end

% Curve fit cubic curve
for i = 1:50:3200
     plot(temp(i:i + 49),cp3(i:i + 49),'k:','linewidth',2);
     Frame(a) = getframe(gcf)
    a = a + 1;
end

% Plot original curve
figure(2)

plot(temp,cp,'r','linewidth',2);
hold on
plot(temp,cp3,':k','linewidth',2);
legend('original curve','cubic curve','location','northwest');
% Creating, Opening, Writing & closing file by making its movie
videofile = VideoWriter('Fitted_curve_corrected.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,Frame)
close(videofile)

First of all we prepared the data by loading the original data file which was already present in the MATLAB directory. Then, for better curve fitting we will use splitwise polynomial method by differentiating the temperatures points from temperature data and to eliminate the error we will use centring and scaling method. Till here we have tried using linear/ cubic method, we also used tool to plot the curve fits using the data of linear and cubic curve fits but as we can see we haven’t reached to a point where our curve have least variation with the original data and to compensate this effect now we will use split wise method so that we can capture the best fit overlapping perfectly with the original data curve.

Here, we will use command [p, ~, mu] which denotes as the variable for centre and scaling and hence error estimation will be resolved.

For a given set of polynomial Polyfit computes least squares polynomials and polyval is used to quickly and easily fit polynomial with a set of data points and these both commands are used to find the predicted C_p.

As the method’s name itself describes its name we have splited the temperatures value and specific heat of oxygen’s values into 32 different parts along with spillting temperatures and specific heat values for linear, cubic, original curve fits and also to make the movie from it so that we can see live capturing of the data we also have attached counter as well as frame commands to the codes and saved the file and dumped it into the AVI format.

YOUTUBE LINK -

https://youtu.be/K8yDeSFMCP4

CONCLUSIONS-

1. For good fits higher degree of polynomial is to be selected.
2. There are several ways to with which we can do curve fitting.
3. Also with the help of curve fitting tool box we can find out SSE, SST, R², SSR R__ADJ
4. After using the splitwise method of curve fitting we see no message in command window regarding poorness of curve fit.

1. Parameters used to measure fitness characteristics of both the curves.

following terms for linear and cubic polynomials

1. The sum of squares due to error.
2. R – Square.
3. Adjusted R – Square.
4. Root mean square error. (RMSE)

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

SSE can be defined as squared differences between each observation and its group mean

Assume the above equation as n^th degree polynomial of the form

F(x) = A₁ + A₂ . x +A₃.x² + A_{n .} x^n-1

When we differentiate the SSE with each and every coefficient of the proposed equation f(x) and equate it to zero, we can find minima and value of that coefficient.

Upon solving the multiple system of equations as we can find the values of all coefficients. This method is known as least square method and most commonly used technique as long as fitting the function of a polynomial.

This method is called sum of squares due to error (SSE).

2. R – Square :-

The quantity that measures how successful the fit is in explaining the variation of the data.

R square is the square of the correlation between response values and predicted response values. In layman terms higher the value of R square better the fit and lower the fit lousier the fit. R square ranges from 0 to 1.

SSR = Sum of squares of the regression.

3. Adjusted R square :-

Adjusted R square is a normalising variation of R – Square which neglects the variation of R square due to insignificant independent variables.

4. Root Mean Square Method

This static is also known as the fit standard error and the standard error of the regression. It’s an estimate the standard deviation of random component in the data and is defined as –

n = total no. of datapoints

Using these four quantities MATLAB effectively deduces the goodness fit typically R – square is greater than 0.95 the fit is assumed to be good.

1. How to make curve fit perfect?

First method of getting perfect curvefit is by interpolation by which errors can be minimized and errors can be reduced to zero.

Second method is by centring and scaling method, the curve fit should match the original fit in that particular order of polynomials.

A perfect fit can always be the best fit but the best fit cannot be the perfect fit always.

In short the difference between actual measured value and predicted value should be minimum.

2. How to make the best fit?

Best fit and perfect fit are completely different phenomenon. To get the best fit we can increase the degree of polynomials. To see the effect of polynomial degree see the function in MATLAB as polyval and polyfit with the random data.

Splitwise method is also best way to get the best curve fit. More the number of splits better the accuracy.

The best fit can be defined by four different methods which are explained above –

3. The sum of squares due to error.
4. R – Square.
5. Adjusted R – Square.
6. Root mean square error. (RMSE)

Best fit is the fit in which the value of R square is 1. And this is only possible when actual curve matches with predicted curve perfectly. Sum of the squares is 0.0.

3. What could be done to improve the cubic fit?

While operating with polynomials when cubic curve was giving error of equation is badly conditioned we can use the centering and scaling method.

It can be also done by dividing the temperature into different ranges as explained above by splitwise function.

It can be done by reducing curve fit value and original dataset.

ERRORS FACED - No major errors faced as such.