2

I'm trying to solve a system of 6 non linear equations with Python. So far I've got:

from scipy.optimize import fsolve 
import math 
from numpy import log

ka1 = 1.045
ka2 = 0.1759
ka3 = 0.3159

def equations(p):
    yh2o, yco2, yh2, ych4, yco = p

    return (yh2o + yco2 + yh2 + ych4 + yco - 1,
        ka1 - (yh2o ** 2)/(yco * (yh2 ** 2),
        ka2 - (yco ** 2) / yco2,
        ka3 - ych4 / (yh2 ** 2),
        0.5 - (2.0 * yco2 + yco + yh2o) / (2 * yh2o + 2 * yh2 + 4 * ych4)))

yh2o, yco2, yh2, ych4, yco = fsolve(equations, [0.2, 0.2, 0.2, 0.2, 0.2])

print(f"yh2o = {yh2o}, yco2 = {yco2}, yh2 = {yh2}, ych4 = {ych4}, yco = {yco}")

When I try to run this I get

TypeError: fsolve: there is a mismatch between the input and output shape of the 'func' argument 'equations'.Shape should be (5,) but it is (2,).

I've tried changing the initial guess but the error is still the same.

I've also tried changing the equations into polynomials like so:

    return (yh2o + yco2 + yh2 + ych4 + yco - 1.0,
        ka1 * (yco * (yh2 ** 2.0)) - (yh2o * yh2o),
        ka2 * (yco2) - (yco ** 2.0),
        ka3 * (yh2 ** 2.0) - ych4,
        0.5 * (2.0 * yh2o + 2.0 * yh2 + 4 * ych4) - (2.0 * yco2 + yco + yh2o)) 

But I get an error saying

TypeError: can't multiply sequence by non-int of type 'float'
  • I don't know about that specific error, but to make it easier to find a solution, you might consider clearing the fractions so you have polynomials instead of rational functions. I just mean to multiply through by the denominators -- note that puts implicit constraints on the results, namely that the roots of the denominators must be filtered out of any solutions. Equations in several variables are kind of touchy; maybe consider solving a special case, by assuming specific values for variables and solving for others, in order to get started. Good luck and have fun. – Robert Dodier Nov 03 '20 at 00:55
  • Thanks, I tried that but got the following error on line 13: `TypeError: can't multiply sequence by non-int of type 'float'` Not quite sure what this means – Carlos Importes Nov 03 '20 at 03:17
  • Thanks for the update. It's not clear what you tried. You can help others help you by updating your question to show exactly what you tried and the error you got. – Robert Dodier Nov 03 '20 at 03:24
  • Can you check [my solution](https://stackoverflow.com/a/64656440/941531), if I'm correct or wrong? – Arty Nov 03 '20 at 03:43
  • You are exactly right, I messed up a parenthesis. It's so frustrating to make these types of errors. Thank you! – Carlos Importes Nov 03 '20 at 03:51
  • @CarlosImportes It is very nice that still you posted a question like this. Because if you made such error it means that other people also may do in same cases, hence this question/answer will be useful for them! – Arty Nov 03 '20 at 03:54

1 Answers1

2

I think you just have problem with parentheses (i.e. placing of ( and )) in your equations (returned tuple).

Below is full corrected working code:

Try it online!

from scipy.optimize import fsolve 
import math 
from numpy import log

ka1 = 1.045
ka2 = 0.1759
ka3 = 0.3159

def equations(p):
    yh2o, yco2, yh2, ych4, yco = p

    return (
        (yh2o + yco2 + yh2 + ych4 + yco - 1),
        ka1 - (yh2o ** 2) / (yco * (yh2 ** 2)),
        ka2 - (yco ** 2) / yco2,
        ka3 - ych4 / (yh2 ** 2),
        0.5 - (2.0 * yco2 + yco + yh2o) / (2 * yh2o + 2 * yh2 + 4 * ych4)
    )

yh2o, yco2, yh2, ych4, yco = fsolve(equations, [0.2, 0.2, 0.2, 0.2, 0.2])

print(f"yh2o = {yh2o}, yco2 = {yco2}, yh2 = {yh2}, ych4 = {ych4}, yco = {yco}")

Output:

yh2o = 0.17829101808889491, yco2 = 0.17513942402710805, yh2 = 0.4163023675952086, ych4 = 0.05474789024926531, yco = 0.1755193000395232
Arty
  • 14,883
  • 6
  • 36
  • 69