1
$\begingroup$

I'm trying to solve the following system of equations and I can't seem to manipulate Mathematica to be get an answer. Any and all help would be appreciated.

Some notes: XNil and XNis should both be between 0 and 1 when solved. The value of 1700 in the equation should be able to range from 1356 to 1726.

HmCu = 13055;
TmCu = 1356;
HmNi = 17472;
TmNi = 1726;

GmCu[T_] = HmCu - T*(HmCu/TmCu)
GmNi[T_] = HmNi - T*(HmNi/TmNi)

R = 8.3144;
Lol = 11100;
Los = 10435;

In[30]= Rationalize[{GmCu[1700] + R*1700*Log[(1 - XNil)/(1 - XNis)] + 
Lol*XNil^2 - Los*XNis^2 == 0 && 
GmNi[1700] + R*1700*Log[(XNil)/(XNis)] + Lol*(1 - XNil)^2 - 
Los*(1 - XNis)^2 == 0}]

Out[30]= {-(1122730/339) + 11100 XNil^2 - 10435 XNis^2 + 
353362/25 Log[(1 - XNil)/(1 - XNis)] == 0 && 
227136/863 + 11100 (1 - XNil)^2 - 10435 (1 - XNis)^2 + 
20786/25 Log[XNil/XNis] == 0}

Solve[Out[30], {XNil, XNis}, Reals]

Solve::nsmet: This system cannot be solved with the methods available to Solve.

EDIT: I was able to find a solution using FindRoot instead of Solve. However, both solutions are above 1 which does not make sense within the context of the problem.

$\endgroup$
  • $\begingroup$ Is In[30] supposed to be the same as Out[30], because they are definitely not. If not, what is the their significance to each other? $\endgroup$ – Bill Watts Dec 6 '17 at 7:58
1
$\begingroup$

Your problem arrises from the numerical difficulty of the solving your equations. To see the problem, let's plot the zero contours of the left hand side of your equations.

HmCu = 13055;
TmCu = 1356;
HmNi = 17472;
TmNi = 1726;
GmCu[T_] = HmCu - T*(HmCu/TmCu);
GmNi[T_] = HmNi - T*(HmNi/TmNi);
R = 8.3144;
Lol = 11100;
Los = 10435;

eqs = 
  Rationalize @
    {GmCu[1700] + R*1700*Log[(1 - XNil)/(1 - XNis)] + Lol*XNil^2 - Los*XNis^2 == 0, 
     GmNi[1700] + R*1700*Log[(XNil)/(XNis)] + Lol*(1 - XNil)^2 - Los*(1 - XNis)^2 == 0}

Quiet @ 
  ContourPlot[Evaluate @ eqs, {XNil, 0, 4}, {XNis, 0, 4}, 
    ContourStyle -> {Red, Blue}]

plot

There appears to be a zero near the upper right corner of the unit square, but the two traces are so nearly parallel that FindRoot is not going to do well when it uses Newton's method (the default method) for this problem. We will tell it to use the secant method by giving it a bracket of starting values rather than a single starting point. Further, we want the bracketing to fit rather closely around the root, and it will be necessary to work with arbitrary precision arithmetic.

To get good starting values, we zoom the plot in on the desired zero.

Quiet @ 
  ContourPlot[Evaluate @ eqs, {XNil, 0, 4}, {XNis, 0, 4},
    PlotRange -> {{.89, .91}, {.91, .93}},
    ContourStyle -> {Red, Blue}]

zoom_in

Now we can use FindRoot and get a good result.

FindRoot[eqs, {{XNil, .89`20, .91`20}, {XNis, .91`20, .93`20}}, 
  WorkingPrecision -> 40]
{XNil -> 0.9010823440563772363797513599485915326375, 
 XNis -> 0.9208285280785948177856956468479419307813}
$\endgroup$
  • $\begingroup$ +1, but this works for me without any fuss: FindRoot[eqs, {{XNil, 0.5}, {XNis, 0.51}}] $\endgroup$ – Michael E2 Dec 6 '17 at 12:41
  • $\begingroup$ @MichaelE2. Interesting finding. It did not occur to me that backing off the starting point toward the lower left corner of unit square, where the two traces are well separated, would cure the problem of using Newton's method with machine arithmetic. Live and learn. $\endgroup$ – m_goldberg Dec 6 '17 at 19:21
1
$\begingroup$

@m_goldberg's comment got me to wondering if I accidentally hit upon a good starting point with the following:

FindRoot[eqs, {{XNil, 0.5}, {XNis, 0.51}}]
(*  {XNil -> 0.901082, XNis -> 0.920829}  *)

It turns out the the OP's problem has an pretty infeasible region (red means the solution did not converge to the solution above):

With[{n = 1024, (* n = 256 is faster and looks ok; 128 takes 4+ sec. for me *)
  s0 = {XNil, XNis} /. FindRoot[eqs, {{XNil, 0.5}, {XNis, 0.51}}]},
 Graphics[{PointSize[0.0001], Red,
   Point[
     Flatten[
      ParallelTable[
       Block[{s = {XNil, XNis} /. Quiet@FindRoot[eqs, {{XNil, x}, {XNis, y}}]},
        If[EuclideanDistance[s, s0] > 0.001 // TrueQ, {x, y}, Nothing]
        ],
       {x, 1./n, 1. - 1./n, 1./n}, {y, 1./n, 1. - 1./n, 1./n}],
      1]
    ]
   },
  Frame -> True]
 ]

Mathematica graphics

$\endgroup$
  • $\begingroup$ +1 for the investigation. It was this kind of behavior that I feared from Newton's method and machine arithmetic. This result makes me feel better about the choices I made in my answer. $\endgroup$ – m_goldberg Dec 7 '17 at 1:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.