Newton-Raphson Method Errors

I need assistance with this problem. I have two functions, where f1 corresponds to the dimensionless van der Waals equation of state, and function 2 corresponds to the dimensionless chemical potential function. X1 corresponds to the liquid mole fraction, and X2 corresponds to the vapor mole fraction. Y is equal to T/Tc, and Z is equal to P/Pv and is also the dimensionless van der Waals equation of state. I have this code, and I need the Newton-Raphson method to solve this equation for values of 0.50 < y < 1.0. As an initial step, I used x = {0.85, 1.3}. I need to obtain an isotherm at the end of the code, but my code is encountering the following issue.

(*Define the GaussEliminationwithPivoting function as you provided*)
GaussEliminationwithPivoting[AA_, bb_] :=
Module[{A = AA, b = bb, n, maxIndex, maxValue, factor, x},
n = Length[b];
For[i = 1, i <= n, i++, maxIndex = i;
maxValue = Abs[A[[i, i]]];
For[j = i + 1, j <= n, j++,
If[Abs[A[[j, i]]] > maxValue, maxIndex = j;
maxValue = Abs[A[[j, i]]];];];
If[maxIndex != i, A[[{i, maxIndex}]] = A[[{maxIndex, i}]];
b[[{i, maxIndex}]] = b[[{maxIndex, i}]];];
For[j = i + 1, j <= n, j++, factor = A[[j, i]]/A[[i, i]];
A[[j]] -= factor*A[[i]];
b[[j]] -= factor*b[[i]];];];
x = Table[0, {n}];
For[i = n, i >= 1, i--, x[[i]] = b[[i]];
For[j = i + 1, j <= n, j++, x[[i]] -= A[[i, j]]*x[[j]];];
x[[i]] /= A[[i, i]];];
{x}];
(*Define the nonlinear system of equations F[x1,x2,\[Lambda]]*)

F[x1_, x2_, \[Lambda]_] := {-(9/x1) + 9/x2 + (
12 *x1* \[Lambda])/(-1 + (3* x1)) + (12 *x2* \[Lambda])/(
1 - (3 *x2)) - (4* \[Lambda] *Log[-1 + (3 *x1)]) + (4 *\[Lambda]*
Log[-1 + 3 x2]), -(3/x1^2) + 3/
x2^2 + (8* (1/(-1 + (3 *x1)) + 1/(1 - (3* x2))) *\[Lambda])}

\[Lambda] = 0.6;
phaseEnvelopeData = {};
\[Lambda]Max = 0.99;

\[Lambda]Increment = 0.01;
maxIter = 1000;
tolerance = 1*10^-6;

While[\[Lambda] <= \[Lambda]Max, iteration = 0;
x = {0.830914, 1.24295};
While[iteration < maxIter,
J = D[F[x1, x2, \[Lambda]], {{x1, x2}}] /. {x1 -> x[[1]],
x2 -> x[[2]]};
If[Det[J] == 0, Print["Singular Jacobian encountered."];
Break[];];
negF = -F[x[[1]], x[[2]], \[Lambda]];
{dx} = GaussEliminationwithPivoting[J, negF];
x = x + dx;
If[Norm[F[x[[1]], x[[2]], \[Lambda]]] < tolerance, Break[]];
iteration++;];
If[iteration < maxIter, AppendTo[phaseEnvelopeData, {\[Lambda], x}];
Print["Convergence achieved for y =", \[Lambda]];
Print["Solution:", x];,
Print["Convergence could not be achieved for y =", \[Lambda]];];
\[Lambda] = \[Lambda] + \[Lambda]Increment;]


The solution that the code is presenting is:

Convergence achieved for y =0.6

Solution:{0.94746,0.94746}

Convergence achieved for y =0.61

Solution:{0.947293,0.947293}

Convergence achieved for y =0.62

Solution:{0.947115,0.947115}

Convergence achieved for y =0.63

Solution:{0.946926,0.946926}

Convergence achieved for y =0.64

Solution:{0.946726,0.946726}

Convergence achieved for y =0.65

Solution:{0.946512,0.946512}

Convergence achieved for y =0.66

Solution:{0.946284,0.946284}

Convergence achieved for y =0.67

Solution:{0.946041,0.946041}

Convergence achieved for y =0.68

Solution:{0.945779,0.945779}

Convergence achieved for y =0.69

Solution:{0.945499,0.945499}

Convergence achieved for y =0.7

Solution:{0.945196,0.945196}

Convergence achieved for y =0.71

Solution:{0.944869,0.944869}

Convergence achieved for y =0.72

Solution:{0.944515,0.944515}

Convergence achieved for y =0.73

Solution:{0.94413,0.94413}

Convergence achieved for y =0.74

Solution:{0.94371,0.94371}

Convergence achieved for y =0.75

Solution:{0.94325,0.94325}

Convergence achieved for y =0.76

Solution:{0.942744,0.942744}

Convergence achieved for y =0.77

Solution:{0.942185,0.942185}

Convergence achieved for y =0.78

Solution:{0.941564,0.941564}

Convergence achieved for y =0.79

Solution:{0.940871,0.940871}

Convergence achieved for y =0.8

Solution:{0.940093,0.940093}

Convergence achieved for y =0.81

Solution:{0.939212,0.939212}

Convergence achieved for y =0.82

Solution:{0.938209,0.938209}

Convergence achieved for y =0.83

Solution:{0.937057,0.937057}

Convergence achieved for y =0.84

Solution:{0.935722,0.935722}

Convergence achieved for y =0.85

Solution:{0.934158,0.934158}

Convergence achieved for y =0.86

Solution:{0.932305,0.932305}

Convergence achieved for y =0.87

Solution:{0.930086,0.930086}

Convergence achieved for y =0.88

Solution:{0.927396,0.927396}

Convergence achieved for y =0.89

Solution:{0.924108,0.924108}

Convergence achieved for y =0.9

Solution:{0.920115,0.920115}

Convergence achieved for y =0.91

Solution:{0.915581,0.915581}

Convergence achieved for y =0.92

Solution:{0.912689,0.912689}

Convergence achieved for y =0.93

Solution:{0.940159,0.940159}

Convergence achieved for y =0.94

Solution:{0.772379,0.768414}

Convergence achieved for y =0.95

Solution:{0.789519,0.784417}

General::munfl: 1/(2.33384*10^102+2.16017*10^100 I)^3 is too small to represent as a normalized machine number; precision may be lost.

General::munfl: 1/(-5.43908*10^104-4.91968*10^102 I)^3 is too small to represent as a normalized machine number; precision may be lost.

General::munfl: 1/(1.29739*10^107-5.62204*10^104 I)^3 is too small to represent as a normalized machine number; precision may be lost.

General::stop: Further output of General::munfl will be suppressed during this calculation.

Singular Jacobian encountered.

Convergence achieved for y =0.96

Solution:{-2.743745564859960*10^324+2.037581766232606*10^320 I,0.585937 +0.217491 I}

Convergence achieved for y =0.97

Solution:{4.32442,4.32442}

Convergence achieved for y =0.98

Solution:{0.775539,1.3761}

Convergence achieved for y =0.99

Solution:{0.830914,1.24295}


And I need to plot this graph

z[x_, y_] := ((8*y)/((3*x) - 1)) - (3/x^2)
resultsXl =
Table[{phaseEnvelopeData[[i, 2, 1]],
z[phaseEnvelopeData[[i, 2, 1]], phaseEnvelopeData[[i, 1]]]}, {i,
Length[phaseEnvelopeData]}];

resultsXv =
Table[{phaseEnvelopeData[[i, 2, 2]],
z[phaseEnvelopeData[[i, 2, 2]], phaseEnvelopeData[[i, 1]]]}, {i,
Length[phaseEnvelopeData]}];
LL = { "X1", "X2"};
LegLL = Table[Style[LL[[i]], FontSize -> 16], {i, 1, Length[LL]}];
ListPlot[{resultsXl, resultsXv}, Frame -> True, Joined -> True,
LabelStyle -> Directive[16], FrameStyle -> Directive[16],
FrameLabel -> {"V/Vc", "P/Pc"},
PlotStyle -> {Directive[16, Red], Directive[16, Black]},
FrameTicksStyle -> Directive[16], FrameTicks -> Automatic,
PlotLegends ->
Placed[LineLegend[LegLL, LegendMarkerSize -> 16,
LegendLayout -> {"Column", 1}], {Right, .30}], PlotRange -> All]


The graph need to seems like that

The left branch of the curve represents the molar volume of coexisting liquid and the right branch represents the molar volume of the vapor. The points defining the curve are generated by solving the equilibrium equation for various temperatures. Note that are two curves (left and right side)

But I'm obtaining that

EDIT: I rewrite this code with was confused and wrong but I'm still having problems as can't obtain a response

f1 = -(3/xl^2) + 3/xv^2 + (8 *(1/(-1 + (3 *xl)) + 1/(1 - (3* xv))) *y);
f2 = (-3 (3 - 9 xl + 4 xl^2 y) + 4 xl (-1 + 3 xl) y Log[-1 + 3 xl])/(
xl (-1 + 3 xl) (3 + Log[16])) - (-3 (3 - 9 xv + 4 xv^2 y) +
4 xv (-1 + 3 xv) y Log[-1 + 3 xv])/(xv (-1 + 3 xv) (3 + Log[16]));

f = {f1, f2};
J = D[f, {{xl, xv}}];
dfdy = D[f, y];
yIter = 0.999;
iterations = 100;
tolerance = 1*^-6;
step = 0.001;

size = Length[Range[0.5, 1, step]];
Z = ConstantArray[0, {2, size}];
X = ConstantArray[0, {2, size}];
xInit = {0.89, 1.20};

For[index = 1, index <= size, index++, yIter = 1 - index*step;
For[iter = 1, iter <= iterations, iter++,
JEval = N[J /. {xl -> xInit[[1]], xv -> xInit[[2]], y -> yIter}];
JInv = Inverse[JEval];
fEval = N[f /. {xl -> xInit[[1]], xv -> xInit[[2]], y -> yIter}];
deltaX = JInv.-fEval;
xInit = xInit + deltaX;
If[Norm[fEval] < tolerance, X[[All, index]] = xInit;
Z[[All, index]] = 8*yIter/(3*xInit[[1]] - 1) - (3/xInit[[1]]^2);
deltaX =
JInv.(step*
D[dfdy /. {xl -> xInit[[1]], xv -> xInit[[2]], y -> yIter}]);
xInit = xInit + deltaX;
Break[];];
If[iter == 99, Print["Did not converge"]];];];

Export["vle.svg",
ListLinePlot[Transpose[{X[[1]], Z[[1]]}], PlotStyle -> Black]];
Export["vle.svg",
ListLinePlot[Transpose[{X[[2]], Z[[2]]}], PlotStyle -> Black]]

• Overlaying two contour plots for fixed values of $\lambda$ is informative: Show[ContourPlot[F[x1, x2, 0.6][[1]], {x1, -10, 10}, {x2, -10, 10}, PlotPoints -> 100, Contours -> {0}, ContourShading -> None, ContourStyle -> {{Thickness[0.01], Red}}], ContourPlot[F[x1, x2, 0.6][[2]], {x1, -10, 10}, {x2, -10, 10}, PlotPoints -> 100, Contours -> {0}, ContourShading -> None]].
– JimB
Commented Nov 7, 2023 at 18:26

The built in function FindRoot can solve your equations. I believe it uses Newton's method by default, though you can choose other options. For your equations, I found singularities were occurring quite often, so I tried a different method, AffineCovariantNewton, which seemed to work better. In your results above, your values for x1 and x2 are equal for all lambda (or y) values. I get similar results, though with no error messages regarding a singular Jacobian (i.e. your failed convergence cases):

data = Table[{Values@
FindRoot[{-(9/x1) +
9/x2 + (12*x1*\[Lambda])/(-1 + (3*x1)) + (12*
x2*\[Lambda])/(1 - (3*x2)) - (4*\[Lambda]*
Log[-1 + (3*x1)]) + (4*\[Lambda]*
Log[-1 + 3 x2]), -(3/x1^2) +
3/x2^2 + (8*(1/(-1 + (3*x1)) +
1/(1 - (3*x2)))*\[Lambda])}, {x1, .83}, {x2, 1.3},
Method -> "AffineCovariantNewton"], \[Lambda]}, {\[Lambda], 0.5,
0.9, .01}]


{{{0.956226, 0.956226}, 0.5}, {{0.956058, 0.956058}, 0.51}, {{0.955882, 0.955882}, 0.52}, {{0.955697, 0.955697}, 0.53}, {{0.955502, 0.955502}, 0.54}, {{0.955296, 0.955296}, 0.55}, {{0.955079, 0.955079}, 0.56}, {{0.95485, 0.95485}, 0.57}, {{0.954606, 0.954606}, 0.58}, {{0.954348, 0.954348}, 0.59}, {{0.954074, 0.954074}, 0.6}, {{0.953782, 0.953782}, 0.61}, {{0.95347, 0.95347}, 0.62}, {{0.953137, 0.953137}, 0.63}, {{0.95278, 0.95278}, 0.64}, {{0.952396, 0.952396}, 0.65}, {{0.951983, 0.951983}, 0.66}, {{0.951537, 0.951537}, 0.67}, {{0.951054, 0.951054}, 0.68}, {{0.950528, 0.950528}, 0.69}, {{0.949955, 0.949955}, 0.7}, {{0.949328, 0.949328}, 0.71}, {{0.948639, 0.948639}, 0.72}, {{0.947877, 0.947877}, 0.73}, {{0.947032, 0.947032}, 0.74}, {{0.94609, 0.94609}, 0.75}, {{0.945032, 0.945032}, 0.76}, {{0.943836, 0.943836}, 0.77}, {{0.942475, 0.942475}, 0.78}, {{0.935061, 0.935061}, 0.79}, {{0.940222, 0.940222}, 0.8}, {{0.945379, 0.945379}, 0.81}, {{0.950496, 0.950496}, 0.82}, {{0.955527, 0.955527}, 0.83}, {{0.960413, 0.960413}, 0.84}, {{0.965074, 0.965074}, 0.85}, {{0.969406, 0.969406}, 0.86}, {{0.973271, 0.973271}, 0.87}, {{0.97648, 0.97648}, 0.88}, {{0.978759, 0.978759}, 0.89}, {{0.97968, 0.97968}, 0.9}}

These do in fact seem to be solutions. Plugging in lambda,x1,x2 for a random case gives {0,0} to within numerical precision:

{-(9/x1) +
9/x2 + (12*x1*\[Lambda])/(-1 + (3*x1)) + (12*
x2*\[Lambda])/(1 - (3*x2)) - (4*\[Lambda]*
Log[-1 + (3*x1)]) + (4*\[Lambda]*Log[-1 + 3 x2]), -(3/x1^2) +
3/x2^2 + (8*(1/(-1 + (3*x1)) +
1/(1 - (3*x2)))*\[Lambda])} /. {\[Lambda] -> .71,
x1 -> .949328, x2 -> .949328}

(* {-8.88178*10^-16, 0.} *)


So it does seem to be the case the x1=x2 for all lambdas in this range, at least for the initial guess you've chosen. Since numerical root finders tend to find local minima, this may be a consequence of the initial guesses you've chosen.

Here's the plot for reference:

  processEnvelopeData[{x1_, x2_},
lambda_] := {{x1, z[x1, lambda]}, {x2, z[x2, lambda]}};
z[x_, y_] := ((8*y)/((3*x) - 1)) - (3/x^2);
ListLinePlot[Transpose[processEnvelopeData @@@ data]]


This is actually two plots. Since x1=x2, the two curves sit right on top of one another.

EDIT: Your notation confused me at first, but after rereading your question more carefully, something occurred to me. Your quantities are dimensionless, but leaving the reference points out, your x's are volume (V), your y (lambda) is temperature (T), and your z[x,y] is pressure (P = P(V,T)). You're trying to plot an isotherm, which implies temperature is held fixed, but you are varying the temperature. What I believe you need to do is choose a fixed temperature (i.e. lambda or y), and vary x over the range of interest to reproduce the plot you showed. Something like this for example, for a generic lambda of 0.8. Other values give similar results. Of course physically, you want to restrict the range to get positive pressure.