CURVE FITTING USING MATLAB

OBJECTIVES:

1. To write codes for linear and cubic polynomial fit for given `C_p` data.

2. To plot the linear and cubic fit curves along with the raw data points.

3. To write a code to show splitwise method.

4. To measure the fitness characteristics for both the curves.

CURVE FITTING:

Curve fitting is one of the most powerful and most widely used analysis tools. Curve fitting examines the relationship between one or more predictors (independent variables) and a response variable (dependent variable), with the goal of defining a "best fit" model of the relationship. Curve fitting can involve either interpolation, where an exact fit to the data is required, or smoothing, in which a smooth function is constructed that approximately fits the data. Fitted curves can be used as an aid for data visualization, to infer values of a function where no data are available, and to summarize the relationships among two or more variables. 

Perfect fit: Perfect fit is when the curve goes through every point. The sum of squared errors `(SSE)` is 0.0 and `R^2` is 1.00.

Best fit: Best fit simply means that the differences between the actual measured Y values and the Y values predicted by the model equation are minimized. It does not mean a perfect fit.

A perfect fit is always a best fit whereas the best fit is not always a perfect fit.

Linear fit:      `y = ax + b`

Cubic fit:       `y = ax^3 + bx^2 + c`

Polynomial of degree 8 fit:         `y = ax^8 + bx^7 + cx^6 + dx^5 + ex^4 + fx^3 + gx^2 + hx + i`

PARAMETERS TO MEASURE FITNESS OF THE CURVES:

Assume data points `(x_1,y_1)`,`(x_2,y_2)`, .....`(x_n,y_n)` and assume the best fit for these set of points wil be `y = f(x)`.

The parameters which help us measure the fitness of the curves are given below

1. The Sum of Squared Errors `(SSE)`:

Error = | Given value of `y` at `i` - fitted value of `f(x)` at `i` |

Error `= | y(i) - f(x(i)) |`

Sum of Squared Errors `(SSE)` `= sum_(i=1) ^n`Error`(i)^2`

2. R - square `(R^2)`:

This quantity measures how successful the fit is in explaining the variation of the data. R - square is the square of the correlation between the response values and the predicted response values.

The value of R - square ranges from 0 to 1.

• 0 for the worst possible fit for the data
• 1 for thee best possible fit for the data

The value of R - square is caluculated as follows

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

if `y = f(x)` is the fit curve then Sum of Squares of Regression `(SSR)` is given as follows

`SSR = sum_(i=1)^n (f(x(i)) - mu)^2`

where `mu` is the mean of the given data set points.

`SST = SSR + SSE`

3. Adjusted R - square `(R^2)`:

Adjusted R - square is a normalizing variation of R - square which neglects the variation of R - square due to insignificant independent variables.

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

where,

n - number of data points

k - number of independent variables

4. Root Mean Squared Error `(RMSE)`:

It is also known as standard error of the regression. It is an estimate of the standard deviation of the random component in the data. It is given by

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

where,

n - total number of the data points

When `R^2` is greater than 0.95, it is assumed to be good fit.

CODE EXPALNATION:

1. Data from the file 'data' is read into cp_data in which first column is temperature and second is Specific Heat 

2.  Linear fit:

i. The data is fitted using linear curve by using 'polyfit' command and the values of coefficients is stored into 'linear_co_effs'

linear_co_effs = polyfit(temperature, cp, 1)

ii. Using these coefficients, fitted values of specific heat are calculated using 'polyval' command and stored in to 'linear_predicted_cp'.

iii. These values are plotted using plot command.

iv. Whether the fit is good or not is checked by calculating `R^2` value using function Rsquare(cp, linear_predicted_cp)

3. Cubic fit: 

i. The data is fitted using cubic curve by using 'polyfit' command and the values of coefficients is stored into 'cubic_co_effs'

[cubic_co_effs, ~, mu] = polyfit(temperature, cp, 3);

ii. Using these coefficients, fitted values of specific heat are calculated using 'polyval' command and stored in to 'cubic_predicted_cp'

cubic_predicted_cp = polyval(cubic_co_effs, temperature, [], mu);

iii. These values are plotted using plot command.

iv. Whether the fit is good or not is checked by calculating `R^2` value using function Rsquare(cp,cubic_predicted_cp)

4. Splitwise Cubic Curve Fit:

i. Given data is split into 4 different ranges. Each range data is fitted using cubic curve by using 'polyfit' command and the values of coefficients is stored into 'co_effs_n' where n is the number of range.

[co_effs_1, ~, mu] = polyfit(temp_1, cp_1, 3);

ii. Using these coefficients, fitted values of specific heat are calculated using 'polyval' command and stored in to 'predicted_cp_n'

predicted_cp_1 = polyval(co_effs_1, temp_1, [], mu);

iii. These values are plotted using plot command for all ranges.

iv. Whether the fit is good or not is checked by calculating `R^2` value using function Rsquare_split finction

R_square = Rsquare_split(cp_1,predicted_cp_1,cp_2,predicted_cp_2,cp_3,predicted_cp_3,cp_4,predicted_cp_4)

5. Polynomial of Degree 8 Curve Fit:

i. The data is fitted using polynomial of degree 8 curve by using 'polyfit' command and the values of coefficients is stored into 'co_effs'

[co_effs, ~, mu] = polyfit(temperature, cp, 8);

ii. Using these coefficients, fitted values of specific heat are calculated using 'polyval' command and stored in to 'predicted_cp'

predicted_cp = polyval(co_effs, temperature, [], mu);

iii. These values are plotted using plot command.

iv. Whether the fit is good or not is checked by calculating `R^2` value using function Rsquare(cp,predicted_cp)

syntax for 'polyfit' command:

`p(x) = p_1x^n + p_2x^(n-1) + .... + p_n x + p_(n+1)`

p = polyfit(x,y,n) returns the coefficients for a polynomial p(x) of degree n that is a best fit (in a least-squares sense) for the data in y.
The coefficients in p are in descending powers, and the length of p is n+1

[p,S] = polyfit(x,y,n) also returns a structure S that can be used as an input to polyval to obtain error estimates.

[p,S,mu] = polyfit(x,y,n) also returns mu, which is a two-element vector with centering and scaling values. mu(1) is mean(x), and mu(2) is standard deviation std(x). Using these values, polyfit centers x at zero and scales it to have unit standard deviation,

`hat x = (x - bar x)/ (sigma_x)`

This centering and scaling transformation improves the numerical properties of both the polynomial and the fitting algorithm.

syntax for 'polyval' command:

CODE:

% CURVE FITTING 
clear all
close all 
clc

cp_data = load('data');         % Loading data from the file 'data'
temperature = cp_data(:,1);     % Getting temperature data from column 1 of cp_data 
cp = cp_data(:,2);              % Getting specific heat data from column 2 of cp_data 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LINEAR FIT FOR CP DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cp = aT + b

% Linear curve fit
% a & b values are obtained for linear polynomial (cp = aT + b) using 'polyfit'
linear_co_effs = polyfit(temperature, cp, 1);       
% At each value of 'temperature', predicted vaule of 'cp' is calculated using 'polyval'
linear_predicted_cp = polyval(linear_co_effs, temperature); 

% Plotting
figure(1)
% Original plot of T vs Cp data from cp_data file
plot(temperature, cp, 'linewidth', 3, 'color', 'b');
hold on 
% Linear curve fit plot of T and predicted cp
plot(temperature, linear_predicted_cp, 'linewidth', 3, 'color', 'r')
hold off
xlabel('Temperature [K]')
ylabel('Specific Heat [KJ/Kmol-K]');
axis([0 3600 900 1350]);
legend({'Original Data','Linear Curve fit'},'Position',[0.3 0.75 0.1 0.1]);
title('C_p vs T - LINEAR CURVE FITTING')

% Caluculation of R_square value for linear fit curve
R_square_linear = Rsquare(cp,linear_predicted_cp) 

if R_square_linear >= 0.95
    disp('LINEAR CURVE FIT IS GOOD')
else
    disp('LINEAR CURVE FIT IS NOT GOOD')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CUBIC FIT FOR CP DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cp = a(T^3) + b(T^2) + cT + d
% a,b,c,d values are obtained for linear polynomial (cp = a(T^3) + b(T^2) + cT + d) using 'polyfit'
% 'polyfit' with three outputs is used to fit a 3rd-degree polynomial using centering and scaling,
% which improves the numerical properties of the problem
% 'polyfit' centers the data in 'temperature' at 0 and scales it to have a standard deviation of 1, 
% which avoids an ill-conditioned Vandermonde matrix in the fit calculation.

[cubic_co_effs, ~, mu] = polyfit(temperature, cp, 3);

% At each value of 'temperature', predicted vaule of 'cp' is calculated using 'polyval'
% Use 'polyval' with four inputs to evaluate cp with the scaled temperatures, (temperature-mu(1))/mu(2)

cubic_predicted_cp = polyval(cubic_co_effs, temperature, [], mu);

figure(2)
% Original plot of T vs Cp data from cp_data file
plot(temperature, cp, 'linewidth', 3, 'color', 'b');
hold on 
% Cubic curve fit plot of T and predicted cp
plot(temperature, cubic_predicted_cp, 'linewidth', 3, 'color', 'r')
hold off
xlabel('Temperature [K]')
ylabel('Specific Heat [KJ/Kmol-K]');
axis([0 3600 900 1350]);
legend({'Original Data','Cubic Curve fit'},'Position',[0.3 0.75 0.1 0.1]);
title('C_p vs T - CUBIC CURVE FITTING')

% Caluculation of R_square value for cubic fit curve
R_square_cubic = Rsquare(cp,cubic_predicted_cp) 

if R_square_cubic >= 0.95
    disp('CUBIC CURVE FIT IS GOOD')
else
    disp('CUBIC CURVE FIT IS NOT GOOD')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SPLIT-WISE CUBIC FIT FOR CP DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cp = aT^3 + bT^2 + cT + d

temp_1 = cp_data(1:500,1);          % Temperature values for 1st split
temp_2 = cp_data(501:1200,1);       % Temperature values for 2nd split
temp_3 = cp_data(1201:2400,1);      % Temperature values for 3rd split
temp_4 = cp_data(2401:3200,1);      % Temperature values for 4th split

cp_1 = cp_data(1:500,2);            % Specific Heat values for 1st split
cp_2 = cp_data(501:1200,2);         % Specific Heat values for 2nd split
cp_3 = cp_data(1201:2400,2);        % Specific Heat values for 3rd split
cp_4 = cp_data(2401:3200,2);        % Specific Heat values for 4th split

% cp = a(T^3) + b(T^2) + cT + d
% a,b,c,d values are obtained for linear polynomial (cp = a(T^3) + b(T^2) + cT + d) using 'polyfit'
% 'polyfit' with three outputs is used to fit a 3rd-degree polynomial using centering and scaling, 
% which improves the numerical properties of the problem
% 'polyfit' centers the data in 'temperature' at 0 and scales it to have a standard deviation of 1,
% which avoids an ill-conditioned Vandermonde matrix in the fit calculation.

% At each value of 'temperature', predicted vaule of 'cp' is calculated using 'polyval'
% Use 'polyval' with four inputs to evaluate cp with the scaled temperatures, (temperature-mu(1))/mu(2)

% 1st split
[co_effs_1, ~, mu] = polyfit(temp_1, cp_1, 3);          % cp = a(T1^3) + b(T1^2) + cT1 + d
predicted_cp_1 = polyval(co_effs_1, temp_1, [], mu);
% 2nd split
[co_effs_2, ~, mu] = polyfit(temp_2, cp_2, 3);          % cp = a(T2^3) + b(T2^2) + cT2 + d
predicted_cp_2 = polyval(co_effs_2, temp_2, [], mu);
% 3rd split
[co_effs_3, ~, mu] = polyfit(temp_3, cp_3, 3);          % cp = a(T3^3) + b(T3^2) + cT3 + d
predicted_cp_3 = polyval(co_effs_3, temp_3, [], mu);
% 4th split
[co_effs_4, ~, mu] = polyfit(temp_4, cp_4, 3);          % cp = a(T4^3) + b(T4^2) + cT4 + d
predicted_cp_4 = polyval(co_effs_4, temp_4, [], mu);

% Plotting
figure(3)
% Original plot of T vs Cp data from cp_data file
plot(temperature, cp, 'linewidth', 3, 'color', 'b');
hold on 
% Split cubic curve fit of T and predicted cp in 1st split
plot(temp_1, predicted_cp_1, 'linewidth', 3, 'color', 'r')
% Split cubic curve fit of T and predicted cp in 2nd split
plot(temp_2, predicted_cp_2, 'linewidth', 3, 'color', 'r')
% Split cubic curve fit of T and predicted cp in 3rd split
plot(temp_3, predicted_cp_3, 'linewidth', 3, 'color', 'r')
% Split cubic curve fit of T and predicted cp in 4th split
plot(temp_4, predicted_cp_4, 'linewidth', 3, 'color', 'r')
hold off
xlabel('Temperature [K]')
ylabel('Specific Heat [KJ/Kmol-K]');
axis([0 3600 900 1350]);
legend({'Original Data','Splitwise Cubic Curve fit'},'Position',[0.3 0.75 0.1 0.1]);
title('C_p vs T - CUBIC CURVE FITTING WITH 4 SPLITS')

% Caluculation of R_square value for split wise cubic fit curve
R_square_split = Rsquare_split(cp_1,predicted_cp_1,cp_2,predicted_cp_2,cp_3,predicted_cp_3,cp_4,predicted_cp_4) 

if R_square_split >= 0.95
    disp('SPLIT-WISE CUBIC CURVE FIT IS GOOD')
else
    disp('SPLIT-WISE CUBIC CURVE FIT IS NOT GOOD')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% POLYNOMIAL OF DEGREE 8 CURVE FIT FOR CP DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cp = a(T^8) + b(T^7) + c(T^6) + d(T^5) + e(T^4) + f(T^3) + g(T^2) + hT + i
% a,b,c,d,e,f,g,h & i values are obtained using 'polyfit'
[co_effs, ~, mu] = polyfit(temperature, cp, 8);
predicted_cp = polyval(co_effs, temperature, [], mu);

figure(4)
% Original plot of T vs Cp data from cp_data file
plot(temperature, cp, 'linewidth', 3, 'color', 'b');
hold on 
% Polynomial (degree = 8) Curve fit plot of T and predicted cp
plot(temperature, predicted_cp, 'linewidth', 3, 'color', 'r')
hold off
xlabel('Temperature [K]')
ylabel('Specific Heat [KJ/Kmol-K]');
axis([0 3600 900 1350]);
legend({'Original Data','Polynomial of Degree 8 Curve fit'},'Position',[0.3 0.75 0.1 0.1]);
title('C_p vs T - CURVE FITTING OF POLYNOMIAL DEGREE 8')

% Caluculation of R_square value for cubic fit curve
R_square_degree8 = Rsquare(cp,predicted_cp) 

if R_square_degree8 >= 0.95
    disp('POLYNOMIAL OF DEGREE 8 CURVE FIT IS GOOD')
else
    disp('POLYNOMIAL OF DEGREE 8 CURVE FIT IS NOT GOOD')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% R-SQUARE FUNCTION FOR CURVE WITHOUT SPLITTING %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function R_Square = Rsquare(originalcp,predictedcp)
    % Calculation of mean specific heat value of given data
    Total_cp =0;
    for i=1:length(originalcp)
        Total_cp = Total_cp + originalcp(i);
    end

    mean_cp = Total_cp/length(originalcp);
    % Calculation of Sum of Squares of Regresion (SSR)
    SSR = 0;
    for j=1:length(originalcp)
        SSR = SSR + (predictedcp(j) - mean_cp)^2;
    end
    % Calculation of SUm of Squares due to Error (SSE)
    SSE = 0;
    for k=1:length(originalcp)
        SSE = SSE + (originalcp(k) - predictedcp(k))^2;
    end

    SST = SSR + SSE;
    R_Square = SSR/SST;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% R-SQUARE FUNCTION FOR CURVE WITH 4 SPLITS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function R_square = Rsquare_split(originalcp1,predictcp1,originalcp2,predictcp2,originalcp3,predictcp3,originalcp4,predictcp4)
    
    Total_cp_1 =0;
    Total_cp_2 =0;
    Total_cp_3 =0;
    Total_cp_4 =0;
    
    % Calculation of mean specific heat value of given data in 1st split
    for i=1:length(originalcp1)
        Total_cp_1 = Total_cp_1 + originalcp1(i);
    end
    mean_cp_1 = Total_cp_1/length(originalcp1);
    
    % Calculation of mean specific heat value of given data in 2nd split
    Total_cp_2 =0;
    for i=1:length(originalcp2)
        Total_cp_2 = Total_cp_2 + originalcp2(i);
    end
    mean_cp_2 = Total_cp_2/length(originalcp2);
    
    % Calculation of mean specific heat value of given data in 3rd split
    Total_cp_3 =0;
    for i=1:length(originalcp3)
        Total_cp_3 = Total_cp_3 + originalcp3(i);
    end
    mean_cp_3 = Total_cp_3/length(originalcp3);
    
    % Calculation of mean specific heat value of given data in 4th split
    Total_cp_4 =0;
    for i=1:length(originalcp4)
        Total_cp_4 = Total_cp_4 + originalcp4(i);
    end
    mean_cp_4 = Total_cp_4/length(originalcp4);

    % Calculation of Sum of Squares of Regresion (SSR)
    % SSR of all splits are added into one. Hence after each loop, SSR is not set to zero.
    % Hence in the loop, all values are added to previous SSR value
    SSR = 0;
    for j=1:length(originalcp1)
        SSR = SSR + (predictcp1(j) - mean_cp_1)^2;
    end
    for j=1:length(originalcp2)
        SSR = SSR + (predictcp2(j) - mean_cp_2)^2;
    end
    for j=1:length(originalcp3)
        SSR = SSR + (predictcp3(j) - mean_cp_3)^2;
    end
    for j=1:length(originalcp4)
        SSR = SSR + (predictcp4(j) - mean_cp_4)^2;
    end

    % Calculation of Sum of Squares due to Error (SSE)
    % SSE of all splits are added into one. Hence after each loop, SSE is not set to zero.
    % Hence in the loop, all values are added to previous SSE value
    SSE = 0;
    for k=1:length(originalcp1)
        SSE = SSE + (originalcp1(k) - predictcp1(k))^2;
    end
    for k=1:length(originalcp2)
        SSE = SSE + (originalcp2(k) - predictcp2(k))^2;
    end
    for k=1:length(originalcp3)
        SSE = SSE + (originalcp3(k) - predictcp3(k))^2;
    end
    for k=1:length(originalcp4)
        SSE = SSE + (originalcp4(k) - predictcp4(k))^2;
    end
    
    SST = SSR + SSE;
    R_square = SSR/SST;
end

OUTPUT:

The value of `R^2` for linear fit is 0.9249 which is less than 0.95. Hence the linear fit for the given data is not good.

The value of `R^2` for cubic fit is 0.9967 which is greater than 0.95, hence it is a better fit for the given data.

There are two methods which can make the cubic fit in to the best fit. They are 

1. Increasing the degree of polynomial while curve fitting the data 

The values of `R^2` for polynomial fit of different degrees is given as follows

Degree             `R^2` Value

   1                        0.9249
   2                        0.9812
   3                        0.9967
   4                        0.9996
   5                        0.9996
   6                        0.9997 
   7                        0.9999
   8                        1.0000

It can be observed that the values of `R^2` reaches 1 when the degree of polynomial for curve fit is increased to 8. It is the best fit for the given data.

2. Split-Wise Method:

By using Split-Wise Method, curve fit for a particular order of polynomial is improved. Here, the data is split into small ranges, and seperate curves are fitted in each range.

Here, the given data is split into 4 ranges and cubic curve fit is used between these ranges

1st Range : 1 to 500

2nd Range : 501 to 1200

3rd Range : 1201 to 2400

4th Range : 2401 to 3200

The value of R^2 is 1 for this fit.

ERRORS:

Function definition is not supported in this context. Create functions in the code file

While writing function for `R^2`, the input variable 'cp' is used. But 'cp' is already used to store specific heat values in the program. Hence this error was showing.

It was corrected by replacing the 'cp' in function with 'originalcp' 

REFERENCES:

1. https://en.wikipedia.org/wiki/Curve_fitting

2. https://skill-lync.com/projects/Curve-Fitting-and-criterion-to-choose-best-fit-41860

3. https://itectec.com/matlab-ref/matlab-function-polyfit-polynomial-curve-fitting/

4. https://terpconnect.umd.edu/~toh/models/CurveFitting.html

5. https://skill-lync.com/projects/program-in-matlab-for-curve-fitting