CODEWITHSUNDEEP
pythonscipysympynon-linear
Thomas Murdoch
Thomas Murdoch

Reputation: 13

Why are some of my equations ignored when solving a non-linear system of equations

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

Upvotes: 1

Views: 124

Answers (1)

smichr
smichr

Reputation: 19135

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.

Upvotes: 1

Related Questions