An issue while applying Newton method for solving nonlinear systems.

22 次查看（过去 30 天）
Fadi Tantish
Fadi Tantish 2021-6-25
评论： Fadi Tantish 2021-6-29
Hello eveyone,
I am facing an issue in applying my code to solve two equations with two variables at a range of different temperatures using Newton's method. My code is as follow:
close all; clc,clear all
% Determine saturated vapor pressure, saturated liquid pressure,&
% their saturation volumes at different temperatures USING NEWTON RAPHSON
% method.
% Condition: " function [f, j] = Psat_fNR(X,T)"
% P_sat == P_l == P_v (f(1,1))&& fugacity_liq == fugacity_vap (f(2,1)).
R=8.3145; % universal gas constant in MPa.cm3/(K.mol)
a = 553661.24; % vdW constant a for water in MPa.cm6/mol2
b = 30.49; % vdW constant b for water in cm3/mol
Tc = 8*a/(27*b*R); % critical temperature of water.
Pc = a/(27*b^2); % critical pressure of water.
T = 400 : 1 : 650; % range of temperatures to be investigated.
maxiter = 200;
tol = 10^-6;
for k = 1 : numel(T)
X0 = [31; 2000];
X = X0;
Xold = X0;
for i = 1:maxiter
[f,j] = Psat_fNR(X,T(k));
X = X - j\f;
err(:,i) = abs(X-Xold);
Xold = X;
if (err(:,i) < tol)
vl = X(1,1);
vv = X(2,1);
P_l = R*T/(vl-b) - a/(vl^2);
P_v = R*T/(vv-b) - a/(vv^2);
break;
end
end
xv(k,:) = X;
end
function [f,j] = Psat_fNR(X,T)
R=8.3145;
a = 553661.24;
b = 30.49;
vl = X(1);
vv = X(2);
f(1,1) = (R*T/(vl-b) - a/(vl^2)) - (R*T/(vv-b) - a/(vv^2));
f(2,1) = ((1/(1-(b/vl))-a/(R*T*vl))-1 - log((1/(1-(b/vl))-a/(R*T*vl))-((b*(R*T/(vl-b) - a/(vl^2)))/(R*T))) - ((a*(R*T/(vl-b) - a/(vl^2)))/(R*T)^2)/(1/(1-(b/vl))-a/(R*T*vl))) - ((1/(1-(b/vv))-a/(R*T*vv))-1 - log((1/(1-(b/vv))-a/(R*T*vv))-((b*(R*T/(vv-b) - a/(vv^2)))/(R*T))) - ((a*(R*T/(vv-b) - a/(vv^2)))/(R*T)^2)/(1/(1-(b/vv))-a/(R*T*vv))) ;
% first derivatives of f(1,1) with respect to vl and vv.
df_vl= (2*a)/vl^3 - (R*T)/(b - vl)^2;
df_vv = (R*T)/(b - vv)^2 - (2*a)/vv^3;
% first derivatives of f(2,1) with respect to vl and vv.
dfug_vl = a/(R*T*vl^2) - b/(vl^2*(b/vl - 1)^2) - (b/(vl^2*(b/vl - 1)^2) + (b*((2*a)/vl^3 - (R*T)/(b - vl)^2))/(R*T) - a/(R*T*vl^2))/(1/(b/vl - 1) - (b*(a/vl^2 + (R*T)/(b - vl)))/(R*T) + a/(R*T*vl)) + (a*((2*a)/vl^3 - (R*T)/(b - vl)^2))/(R^2*T^2*(1/(b/vl - 1) + a/(R*T*vl))) + (a*(b/(vl^2*(b/vl - 1)^2) - a/(R*T*vl^2))*(a/vl^2 + (R*T)/(b - vl)))/(R^2*T^2*(1/(b/vl - 1) + a/(R*T*vl))^2);
dfug_vv = (b/(vv^2*(b/vv - 1)^2) + (b*((2*a)/vv^3 - (R*T)/(b - vv)^2))/(R*T) - a/(R*T*vv^2))/(1/(b/vv - 1) - (b*(a/vv^2 + (R*T)/(b - vv)))/(R*T) + a/(R*T*vv)) + b/(vv^2*(b/vv - 1)^2) - a/(R*T*vv^2) - (a*((2*a)/vv^3 - (R*T)/(b - vv)^2))/(R^2*T^2*(1/(b/vv - 1) + a/(R*T*vv))) - (a*(b/(vv^2*(b/vv - 1)^2) - a/(R*T*vv^2))*(a/vv^2 + (R*T)/(b - vv)))/(R^2*T^2*(1/(b/vv - 1) + a/(R*T*vv))^2);
j = [df_vl, df_vv;
dfug_vl, dfug_vv];
end
I keep receiving a warning message for line 30. I don't know what is the issue!
warning message:
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.
> In NR_Psat_fu (line 30)
line 30 : err(:,i) = abs(X-Xold);
I would be glad if helped me to solve this issue and teach me how to create a table as a final outcome of the code that include T, vl, vv, P_l, and P_v .
Many thanks in advance.
  2 个评论
Fadi Tantish
Fadi Tantish 2021-6-25
line 30 :
err(:,i) = abs(X-Xold);
warning :
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.
> In NR_Psat_fu (line 30)

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2021-6-25
The actual warning is about the \ on the line above. Your j matrix includes a NAN or else your j matrix includes infinities.
  5 个评论
显示 3更早的评论
Walter Roberson
Walter Roberson 2021-6-26
Ran in:
Let us check:
T = 400;
syms X [1 2]
[F,J] = Psat_fNR(X,T);
d11 = diff(F(1),X(1));
d12 = diff(F(1),X(2));
d21 = diff(F(2),X(1));
d22 = diff(F(2),X(2));
simplify(J(1,1)-d11)
ans = 
0
simplify(J(1,2)-d12)
ans = 
0
simplify(J(2,1)-d21)
ans = 
simplify(J(2,2)-d22)
ans = 
%the following code has been copied without change
function [f,j] = Psat_fNR(X,T)
R=8.3145;
a = 553661.24;
b = 30.49;
vl = X(1);
vv = X(2);
f(1,1) = (R*T/(vl-b) - a/(vl^2)) - (R*T/(vv-b) - a/(vv^2));
f(2,1) = ((1/(1-(b/vl))-a/(R*T*vl))-1 - log((1/(1-(b/vl))-a/(R*T*vl))-((b*(R*T/(vl-b) - a/(vl^2)))/(R*T))) - ((a*(R*T/(vl-b) - a/(vl^2)))/(R*T)^2)/(1/(1-(b/vl))-a/(R*T*vl))) - ((1/(1-(b/vv))-a/(R*T*vv))-1 - log((1/(1-(b/vv))-a/(R*T*vv))-((b*(R*T/(vv-b) - a/(vv^2)))/(R*T))) - ((a*(R*T/(vv-b) - a/(vv^2)))/(R*T)^2)/(1/(1-(b/vv))-a/(R*T*vv))) ;
% first derivatives of f(1,1) with respect to vl and vv.
df_vl= (2*a)/vl^3 - (R*T)/(b - vl)^2;
df_vv = (R*T)/(b - vv)^2 - (2*a)/vv^3;
% first derivatives of f(2,1) with respect to vl and vv.
dfug_vl = a/(R*T*vl^2) - b/(vl^2*(b/vl - 1)^2) - (b/(vl^2*(b/vl - 1)^2) + (b*((2*a)/vl^3 - (R*T)/(b - vl)^2))/(R*T) - a/(R*T*vl^2))/(1/(b/vl - 1) - (b*(a/vl^2 + (R*T)/(b - vl)))/(R*T) + a/(R*T*vl)) + (a*((2*a)/vl^3 - (R*T)/(b - vl)^2))/(R^2*T^2*(1/(b/vl - 1) + a/(R*T*vl))) + (a*(b/(vl^2*(b/vl - 1)^2) - a/(R*T*vl^2))*(a/vl^2 + (R*T)/(b - vl)))/(R^2*T^2*(1/(b/vl - 1) + a/(R*T*vl))^2);
dfug_vv = (b/(vv^2*(b/vv - 1)^2) + (b*((2*a)/vv^3 - (R*T)/(b - vv)^2))/(R*T) - a/(R*T*vv^2))/(1/(b/vv - 1) - (b*(a/vv^2 + (R*T)/(b - vv)))/(R*T) + a/(R*T*vv)) + b/(vv^2*(b/vv - 1)^2) - a/(R*T*vv^2) - (a*((2*a)/vv^3 - (R*T)/(b - vv)^2))/(R^2*T^2*(1/(b/vv - 1) + a/(R*T*vv))) - (a*(b/(vv^2*(b/vv - 1)^2) - a/(R*T*vv^2))*(a/vv^2 + (R*T)/(b - vv)))/(R^2*T^2*(1/(b/vv - 1) + a/(R*T*vv))^2);
j = [df_vl, df_vv;
dfug_vl, dfug_vv];
end
This tells us that you have correctly coded df_vl and df_vv in your Psat_fNR, but that your dfug_vl and dfug_vv are not correct partial derivatives of f(2,1)
Fadi Tantish
Fadi Tantish 2021-6-29
I appreciate your explanation on how to check the structure of the matrix. I fixed the derivatives of f(2,1) and now all answers equal to 0. However, I still face the same warning
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.
It seems that there is something going worng in replacing the variables and make the Jacobian matrix closes to infinite after 2 or 3 iterations instead of closing to zero.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Nonlinear Optimization 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by