3
$\begingroup$

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:

  1. $\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=\frac{3}{v_{r,3}} - \frac{3}{v_{r,1}} + \frac{8*0.85}{3}\, \, 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}]

enter image description here

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?

$\endgroup$
  • 3
    $\begingroup$ The difference is that Plot holds it's first argument, while FindRoot does not. The simplest solution is to add an argument restriction, i.e., f[pr_?NumericQ, tr_] := Module[...]. $\endgroup$ – Carl Woll Sep 19 '19 at 16:14
  • 1
    $\begingroup$ Guess what? It works! Thank you. Why did that happen though? $\endgroup$ – Assaad Mrad Sep 19 '19 at 16:16
  • 4
    $\begingroup$ support.wolfram.com/kb/12502 $\endgroup$ – user6014 Sep 19 '19 at 16:20
12
$\begingroup$

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.

$\endgroup$

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.