Problems using Linear Regression and syntax

Can I get some help getting this code to run?

I am trying to use linear regression and the golden search method to find the constants in the Antoine equation. I think my logic is correct -- I just have issues with MATLAB's syntax. Any help is appreciated.

*************************************

function sumSquares = residuals(T,P,C)
Y = log(P);
X = [ones(size(T)), 1./(T+C)]; 
AB = X\Y; % complete the code
A = AB(1);
B = AB(2);
sumSquares = sum((log(P) - A - (B./(T+C))).^2);
end
********************************************
function ABC = fitAntoine(T,P)
T = [273.15; 300.00;   310.00;   320.00;   330.00;  340.00;  350.00;    360.00;    370.00;    380.00;    390.00;   400.00];
P = [728.56; 3911.92; 6820.02; 11205.84; 18049.22; 28158.07; 42655.42; 63062.41; 91129.51; 129027.61; 179180.44; 244541.42];
Ca = -10;
Cb = 10;
Cbeta = Cb - 0.8*(Cb-Ca);
Calpha = Ca + 0.8 * (Cb-Ca);
fCalpha = residuals(T, P, Calpha);
fCbeta = residuals(T, P, Cbeta);
while abs(fCalpha-fCbeta) >sqrt(eps('double'))
    
    Cbeta = Cb - 0.8*(Cb-Ca);
    Calpha = Ca + 0.8 * (Cb-Ca);
    fCalpha = residuals(T, P, Calpha);
    
    fCbeta = residuals(T, P, Cbeta);
    
    if fCbeta >= fCalpha
        Ca = Cbeta;
        C = min(Ca, Cb);
    else 
        Cb = Calpha;
        C = min(Ca, Cb);
    end
    
end
Y = log(P);
X = [ones(size(T)), 1./(T+C)]; 
AB = X\Y; % complete the code
A = AB(1);
B = AB(2);
ABC = [A, B, C]; % complete the code
end

[John D'Errico]

You are not doing something overtly bad, but it is a poor idea to write your own code to do that which is already done better using existing tools. For example, use polyfit to perform the linear regression, combined with a tool like fminbnd to minimize the sum of squares, and therefore to solve for C.

Essentially, you could write much of what you did here, in just a few short lines.

[Dylan]

I appreciate the suggegtions, Ill try fiddling around with those functions. Currently I'm a student so I'm working off the professors example, and many of the students have less experience with MatLAB than myself, so I'm sure he wants us to understand the workflow on a more fundamental level. As it stands though I cannot get my output to match the professors. That being

ABC =

1.0e+03 *

0.0236

-4.1789

-0.0281

Answer 1

[Askic V]

Can you clarify what exactly is the problem with Matlab syntax?

Here is what is obtained if your code is executed:

clear
clc
T = [273.15; 300.00;   310.00;   320.00;   330.00;  340.00;  350.00;    360.00;    370.00;    380.00;    390.00;   400.00];
P = [728.56; 3911.92; 6820.02; 11205.84; 18049.22; 28158.07; 42655.42; 63062.41; 91129.51; 129027.61; 179180.44; 244541.42];
Ca = -10;
Cb = 10;
Cbeta = Cb - 0.8*(Cb-Ca);
Calpha = Ca + 0.8 * (Cb-Ca);
fCalpha = residuals(T, P, Calpha);
fCbeta = residuals(T, P, Cbeta);
    while abs(fCalpha-fCbeta) >sqrt(eps('double'))
  
        Cbeta = Cb - 0.8*(Cb-Ca);
        Calpha = Ca + 0.8 * (Cb-Ca);
        fCalpha = residuals(T, P, Calpha);
    
        fCbeta = residuals(T, P, Cbeta);
   
        if fCbeta >= fCalpha
            Ca = Cbeta;
            C = min(Ca, Cb);
        else 
            Cb = Calpha;
            C = min(Ca, Cb);
        end
        
    end
    Y = log(P);
    X = [ones(size(T)), 1./(T+C)]; 
    AB = X\Y; % complete the code
    A = AB(1);
    B = AB(2);
    ABC = [A, B, C] % complete the code

ABC = 1×3

1.0e+03 * 0.0245 -4.6973 -0.0100

    
    p1 = exp(A+B./(C+T)); % calculate new p values
    
    plot(T,p1)
    hold on
    plot(T,P,'o')
    
    legend('Calculated coefficients', 'Original data')

    
function sumSquares = residuals(T,P,C)
    Y = log(P);
    X = [ones(size(T)), 1./(T+C)]; 
    AB = X\Y; % complete the code
    A = AB(1);
    B = AB(2);
    sumSquares = sum((log(P) - A - (B./(T+C))).^2);
end

What issues are you experiencing?

[Dylan]

Well it looks like C is just approaching the bounds initially set by Cb whenever I change it. The syntax issues I'm having are related to the examples I'm using to build this code. Specifically naming conventions when dont you have to have some sort of header at the top of the code. Like what is the difference b/t

function [m b] = regressionPractice(x, y)

and

function sumsquares = residualsPractice(T, P, C)

Also why werent the variables I had saved in the workspace not showing up when I entered the loop?

I am defiitely missing some base knowledge that resricts me from fully utilizing MatLAB. Any good resources to get practice with that kind of stuff.

This is the professors example output.

ABC =

1.0e+03 *

0.0236

-4.1789

-0.0281

Answer 2

[Torsten]

For the Antoine equation, you usually work with log10:

T = [273.15; 300.00;   310.00;   320.00;   330.00;  340.00;  350.00;    360.00;    370.00;    380.00;    390.00;   400.00];
P = [728.56; 3911.92; 6820.02; 11205.84; 18049.22; 28158.07; 42655.42; 63062.41; 91129.51; 129027.61; 179180.44; 244541.42];
Ca = -10;
Cb = 10;
Cbeta = Cb - 0.8*(Cb-Ca);
Calpha = Ca + 0.8 * (Cb-Ca);
fCalpha = residuals(T, P, Calpha);
fCbeta = residuals(T, P, Cbeta);
while abs(fCalpha-fCbeta) >sqrt(eps('double'))
    
    Cbeta = Cb - 0.8*(Cb-Ca);
    Calpha = Ca + 0.8 * (Cb-Ca);
    fCalpha = residuals(T, P, Calpha);
    
    fCbeta = residuals(T, P, Cbeta);
    
    if fCbeta >= fCalpha
        Ca = Cbeta;
        C = min(Ca, Cb);
    else 
        Cb = Calpha;
        C = min(Ca, Cb);
    end
    
end
Y = log10(P);
X = [ones(size(T)), 1./(T+C)]; 
AB = X\Y; % complete the code
A = AB(1);
B = AB(2);
% Do nonlinear regression
fun = @(x,xdata) 10.^(x(1)+x(2)./(xdata+x(3)));
sol = lsqcurvefit(fun,[A B C],T,P)

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

sol = 1×3

1.0e+03 * 0.0103 -1.8404 -0.0258

A = sol(1)

A = 10.3063

B = sol(2)

B = -1.8404e+03

C = sol(3)

C = -25.7850

hold on
plot(T,P,'o')
plot(T,10.^(A+B./(T+C)))
hold off
grid on

%*************************************
function sumSquares = residuals(T,P,C)
Y = log10(P);
X = [ones(size(T)), 1./(T+C)]; 
AB = X\Y; % complete the code
A = AB(1);
B = AB(2);
sumSquares = sum((log10(P) - A - (B./(T+C))).^2);
end
%********************************************

[Dylan]

This is the professors example output.

ABC =

1.0e+03 *

0.0236

-4.1789

-0.0281

So I believe I am making a mistake somewhere

[Torsten]

Here is the result that you should get:

T = [273.15; 300.00;   310.00;   320.00;   330.00;  340.00;  350.00;    360.00;    370.00;    380.00;    390.00;   400.00];
P = [728.56; 3911.92; 6820.02; 11205.84; 18049.22; 28158.07; 42655.42; 63062.41; 91129.51; 129027.61; 179180.44; 244541.42];
fun = @(x,xdata) x(1)+x(2)./(xdata+x(3));
A = 23.6;
B = -4180;
C = -28.1;
format long
sol = lsqcurvefit(fun,[A B C],T,log(P))

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

sol = 1×3

1.0e+03 * 0.023642583870501 -4.179436973254965 -0.028039067448448

And it is identical to what your professor got.

[Torsten]

The usual way to solve the fitting problem:

syms A B C
T = sym('T',[12 1]);
P = sym('P',[12,1]);
f = sum((log(P)-(A+B./(T+C))).^2);
dfdA = diff(f,A);
dfdB = diff(f,B);
dfdC = diff(f,C);
Tnum = [273.15; 300.00;   310.00;   320.00;   330.00;  340.00;  350.00;    360.00;    370.00;    380.00;    390.00;   400.00];
Pnum = [728.56; 3911.92; 6820.02; 11205.84; 18049.22; 28158.07; 42655.42; 63062.41; 91129.51; 129027.61; 179180.44; 244541.42];
%sol = vpasolve([subs(dfdA,[T P],[Tnum,Pnum])==0,subs(dfdB,[T P],[Tnum,Pnum])==0,subs(dfdC,[T P],[Tnum,Pnum])==0],[A B C])
dfdA_num = subs(dfdA,[T P],[Tnum,Pnum]);
dfdB_num = subs(dfdB,[T P],[Tnum,Pnum]);
dfdC_num = subs(dfdC,[T P],[Tnum,Pnum]);
fun = matlabFunction([dfdA_num,dfdB_num,dfdC_num]);
format long
sol = fsolve(@(x)fun(x(1),x(2),x(3)),[23.6 -4180 -28.1])

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

sol = 1×3

1.0e+03 * 0.023642583992292 -4.179437047451260 -0.028039064786185

[Askic V]

@Torsten Can you please why it is important to include line:

format long

before calling the solve function?

[Torsten]

Can you please why it is important to include line:

format long

before calling the solve function?

Not important - I just wanted to see the solution with more digits to compare to the solution of the OP's professor.