I am getting the following error when I use fsolve "Equation solved, fsolve stalled....."

6 ビュー (過去 30 日間)
Edwin Williams
Edwin Williams 2018 年 12 月 2 日
コメント済み: Alan Weiss 2018 年 12 月 3 日
Hello,
I am ChemE student and I am trying to use fsolve to solve something called an SRK equation of state. I don't know why this is happening. It does say that the equation is solved but it gives the wrong answer. The answer should be 430.97. I double checked my code to make sure the wrong answer does not have another cause.
This is what appears in my command window:
>> xSol = fsolve(@(T) SRK(T),435)
Equation solved, fsolve stalled.
fsolve stopped because the relative size of the current step is less than the
default value of the step size tolerance squared and the vector of function values
is near zero as measured by the default value of the function tolerance.
<stopping criteria details>
xSol =
410.8556
stopping criteria details:
fsolve stopped because the relative norm of the current step, 7.716546e-17, is less than
max(options.StepTolerance^2,eps) = 1.000000e-12. The sum of squared function values,
r = 1.431486e-19, is less than sqrt(options.FunctionTolerance) = 1.000000e-03.
Optimization Metric Options
relative norm(step) = 7.72e-17 max(StepTolerance^2,eps) = 1e-12 (default)
r = 1.43e-19 sqrt(FunctionTolerance) = 1.0e-03 (default)
function y = SRK(T)
% Absolute Pressure in atm %
p = 51.0;
% Absolute pressure converted to Pascals %
P = p*101325;
% Volume in Liters%
v = 35.0;
% Volume converted to cubic meters%
V = v/1000;
% Number of moles %
n = 50.0;
% Molar Volume %
Vh = V/n;
% Gas Constant %
R = 8.314;
% Critical Temperature for Carbon Dioxide in Kelvin%
T1 = 304.2;
% Critical Pressure for Carbon Dioxide in atm%
p1 = 72.9;
% Critical Pressure for Carbon Dioxide converted to pascal%
P1 = p1*101325;
% Critical Temperature for Argon in Kelvin%
T2 = 150.87;
% Critical Pressure for Argon in atm%
p2 = 48.34;
% Critical Pressure for Argon converted to pascal%
P2 = p2*101325;
% Pitzer Acentric Factor for Carbon Dioxide %
w1 = 0.225;
% Pitzer Acentric Factor for Argon %
w2 = -0.002;
% a variable for SRK equation %
a1 = 0.42747*((R*T1)^2)/P1;
a2 = 0.42747*((R*T2)^2)/P2;
% b variable for SRK equation %
b1 = 0.08664*(R*T1/P1);
b2 = 0.08664*(R*T2/P2);
% m variable for SRK equation %
m1 = 0.48508 + 1.55171*w1 - 0.1561*(w1)^2;
m2 = 0.48508 + 1.55171*w2 - 0.1561*(w2)^2;
% Reduced Temperature for Carbon Dioxide %
Tr1 = T/T1;
% Reduced Temperature for Argon %
Tr2 = T/T2;
% alpha variable for SRK equation %
alpha1 = (1 + m1*(1-sqrt(Tr1)))^2;
alpha2 = (1 + m2*(1-sqrt(Tr2)))^2;
% SRK Equation %
y = P - R*T/(Vh - b1) - alpha1*a1/(Vh*(Vh + b1));
y = P - R*T/(Vh - b2) - alpha2*a2/(Vh*(Vh + b2));
end

回答 (1 件)

Alan Weiss
Alan Weiss 2018 年 12 月 3 日
編集済み: Alan Weiss 2018 年 12 月 3 日
Well, fsolve arrived at a solution. If you plug the value that it returns, you see that it is a solution of the equation.
If you plug in the value 430.97 does it give a solution as well?
For a 1-D problem such as this, I would use fzero instead of fsolve, and if I wanted only the solution near 430.97 I would give an initial interval such as [430 432] as the x0 argument for fzero.
Alan Weiss
MATLAB mathematical toolbox documentation
  1 件のコメント
Alan Weiss
Alan Weiss 2018 年 12 月 3 日
You might want to try this experiment:
pp = linspace(430,432);
yy = zeros(size(pp));
for ii = 1:length(pp)
yy(ii) = SRK(pp(ii));
end
plot(pp,yy)
I get a plot showing that the root is not in the interval [430,432].
Alan Weiss
MATLAB mathematical toolbox documentation

サインインしてコメントする。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by