Problems using Linear Regression and syntax
4 次查看(过去 30 天)
显示 更早的评论
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
2 个评论
John D'Errico
2023-2-24
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.
回答(2 个)
Askic V
2023-2-24
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
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?
Torsten
2023-2-24
编辑:Torsten
2023-2-25
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)
A = sol(1)
B = sol(2)
C = sol(3)
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
%********************************************
5 个评论
Torsten
2023-2-25
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Linear and Nonlinear Regression 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!