How to define and solve mutually referencing expressions in SymPy

Question

Suppose I have the following system of equations:

x = 1/y + b
y = 2*x

It follows that b = x-1/(2*x)

If I try to define these expressions in SymPy though,

from sympy import *
b = symbols('b')
x = 1/y + b
y = 2*x

it gives me

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-4-5a9b5f6b923e> in <module>()
      1 from sympy import *
      2 b = symbols('b')
----> 3 x = 1/y + b
      4 y = 2*x
      5 

NameError: name 'y' is not defined

What's the proper way to set up this set of equations so that I can solve it for x? Extending that, how to do an evalf on it?

Answer 1

You can use solve to solve your equations.

In [1]: from sympy import *

In [2]: from sympy.abc import x, y, b

In [3]: eq1 = Eq(x - 1/y - b, 0)

In [4]: eq2 = Eq(2*x - y, 0)

In [5]: solve([eq1, eq2], b)
Out[5]: {b: x - 1/y}

Answer 2

= expression is not does what you think. In regular python syntax the = operator used for variable assignment. You should define everything in symbols then convert equation to be equal to zero and solve it. The main problem is the lack of general procedural programming knowledge.

You should use solvers which requires the equation to be expressed as the following.

x = 1/y + b
x - 1/y - b = 0

The code could be written as following.

from sympy.solvers import solve
from sympy import symbols

b = symbols('b')
x = symbols('x')
y = 2 * x

solve(x - 1/y - b, b)

Then the output will be what you want to see.

>>> [x - 1/(2*x)]

If you want to express b in terms of y the following code could be used.

from sympy.solvers import solve
from sympy import symbols

b = symbols('b')
y = symbols('y')
x = y/2

solve(x - 1/y - b, b)

Then the output will be:

>>> [y/2 - 1/y]

If you want to express b in terms of both x and y you could use the following (thanks to Sudhanshu Mishra's answer):

from sympy.solvers import solve
from sympy import symbols, Eq

b = symbols('b')
y = symbols('y')
x = symbols('x')

eq1 = Eq(1/y + b, x)
eq2 = Eq(2 * x, y)

solve([eq1, eq2], b)

The result will be:

>>> {b: x - 1/y}