Sympy nsolve for simultaneous equation chemical equilibria

Question

I am trying to use Sympy's nsolve as a numerical method to solve reaction progress variables in simultaneous chemical equilibria. Unfortunately, for the present case no amount of initial guessing appears to arrive at a satisfactory solution. I believe this system does have a suitable solution where all four progress variables (the unknowns) have real solutions between 0 and 1.

Here is the simplified code:

from sympy import *
var('x1,x2,x3,x4')
P = 1.          # Pressure [bar]
n1 = 2.*x1      # x1 is progress variable for reaction 1
n2 = x2         # x2 is progress variable for reaction 2
n3 = 2.*x3      # x3 is progress variable for reaction 3
n4 = x4         # x4 is progress variable for reaction 4
n5 = 1.-x1-x2-3.*x3-2.*x4
nTot = n1+n2+n3+n4+n5

rxn1 = Eq(P*n1**2/(nTot*n5),3.06e-3)
rxn2 = Eq(n2/n5,8.02e5)
rxn3 = Eq(nTot*n3**2/(P*n5**3),1.07e21)
rxn4 = Eq(nTot*n4/(P*n5**2),8.25e14)

reactions = (rxn1,rxn2,rxn3,rxn4)
unknowns = (x1,x2,x3,x4)
answer = nsolve(reactions,unknowns,(1e-4,.1,.3,.2))
print (answer)

Here is the output from Sympy (I am using Sympy v1.1.1 in a Spyder3 environment with 64 bit Python 2.7):

Traceback (most recent call last):

  File "<ipython-input-9-43632f4fa60f>", line 18, in <module>
    answer = nsolve(reactions,unknowns,(1e-4,.1,.3,.2))

  File "C:\Users\<user>\AppData\Local\Continuum\anaconda2\lib\site-packages\sympy\utilities\decorator.py", line 91, in func_wrapper
    return func(*args, **kwargs)

  File "C:\Users\<user>\AppData\Local\Continuum\anaconda2\lib\site-packages\sympy\solvers\solvers.py", line 2847, in nsolve
    x = findroot(f, x0, J=J, **kwargs)

  File "C:\Users\<user>\AppData\Local\Continuum\anaconda2\lib\site-packages\mpmath\calculus\optimization.py", line 975, in findroot
    % (norm(f(*xl))**2, tol))

ValueError: Could not find root within given tolerance. (1.1449e+42 > 2.1684e-19)
Try another starting point or tweak arguments.

I have attempted to use various solve (unsurprisingly, no analytic solution is available) and nsolve methods, none of which have thusfar been successful for this system of non-linear equations. I have also tried recasting the equations in logarithmic form with similar results.

Are there other solver settings that you recommend I attempt to employ?

Thank you for any assistance! Dave

Answer 1

Part of the difficulty in finding a solution is due to the wide variation in magnitudes of the constants. For example, rxn3 has a value on the right-hand side of 1070000000000000000000, while rxn1 has a right-hand side of 0.00306.
Numeric algorithms often have trouble with equations with extreme floating-point values.

Nevertheless, we can get somewhat close by symbolically solving the 4 equations for n1, n2, n3, n4. (Finding a solution in terms of the n-variables is essentially the same as solving for the x-variables, but allows us to address the main problem with fewer variables/equations).

soln = sym.solve([rxn1, rxn2, rxn3, rxn4], [n1, n2, n3, n4], dict=True)[0]
# {n1: -0.0553172667437573*sqrt(n5*nTot),
#  n2: 802000.0*n5,
#  n3: -32710854467.5922*sqrt(n5**3/nTot),
#  n4: 825000000000000.0*n5**2/nTot}

This gives us n1, n2, n3 and n4 in terms of n5 and nTot. We can now attempt to solve the two remaining equations which define n5 and nTot -- rxn5 and rxn6, (see below) -- for n5 and nTot.

Unfortunately, sym.solve fails to find a symbolic solution:

In [114]: sym.solve([rxn5, rxn6], [n5, nTot])
Out[114]: []

Note that when sym.solve returns an empty list it does not necessarily imply there are no solutions, only that sym.solve failed to find one.

So let's try to find a numeric solution. Per the docs for sym.nsolve,

mpmath.findroot is used and you can find there more extensive documentation, especially concerning keyword parameters and available solvers.

and the docs for mpmath.findroot describe a maxstep parameter. By increasing the maxstep parameter, it turns out sym.nsolve can find a solution for n5 and nTot:

In [148]: sym.nsolve([rxn5, rxn6], [n5, nTot], [0.5, 0.5])
ValueError: Could not find root within given tolerance. (19920803.9870484600143 > 2.16840434497100886801e-19)
Try another starting point or tweak arguments.

In [149]: sym.nsolve([rxn5, rxn6], [n5, nTot], [0.5, 0.5], maxsteps=50)
Out[149]: 
Matrix([
[1.83437114233209e-8],
[  0.477964294180007]])

Now we can check if the solution makes sense:

import sympy as sym

P = 1.0          # Pressure [bar]
n1, n2, n3, n4, n5, nTot = sym.symbols('n1,n2,n3,n4,n5,nTot', positive=True, real=True)
rxn1 = sym.Eq(P*n1**2/(nTot*n5),3.06e-3)
rxn2 = sym.Eq(n2/n5,8.02e5)
rxn3 = sym.Eq(nTot*n3**2/(P*n5**3),1.07e21)
rxn4 = sym.Eq(nTot*n4/(P*n5**2),8.25e14)
rxn5 = sym.Eq(n5, 1.0-(n1/2.0)-n2-3.0*(n3/2.0)-2.0*n4)
rxn6 = sym.Eq(nTot, n1+n2+n3+n4+n5)

soln = sym.solve([rxn1, rxn2, rxn3, rxn4], [n1, n2, n3, n4], dict=True)[0]
rxn5 = rxn5.subs(soln)
rxn6 = rxn6.subs(soln)

expr5 = rxn5.rhs - rxn5.lhs
expr6 = rxn6.rhs - rxn6.lhs
print(expr5)
print(expr6)
# -1.65e+15*n5**2/nTot - 802001.0*n5 + 0.0276586333718787*sqrt(n5*nTot) + 49066281701.3884*sqrt(n5**3/nTot) + 1.0
# 825000000000000.0*n5**2/nTot + 802001.0*n5 - nTot - 0.0553172667437573*sqrt(n5*nTot) - 32710854467.5922*sqrt(n5**3/nTot)

root = sym.nsolve([rxn5, rxn6], [n5, nTot], [0.5, 0.5], maxsteps=100)

# or using scipy.optimize.fsolve
# import scipy.optimize as optimize
# f = sym.lambdify([n5, nTot], (expr5, expr6))
# root = optimize.fsolve(lambda x: f(*x), [[0.5, 1e7]], xtol=1e-16)
# print('f({}) = {}'.format(root, f(*root)))

# or using bisect with verify=False -- mentioned in the `sym.nsolve` docs
# root = sym.nsolve([rxn5, rxn6], [n5, nTot], [0.5, 0.5], solver='bisect', verify=False)

print(root)
# Matrix([[1.64738651107199e-8], [0.530353285064842]])

root = {n5:root[0], nTot:root[1]}
nsoln = {var:val.subs(root) for var,val in soln.items()}
nsoln.update(root)
for var in [n1, n2, n3, n4, n5, nTot]:
    print('{} = {}'.format(var, nsoln[var]))
for i, rxn in enumerate((rxn1, rxn2, rxn3, rxn4, rxn5, rxn6), 1):
    print('rxn{}: {}'.format(i, (rxn.lhs-rxn.rhs).subs(nsoln)))

which yields this solution:

n1 = 0.00000517060185532663
n2 = 0.0132120398187973
n3 = 0.0949735158684791
n4 = 0.422162542301846
n5 = 1.64738651107199E-8
nTot = 0.530353285064842

and these errors (differences between the left- and right-hand sides of the respective equations):

rxn1: -8.67361737988404E-19
rxn2: 0
rxn3: 131072.000000000
rxn4: -0.125000000000000
rxn5: -2.22044604925031E-16
rxn6: 5.55111512312578E-17

Note that while the absolute error for rxn3 is rather large, the relative error is < 2e-16:

In [13]: 131072/1.07e21
Out[13]: 1.2249719626168225e-16

---

By the way, the sym.nsolve docs also say:

Note, however, that functions which are very steep near the root the verification of the solution may fail. In this case you should use the flag verify=False and independently verify the solution.

and

One might safely skip the verification if bounds of the root are known and a bisection method is used:

In [155]: sym.nsolve([rxn5, rxn6], [n5, nTot], [0.5, 0.5], verify=False, solver='bisect')
Out[155]: 
Matrix([
[2.77113579502186e-9],
[2.83553974806881e-6]])

but this solution does not turn out to be quite as good as the solution found with maxsteps=50 (or greater).