Maxwell construction problem with FindRoot

Question

A Maxwell construction is just an integral of area and then solving for a root to find vapor pressure.

This must be exactly the same value as where the Gibbs energy plot crosses at a point. I found that to be p = 0.81.

The problem is that my areadifferential function doesn't work properly; it doesn't return p = 0.81;

What mistake I have made in my coding?

t = 0.95; 

p[v_] := (8*t)/(3*v - 1) - 3/v^2;

Plot[p[v], {v, 0.5, 3}, 
  PlotRange -> {{0, 3}, {0.6, 1}}, 
  AxesLabel -> {V/Vc, P/Pc}]

g[v_] := (-t)*Log[3*v - 1] + 0.95/(3*v - 1) - 9/(4*v); 

ParametricPlot[{p[v], g[v]}, {v, .65, 2.25},
  AxesLabel -> {P/Pc, G/NKT},
  PlotLabel -> "Gibbs Free Energy Vs. P/Pc"]

That returns pressure p = P/Pc = 0.81 as a plot which is correct

This is the error Maxwell construction "area differential" part

pint[v_] := (8/3)*t*Log[3*v - 1] + 3/v;

areadifferential[p0_, v1guess_, v2guess_] := 
  (v1 = FindRoot[p[v] == p0, {v, v1guess}][[1,2]]; 
   v2 = FindRoot[p[v] == p0, {v, v2guess}][[1,2]];
pint[v2] - pint[v1] - p0*(v2 - v1))

FindRoot[areadifferential[p0, 0.7, 2] == 0, {p0, 0.8, 0.82}]

Answer 1

Good job on describing your problem.

Copying your definitions from your question

p[v_] := (8 t)/(3 v - 1) - 3/v^2
g[v_] := t/(3 v - 1) - t Log[-1 + 3 v] - 9/(4 v)
pint[v_] := 8/3 t Log[-1 + 3 v] + 3/v
t = 0.95

I found the value 0.81188 to be the value that minimized your areadifferential function (shown at the bottom).

Look at the plot of P/Pc vs V/Vc and note the three points corresponding to three different V/Vc values resulting in p[v] = 0.8118.

Module[
 {
  v1 = v /. FindRoot[p[v] == 0.8118, {v, 0.5}],
  p1,
  v2 = v /. FindRoot[p[v] == 0.8118, {v, 1}],
  p2,
  v3 = v /. FindRoot[p[v] == 0.8118, {v, 2}],
  p3
  },
 p1 = p[v1];
 p2 = p[v2];
 p3 = p[v3];
 Show[
  Plot[p[v], {v, 0.65, 2.25}, PlotRange -> {{0, 3}, {0.6, 1}}, 
   AxesLabel -> {"V/Vc", "P/Pc"}],
  Graphics[
   {
    PointSize[0.03],
    Red,
    Point[{v1, p1}],
    Green,
    Point[{v2, p2}],
    Black,
    Point[{v3, p3}]
    }
   ]
  ]
 ]

The problem that you are experiencing is because the value of 0.7 for v1guess is just to the right of the valley so it converges on the wrong point (i.e., the green point).

Change v1guess to 0.5 and you will be fine, it will converge on the red point. Using v2guess of 2.0 works fine causing convergence on the black point.

I slightly modified the areadifferential function to use Module rather than parenthesis. One needs to constrain the input arguments to be numeric so FindRoot doesn't complain (FindFoot first tries to work with symbolic arguments which doesn't work with your problem).

areadifferential[
p0_?NumericQ,
v1guess_?NumericQ,
v2guess_?NumericQ] := Module[
  {
   v1 = FindRoot[p[v] == p0, {v, v1guess}][[1, 2]],
   v2 = FindRoot[p[v] == p0, {v, v2guess}][[1, 2]]
   },
  pint[v2] - pint[v1] - p0*(v2 - v1)
  ]

and then

FindRoot[areadifferential[p0, 0.5, 2] == 0, {p0, 0.8, 0.82}]
(* {p0 -> 0.811879} *)