sympy solve vs. solveset vs. nsolve

Question

I am trying to solve the following equation for r:

from sympy import pi, S, solve, solveset, nsolve, symbols

(n_go, P_l, T, gamma_w, P_g, r, R_mol) = symbols(
    'n_go, P_l, T, gamma_w, P_g, r, R_mol', real=True)

expr = -P_g + P_l - 3*R_mol*T*n_go/(4*r**3*pi) + 2*gamma_w/r
soln = solveset(expr, r, domain=S.Reals)
soln1 = solve(expr, r)

soln is of the form Complement(Intersection(FiniteSet(...))), which I really don't know what to do with. soln1 is a list of 3 expressions, two of which are complex. In fact, if I substitute values for the symbols and compute the solutions for soln1, all are complex:

vdict = {n_go: 1e-09, P_l: 101325, T: 300, gamma_w: 0.07168596252716256, P_g: 3534.48011713030, R_mol: 8.31451457896800}
for result in soln1:
    print(result.subs(vdict).n())

returns:

-9.17942953565355e-5 + 0.000158143657514283*I
-9.17942953565355e-5 - 0.000158143657514283*I
0.000182122477993494 + 1.23259516440783e-32*I

Interestingly, substituting values first and then using solveset() or solve() gives a real result:

solveset(expr.subs(vdict), r, domain=S.Reals).n()

{0.000182122477993494}

Conversely, nsolve fails with this equation, unless the starting point contains the first 7 significant digits of the solution(!):

nsolve(expr.subs(vdict), r,0.000182122 )

ValueError: Could not find root within given tolerance. (9562985778.9619347103 > 2.16840434497100886801e-19)

It should not be that hard, here is the plot:

My questions:

1. Why is nsolve so useless here?
2. How can I use the solution returned from solveset to compute any numerical solutions?
3. Why can I not obtain a real solution from solve if I solve first and then substitute values?

Answer 1

Your expr is essentially a cubic equation.

Applying the subs before or after solving should not substantially change anything.

soln

soln is of the form Complement(Intersection(FiniteSet(<3 cubic solutions>), Reals), FiniteSet(0)) i.e. a cubic solution on a real domain excluding 0.

The following should give you a simple FiniteSet back but evalf does not seem to be implemented well for sets.

print(soln.subs(vdict).evalf())

Hopefully something will be done about it soon.

1

The reason why nsolve is not useful is because the graph is almost asymptotically vertical. According to your graph, the gradient is roughly 1.0e8. I don't think nsolve is useful for such steep graphs.

Plotting your substituted expression we get: Zooming out we get:

This is a pretty wild function and I suspect nsolve uses an epsilon that is to large to be useful in this situation. To fix this, you could provide more reasonable numbers that are closer to 1 when substituting. (Consider providing different units of measurement. eg. instead of meters/year consider km/hour)

2

It is difficult to tell you how to deal with the output of solveset in general because every type of set needs to be dealt with in different ways. It's also not mathematically sensible since soln.args[0].args[0].args[0] should give the first cubic solution but it forgets that this must be real and nonzero.

You can use args or preorder_traversal or things to navigate the tree. Also reading documentation of various sets should help. solve and solveset need to be used "interactively" because there are lots of possible outputs with lots of ways to understand it.

3

I believe soln1 has 3 solutions instead of 4 as you state. Otherwise, your loop would print 4 lines instead of 3. All of them are technically of complex (as is the nature with floats in Python). However, the third solution you provide has a very small imaginary component. To remove these kinds of finicky things, there is an argument called chop which should help:

for result in soln1:
    print(result.subs(vdict).n(chop=True))

One of the results is 0.000182122477993494 which looks like your root.

Answer 2

The answer from Maelstrom is good but I just want to add a few points.

The values you substitute are all floats and with those values the polynomial is ill-conditioned. That means that the form of the expression that you substitute into can affect the accuracy of the returned results. That is one reason why substituting values into the solution from solve does not necessarily give exactly the same value that you get from substituting before calling solve.

Also before you substitute the symbols it isn't possible for solve to know which of the three roots is real. That's why you get three solutions from solve(expr, r) and only one solution from solve(expr.subs(vdict), r). The third solution which is real after the substitution is the same (ignoring the tiny imaginary part) as returned by solve after the substitution:

In [7]: soln1[2].subs(vdict).n()                                                                                                                              
Out[7]: 0.000182122477993494 + 1.23259516440783e-32⋅ⅈ

In [8]: solve(expr.subs(vdict), r)                                                                                                                            
Out[8]: [0.000182122477993494]

Because the polynomial is ill-conditioned and has a large gradient at the root nsolve has a hard time finding this root. However nsolve can find the root if given a narrow enough interval:

In [9]: nsolve(expr.subs(vdict), r, [0.0001821, 0.0001823])                                                                                                   
Out[9]: 0.000182122477993494

Since this is essentially a polynomial your best bet is actually to convert it to a polynomial and use nroots. The quickest way to do this is using as_numer_denom although in this case that introduces a spurious root at zero:

In [26]: Poly(expr.subs(vdict).as_numer_denom()[0], r).nroots()                                                                                               
Out[26]: [0, 0.000182122477993494, -9.17942953565356e-5 - 0.000158143657514284⋅ⅈ, -9.17942953565356e-5 + 0.000158143657514284⋅ⅈ]

Answer 3

Here is an unswer to the underlying question: How to compute the roots of the above equation efficiently?

Based on the suggestion by @OscarBenjamin, we can do even better and faster by using Poly and roots instead of nroots. Below, sympy computes in no time the roots of the equation for 100 different values of P_g, while keeping everything else constant:

from sympy import pi, Poly, roots, solve, solveset, nsolve, nroots, symbols
(n_go, P_l, T, gamma_w, P_g, r, R_mol) = symbols(
    'n_go, P_l, T, gamma_w, P_g, r, R_mol', real=True)
vdict = {pi:pi.n(), n_go:1e-09, P_l:101325, T:300, gamma_w:0.0717, R_mol: 8.31451457896800}
expr = -P_g + P_l - 3*R_mol*T*n_go/(4*r**3*pi) + 2*gamma_w/r
expr_poly = Poly(expr.as_numer_denom()[0], n_go, P_l, T, gamma_w, P_g, r, R_mol, domain='RR[pi]')
result = [roots(expr_poly.subs(vdict).subs(P_g, val)).keys() for val in range(4000,4100)]

All that remains is to check if the solutions are fulfilling our conditions (positive, real). Thank you, @OscarBenjamin!

PS: Should I expand the topic above to include nroots and roots?