An issue while applying Newton method for solving nonlinear systems.

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.

[Walter Roberson]

What is the warning message and which is line 30?

[Fadi Tantish]

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)

Answer 1

[Walter Roberson]

The actual warning is about the \ on the line above. Your j matrix includes a NAN or else your j matrix includes infinities.

[Fadi Tantish]

the j matrix (Jacobian matrix) was obtained using diff function on f(1,1) and f(2,1) with respect to the variables. Do you have any tip on how to check if it has a value or not?

[Walter Roberson]

dbstop if naninf

[Fadi Tantish]

I am sorry. I am new to matlab I didn't understand what you say? is there a function to check for this nan or inf?

[Walter Roberson]

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]

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.