sympy.solve(), m equations, n unknowns, n > m, find one unknown expressed in terms of n - m - 1 other unknowns

Question

This is a question on using the Sympy solve() function in Python.

I am limited by my understanding of all the options for sympy.solve(). I am trying to use the symbolic solver to do the following:

Suppose I have 5 equations in 8 variables (unknowns). I would like to find an expression for one variable in terms of 2 other variables (using the 5 equations to eliminate 5 of the variables). I.e., given:

f1(x,a,b,c,d,e,u,w) = 0
f2(x,a,b,c,d,e,u,w) = 0
f3(x,a,b,c,d,e,u,w) = 0
f4(x,a,b,c,d,e,u,w) = 0
f5(x,a,b,c,d,e,u,w) = 0

find x = g(u,w) as an expression, if possible. In fact, I would be happy with an expression of the form: g(u,w,x) = 0, if Sympy can basically eliminate the other variables and return an expression that could be solved numerically for a variable when two others are specified.

My functions are a mix of linear and non-linear equations, fairly simple ones. I believe there is a symbolic solution based on manual attempts to solve it, but the substitutions became too unwieldy to complete accurately. I think Sympy.solve() should be able to do what I want, but I'm not sure how to set it up. Is there a way to do this? Also, is there a "verbose" mode for Sympy.solve() to understand why it returns nothing sometimes (not even empty list).

Answer 1

This is not an answer, rather a version of the original question. (I cannot find a way to format comments).

Here are the actual equations I wish to manipulate with Sympy:

z1 = (ha*hp)/hha - ka1
z2 = (a*hp)/ha - ka2
z3 = oh*hp - kw
z4 = hha + ha + a - I
z5 = oh + ha + 2*a - hp

I believe it is possible to derive an equation with just hp, ka1, ka2, kw, and "I" variables remaining. Ka1, ka2, kw, and "I" eventually can be given numeric values, resulting in an expression just in hp, requiring graphic/numeric solution. If I try this:

smp.solve([z1,z2,z3,z4,z5],set=True,manual=True,warn=True,check=False)

I get nothing, but I don't know the right way to set it up. My goal is to obtain something of the form:

g(hp, ka1, ka2, kw, I) = 0

and ultimately:

g(hp) = 0

after substituting numbers for the other 4 variables.

Occasionally Sympy does spit out "solutions". When I substituted numbers for the 4 variables mentioned, a solution appeared once, but I could not repeat it. Perhaps some small change caused it, I can't reliably return it to the state where it gave some answers.

Answer 2

If you solve((z1,z2,z3,z4),(a,ha,hha,oh)) and substitute into z5 and re-arrange for the numerator that should be 0, you end up with a quartic in hp:

>>> sol = solve((z1,z2,z3,z4),(a,ha,hha,oh))
>>> collect(z5.subs(sol).as_numer_denom()[0], hp)
-hp**4 - hp**3*ka1 + hp**2*(I*ka1 - ka1*ka2 + kw) + hp*(2*I*ka1*ka2 + ka1*kw) + ka1*ka2*kw

There is not a simple form for the symbolic solution of the quartic, but if you plug in the values for the other variables you will find a solution, preferably with real_root or nsolve.

I think the basic approach, on sets of equations like this, is to not get too far ahead of yourself in terms of requesting a solution and identifying what you do and don't want to solve for. You are interested in ha so you will need one equation for that -- and maybe you don't know if it will be complicated or not -- so solve 4 of the equations for symbols other than ha and other input parameters like ka1, ka2, and kw.

from sympy import *
...
eqs = (z1,z2,z3,z4,z5)
for s in subsets(eqs,4):
 sol = solve(s, exclude=(ha,kw,ka1,ka2), dict=True)
 print([eqs.index(i) for i in s])
 e = [i for i in eqs if i not in s][0]
 from sympy.solvers.solvers import unrad
 for si in sol:
        n = e.subs(si).as_numer_denom()[0].expand()
        u = unrad(n, ha)
        if u is not None:
            n = collect(u[0], ha)
        print('\t',degree(Poly(n,ha)),list(ordered(n.free_symbols-{ha})))

This gives a list of the equations selected and the order of ha in the final equation after substituting in the answer obtained from the 4 equations (potentially needing to use unrad to remove a square root), and the list of free symbols other than ha in that equation.

[0, 1, 2, 3]
     1 [ka2, kw, oh]
[0, 1, 2, 4]
     5 [I, ka1, ka2, kw]
     5 [I, ka1, ka2, kw]
[0, 1, 3, 4]
     1 [ka2, kw, oh]
     1 [ka2, kw, oh]
[0, 2, 3, 4]
     1 [ka2, kw, oh]
[1, 2, 3, 4]
     3 [hha, ka1, ka2, kw]
     3 [hha, ka1, ka2, kw]
​```

So it's up to you to pick the 4 you want to use to get to an equation with `ha` and the desired parameters that you want to be able to use as input variables. The Groebner approach does this picking on its own based on the variables that you pass to it.

Answer 3

You can do this with a Groebner basis. Here are your equations:

In [22]: from sympy import *

In [23]: a, ha, hp, hha, ka1, ka2, oh, I = symbols('a, ha, hp, hha, ka1, ka2, oh, I')

In [24]: z1 = (ha*hp)/hha - ka1
    ...: z2 = (a*hp)/ha - ka2
    ...: z3 = oh*hp - kw
    ...: z4 = hha + ha + a - I
    ...: z5 = oh + ha + 2*a - hp

In [25]: eqs = [z1, z2, z3, z4, z5]

In [26]: eqs
Out[26]: 
⎡ha⋅hp        a⋅hp                                                         ⎤
⎢───── - ka₁, ──── - ka₂, hp⋅oh - kw, -I + a + ha + hha, 2⋅a + ha - hp + oh⎥
⎣ hha          ha                                                          ⎦

This is almost a polynomial system but not quite so let's convert to polynomial:

In [27]: eqs = [eq.as_numer_denom()[0] for eq in eqs]

In [28]: eqs
Out[28]: [ha⋅hp - hha⋅ka₁, a⋅hp - ha⋅ka₂, hp⋅oh - kw, -I + a + ha + hha, 2⋅a + ha - hp + oh]

Now compute a Groebner basis specifying the symbols you want to eliminate and the symbol that you want to keep:

In [29]: gb = groebner(eqs, [hha, ha, oh, a, hp])

In [30]: for eq in gb: pprint(eq)
    ⎛   2            ⎞       3     2                                    
hha⋅⎝ka₁  - 4⋅ka₁⋅ka₂⎠ + 2⋅hp  + hp ⋅ka₁ + hp⋅(-2⋅I⋅ka₁ - 2⋅kw) - ka₁⋅kw
   ⎛               2⎞     3     2                                                     
ha⋅⎝ka₁⋅ka₂ - 4⋅ka₂ ⎠ - hp  + hp ⋅(-ka₁ + 2⋅ka₂) + hp⋅(I⋅ka₁ + kw) + ka₁⋅kw - 2⋅ka₂⋅kw
                3     2                                                      
2⋅I⋅ka₁⋅ka₂ - hp  - hp ⋅ka₁ + hp⋅(I⋅ka₁ - ka₁⋅ka₂ + kw) + ka₁⋅ka₂⋅oh + ka₁⋅kw
       2                  2     ⎛   2                2⎞     3                   2 ⎛   2            
- I⋅ka₁ ⋅ka₂ + 4⋅I⋅ka₁⋅ka₂  + a⋅⎝ka₁ ⋅ka₂ - 4⋅ka₁⋅ka₂ ⎠ + hp ⋅(ka₁ - 2⋅ka₂) + hp ⋅⎝ka₁  - 3⋅ka₁⋅ka₂

⎞      ⎛       2                                  ⎞      2                  
⎠ + hp⋅⎝- I⋅ka₁  + 2⋅I⋅ka₁⋅ka₂ - ka₁⋅kw + 2⋅ka₂⋅kw⎠ - ka₁ ⋅kw + 3⋅ka₁⋅ka₂⋅kw
  4     3         2                                                                  
hp  + hp ⋅ka₁ + hp ⋅(-I⋅ka₁ + ka₁⋅ka₂ - kw) + hp⋅(-2⋅I⋅ka₁⋅ka₂ - ka₁⋅kw) - ka₁⋅ka₂⋅kw

In [31]: gb[-1]
Out[31]: 
  4     3         2                                                                  
hp  + hp ⋅ka₁ + hp ⋅(-I⋅ka₁ + ka₁⋅ka₂ - kw) + hp⋅(-2⋅I⋅ka₁⋅ka₂ - ka₁⋅kw) - ka₁⋅ka₂⋅kw

The final equation gb[-1] has the other symbols eliminated and gives a polynomial for hp where the coefficients depend on the other variables listed as parameters. It is possible to get an explicit formula for hp since this is a quartic but the formula will be messy. Instead as you suggest we can substitute numeric values and solve:

In [34]: gb[-1].subs({ka1:1,ka2:2,I:3,kw:4})
Out[34]: 
  4     3       2            
hp  + hp  - 5⋅hp  - 16⋅hp - 8

In [35]: nroots(gb[-1].subs({ka1:1,ka2:2,I:3,kw:4}))
Out[35]: 
[-0.629686400933711, 2.91693500045315, -1.64362429975972 - 1.28608250861107⋅ⅈ, -1.64362429975972 + 
1.28608250861107⋅ⅈ]

In [36]: [r.n(3) for r in real_roots(gb[-1].subs({ka1:1,ka2:2,I:3,kw:4}))]
Out[36]: [-0.63, 2.92]

You could also use lambdify and numpy's roots function to get faster but less accurate roots.