Solve a system of symbolic equations in python sympy

Question

Being new to sympy, I would benefit from a little help to put me on the rails... In the example below, how can I obtain the result for G_S from the system of two symbolic equations? The answer should be G_S = S11 - S12 S21 / (1 + S22).

S11,S12,S21,S22 = symbols("S11 S12 S21 S22")
S = Matrix([[S11, S12], [S21,S22]])
a1,b1,a2,b2 = symbols("a1 b1 a2 b2")
In = Matrix([a1, b2])
Out = Matrix([b1, a2])

system = [Eq(Out,S*In), Eq(a2+b2,0)]
G_S = b1/a1

Olivier

Answer 1

Based on the help given in the other answers, I summarize a satisfactory solution:

sol = solve(system, [a2,b2,b1])
G_S = sol[b1]/a1
G_S = factor_terms(collect(expand(G_S),'S11'))   # rearrange the expression

The Sij are viewed as coefficients while the ai and bi are viewed as variables. The system represents 3 equations and allows to solve for that number of unknow variables. We thus solve for a2,b2,b1 and divide the b1 solution by the a1 unknown.

Answer 2

Disclaimer: there is probably a more elegant way to do this. I'm pretty new to SymPy myself.

You just need to tell SymPy to actually solve the system of equations. Do all the steps just like you did, but right after you define system you need to tell SymPy to solve for b1 and a1, then use those solutions to calculate G_S. (Note that solve returns a dict.)

Edit: even though we're only using the answers for b1 and a1 when we calculate G_S, we still need to tell solve to solve for all four variables a1, b1, a2, b2 for it to give us the correct answer.

system = [Eq(Out,S*In), Eq(a2+b2,0)]
soln = solve(system, a1, b1, a2, b2)
G_S = soln[b1]/soln[a1]

When I do this SymPy gives me the correct-but-still-ugly answer of -(-S11*(S22 + 1) + S12*S21)/(S22 + 1). If I call simplify(G_S) I get the less-ugly-but-still-ugly answer of (S11*(S22 + 1) - S12*S21)/(S22 + 1). That's the fun of doing symbolic math with a computer -- its idea of "simple" is never quite the same as a human's.

Answer 3

@Mr-Snrub has shown how to solve the equations. To get a nicer simplification you can try

>>> eq
-(-S11*(S22 + 1) + S12*S21)/(S22 + 1)
>>> factor_terms(collect(expand(eq),'S11'))
S11 - S12*S21/(S22 + 1)

The problem with non-target simplification is that the sum ends up getting expanded and lost. By letting cse capture that repeated expression of S22 + 1 and expanding the cse expression, e, and backsubstituting the common expressions -- in this case, only 1 -- the expression you prefer is obtained:

>>> r,e=cse(eq)
>>> e[0].expand().subs(reversed(r))
S11 - S12*S21/(S22 + 1)