1

I am trying to solve a system of non-linear equations. The issue is that my input values are being returned to me as a valid solution. In order for this to happen, some equations have to be ignored. I have tried both sympy.solvers.nsolve() and scipy.optimize.fsolve(). Both give the same answers. I have posted both codes starting with scipy and then sympy. Lastly, I have posted example results.

from scipy.optimize import fsolve  # required library
import sympy as sp
import scipy
# known values hard coded for testing
a = 1
b = 3
c = a + b
d = 0
kp_NO = 0.00051621
kp_H2O = 0.0000000127
kp_CO2 = 0.00000001733
p = 5 # pressure
dec = 3 # decimal point precision

# Solving the system of equations
def equations(vars):
    (e, f, g, h, j, k, l) = vars
    f1 = e + j - a
    f2 = e + 2*g + 2*j + k + l - a - 2*c
    f3 = f + k - b
    f4 = 2*h + l - 2*d - (2*79/21)*c
    f5 = kp_NO - (l/(c*(79/21)*c))
    f6 = kp_H2O - (k/(b*sp.sqrt((c*p)/(e + f+ g + h + j + k + l))))
    f7 = kp_CO2 - (j/(a*sp.sqrt((c*p)/(e + f + g + h + j + k + l))))
    return[f1, f2, f3, f4, f5, f6, f7]
e, f, g, h, j, k, l = scipy.optimize.root_scalar(equations, (0, 0, 0, 0, 0, 0, 0))
# CO, H2, O2, N2, CO2, H2O, NO
print(e, f, g, h, j, k, l)

import sympy as sp # required library
import math as m
# known values hard coded for testing
a = 1
b = 3
c = a + b
d = 0
kp_NO = 0.0000000015346
kp_H2O = 1.308*10**-23
kp_CO2 = 9.499*10**-46
p = 5 # pressure
dec = 10 # decimal point precision

# Solving the system of equations
e, f, g, h, j, k, l = sp.symbols('e, f, g, h, j, k, l')
f1 = e + j - a
f2 = e + 2*g + 2*j + k + l - a - 2*c
f3 = f + k - b
f4 = 2*h + l - 2*d - (2*79/21)*c
f5 = kp_NO - (l / (c * (79 / 21) * c))
f6 = kp_H2O - (k / (b * sp.sqrt((c * p) / (e + f + g + h + j + k + l))))
f7 = kp_CO2 - (j / (a * sp.sqrt((c * p) / (e + f + g + h + j + k + l))))
# f5 = m.log(kp_NO, p) + h/2 + g/2 - l
# f6 = m.log(kp_H2O, p) + g/2 + f - k
# f7 = m.log(kp_CO2, p) + e + g/2 - j
results = sp.solvers.nsolve([f1, f2, f3, f4, f5, f6, f7], [e, f, g, h, j, k, l], [0.00004, 0.00004, 0.49, 3.76, 0.25, 0.75, 0.01], manual=True)

e = round(results[0], dec)
f = round(results[1], dec)
g = round(results[2], dec)
h = round(results[3], dec)
j = round(results[4], dec)
k = round(results[5], dec)
l = round(results[6], dec)
# CO, H2, O2, N2, CO2, H2O, NO
print(e, f, g, h, j, k, l)

1.00000000000000 3.00000000000000 3.9999999538 15.0476190014 0.0 0.0 9.24e-8
Thomas Murdoch
  • 13
  • 2
  • Since both SymPy and SciPy apparently give the same result I'm inclined to think that both are correct for the given input. Maybe the input isn't what you intended it to be? What makes you think that the output is incorrect? – Oscar Benjamin Dec 03 '21 at 22:10
  • 1
    The code that you show that uses scipy has errors, so it can not be run as it is to reproduce the behavior that you report. It will be easier for someone to help you if you provide a [minimal *reproducible* example](https://stackoverflow.com/help/minimal-reproducible-example). – Warren Weckesser Dec 03 '21 at 22:29
  • 1
    The solution the SymPy code produces indeed is a zero for the 7 equations you have implemented there. Why do you believe some equations are being ignored? – nonDucor Dec 03 '21 at 22:37
  • Hey all, thanks for the responses. I think I wasn't clear in stating that there should be multiple answers to these equations. Some of the Kp values are very small and as such may be rounded off? The smallest being on the order of 10e-46. As for a reproducible error, I am not sure why the scipy isn't working. sympy gives 0.2500000000 0.7500000000 0.9999999971 3.7619047590 0.0 0.0 5.8e-9 – Thomas Murdoch Dec 04 '21 at 02:17
  • I believe that some may be ignored because it is giving the inputs back to me as solutions. This may be because the kp data is rounded off? Also, this is supposed to describe the chemical reaction of syngas at varying compositions. The number of moles is supposed to change for things like CO and H2. However, my inputs of 0.75 and 0.25 (respectively) are being given back as answers – Thomas Murdoch Dec 04 '21 at 02:27
  • The scipy errors were solved, I was guessing using zeros but changing this to 0.1 removed the errors and gave the result: 0.2499999959636201 0.7499999911281063 0.9990290271172916 3.7609337954761903 4.036379920862849e-09 8.871893602197363e-09 0.0019419328571428573 – Thomas Murdoch Dec 04 '21 at 02:30

1 Answers1

1

You will also get no change from the initial guess if each of the equations is within the tolerance of the solver. And with the scaling of your equations, I suspect this is the case. Here is another approach:

Let's work with rationals instead of decimals:

>>> from sympy import nsimplify
>>> eqs = [nsimplify(i, rational=True) for i in [f1, f2, f3, f4, f5, f6, f7]]
>>> v = [e, f, g, h, j, k, l]

All but the last equation is easy to solve for all but g; there is 1 solution:

 >>> lpart = solve(eqs[:-1], set(v) - {g}, dict=True)[0]

When radicals are removed from the numerator of the last equation, it yields a cubic for which roots can be found

>>> from sympy import real_roots
>>> from sympy.solvers.solvers import unrad
>>> u, cov = unrad(eqs[-1].subs(lpart).as_numer_denom()[0]); assert cov == []
>>> gsol = list(real_roots(u))

It looks like the first two real roots are solutions of the last equation:

>>> [abs(eqs[-1].simplify().subs(g, i).n()) for i in gsol]
[0.e-145, 0.e-145, 7.84800000000000e-23]

You can check to see if they are also solutions of the others. (The real roots are not compact solutions in terms of radicals.)

You might leave the constants as symbols and repeat the above and only solve the last equation with values when needed -- if needed.

smichr
  • 16,948
  • 2
  • 27
  • 34