Mfcttrf

Slight background: I am trying to find the coexistence reduced pressure ($p_r$) of a Van der Waals fluid at a reduced temperature of 0.85, this requires finding the solution to the following system of equations where $p_r$ is the only control variable:

1. $p_r = -(3/v_r^2) + (8 * 0.85)/(-1 + 3 v_r)$

Of course, the above equation is a cubic in vr, and you get three positive roots for certain values of $0<pr<1$. I label the smallest of these roots $v_r,1$ up
to $v_r,3$. Then, the second equation to be satisfied is the following integral equation:

2. $int_v_r,1^v_r,3 left[-(3/v_r^2) + frac(8 * 0.85)(-1 + 3 v_r)right] dv_r = p_r(v_r,3-v_r,1)$

The integral turns out to have an analytical solution, so to make things easier, I write:

$int_v_r,1^v_r,3 left[-(3/v_r^2) + frac(8 * 0.85)(-1 + 3 v_r)right] dv_r=frac3v_r,3 - frac3v_r,1 + frac8*0.853, , Log[frac-1 + 3 v_r,3-1 + 3 v_r,1]$

To do that on Mathematica, I introduce the following module:

f[pr_, tr_] := Module[vr1, vr2, vr3,
 sol = Sort[vr /. NSolve[pr == -(3/vr^2) + (8*tr)/(-1 + 3 vr), vr]];
 vr1 = sol[[1]]; vr2 = sol[[2]]; vr3 = sol[[3]];
 Re[(3/vr + 8/3 tr Log[-1 + 3 vr] /. 
 vr -> vr3) - (3/vr + 8/3 tr Log[-1 + 3 vr] /. vr -> vr1) - 
 pr*(vr3 - vr1) ]
 ]

Then I plot f[pr,0.85] to find a graphical solution:

Plot[f[pr, 0.85], pr, 0.4, 0.55]

So it is clear that the solution is pr = 0.505 (also independently checked).

However, watch what happens when I use FindRoot:

FindRoot[f[pr, 0.85], pr, 0.4, 0.1, 0.6, Method -> "Newton"]

I get:

pr -> 0.530736

I've tried changing the Method and the search intervals to no avail.

Can someone please help with this?

The problem is that FindRoot evaluates f[pr, .8] where pr is just a symbol, it is not a number. Then, it replaces pr with a number. However, the order of the roots when pr is a symbol differs from the order of the roots when pr is a number. Here is a function that just solves your first equation:

g[pr_, tr_] := Sort[vr/.NSolve[pr==-(3/vr^2)+(8*tr)/(-1+3 vr),vr]];

Now, let's evaluate g with a symbolic pr:

g[pr, .85] //Short

-((0.0222222 (-34.-<<1>>))/pr)+((<<21>> -<<1>>) <<1>>)/(pr <<1>>^<<1>>)-((<<1>>) <<1>>)/pr,<<1>>,<<1>>

It's rather complicated, so I truncated the output. Now compare substituting a value for pr into the above expression versus just evaluating g with a numeric pr:

g[pr, .85] /. pr -> .5 //Chop
g[.5, .85]

1.13734, 0.553756, 3.17557

0.553756, 1.13734, 3.17557

The order is different. This is why your FindRoot code doesn't work. Defining f so that it doesn't evaluate for symbolic pr fixes your problem:

Clear[f]
f[pr_?NumericQ, tr_] := Module[vr1, vr2, vr3,
 sol = Sort[vr /. NSolve[pr == -(3/vr^2) + (8*tr)/(-1 + 3 vr), vr]];
 vr1 = sol[[1]]; vr2 = sol[[2]]; vr3 = sol[[3]];
 Re[
 (3/vr + 8/3 tr Log[-1 + 3 vr] /. vr -> vr3) - 
 (3/vr + 8/3 tr Log[-1 + 3 vr] /. vr -> vr1) - 
 pr*(vr3 - vr1)
 ]
]
FindRoot[f[pr, .85], pr, .5]

pr -> 0.504492

Another possibility is to use Solve with the Cubics->False option instead of NSolve. This is because Solve will return Root objects, and Mathematica knows how to sort Root objects even for symbolic arguments:

Clear[f]
f[pr_, tr_] := Module[vr1, vr2, vr3,
 sol = vr /. Quiet @ Solve[pr == -(3/vr^2) + (8*tr)/(-1 + 3 vr), vr, Cubics->False];
 vr1 = sol[[1]]; vr2 = sol[[2]]; vr3 = sol[[3]];
 Re[
 (3/vr + 8/3 tr Log[-1 + 3 vr] /. vr -> vr3) - 
 (3/vr + 8/3 tr Log[-1 + 3 vr] /. vr -> vr1) - 
 pr*(vr3 - vr1)
 ]
]
FindRoot[f[pr, .85], pr, .5]

pr -> 0.504492

(the Quiet is to avoid messages related to using exact methods with inexact coefficients)

Let's examine the output of Solve:

Solve[pr == -(3/vr^2) + (8*tr)/(-1 + 3 vr), vr, Cubics->False]

vr -> Root[-3 + 9 #1 + (-pr - 8 tr) #1^2 + 3 pr #1^3 &, 1], vr ->
Root[-3 + 9 #1 + (-pr - 8 tr) #1^2 + 3 pr #1^3 &, 2], vr ->
Root[-3 + 9 #1 + (-pr - 8 tr) #1^2 + 3 pr #1^3 &, 3]

Notice how each solution is just a Root object with an index. By default, the indices correspond to roots with increasing real part. This is what you were trying to achieve with your Sort, so it isn't necessary in the above code.