Solving symbolic equations with SymPy

SymPy is a Python library for symbolic mathematics. It is one of the layers used in SageMath, the free open-source alternative to Maple/Mathematica/Matlab. When you have simple but big calculations that are tedious to be solved by hand, feed them to SymPy, and at least you can be sure it will make no calculation mistake ;-)

The basic functionalities of SymPy are expansion/factorization/simplification of symbolic expressions, limit calculations, differentiation, integration, algebraic equation solving, and some simple differential equation solving. See the SymPy tutorial by Fabian Pedregosa for a proper introduction.

Solving a system of quadratic equations

In the code snippet below, we want to solve a system of quadratic polynomials:

\begin{equation*} \begin{array}{rcl} z\,\ddot{x} + (x_Z - x)(g + \ddot{z}) & = & 0 \\ z\,\ddot{y} + (y_Z - y)(g + \ddot{z}) - \dot{L}_y & = & 0 \\ - y\,\ddot{x} + x\,\ddot{y} - z_Z\,(g + \ddot{z}) + \dot{L}_z & = & 0 \end{array} \end{equation*}

(The notations come from physics, where these equations are used to calculate the zero-moment point.) First, declare variables using the var() construct:

from sympy import var

Ldy, Ldz = var('Ldy Ldz')
g, x, y, z = var('g x y z')
xZ, yZ, zZ = var('xZ yZ zZ')
xdd, ydd, zdd = var('xdd ydd zdd')

You can then use them directly as Python variables, performing all common operations such as addition or multiplication. Next, define the expressions to be zeroed and pass them to the solve() function:

from sympy import solve

E1 = z * xdd + (xZ - x) * (g + zdd)
E2 = z * ydd + (yZ - y) * (g + zdd) - Ldy
E3 = -y * xdd + x * ydd - zZ * (g + zdd) + Ldz

sols = solve([E1, E2, E3], [xdd, ydd, Ldy])

print "xdd = ", (sols[xdd]).factor()
print "ydd = ", (sols[ydd]).factor()
print "Ldy = ", (sols[Ldy]).factor()

The second argument of solve() indicates the set of "output" variables. Indeed, we have three equations for twelve variables. Each equation can be used to express one variable as function of the others. Thus, we can pick three variables and express them as functions of the remaining nine. This is what we do here with \(\ddot{x}\), \(\ddot{y}\) and \(\dot{L}_y\).

Reading the solution

The output from this code is:

xdd =  (g + zdd)*(x - xZ)/z
ydd =  (-Ldz*z + g*x*y - g*xZ*y + g*z*zZ + x*y*zdd - xZ*y*zdd + z*zZ*zdd)/(x*z)
Ldy =  (-Ldz*z + g*x*yZ - g*xZ*y + g*z*zZ + x*yZ*zdd - xZ*y*zdd + z*zZ*zdd)/x

Which corresponds to the solution:

\begin{equation*} \begin{array}{rcl} \ddot{x} & = & \frac{1}{z} \left(g + \ddot{z}\right) \left(x - x_Z\right) \\ \ddot{y} & = & \frac{1}{xz} \left(- \dot{L}_z z + g x y - g x_Z y + g z z_Z + x y \ddot{z} - x_Z y \ddot{z} + z z_Z \ddot{z}\right) \\ \dot{L}_y & = & \frac{1}{x} \left(- \dot{L}_z z + g x y_Z - g x_Z y + g z z_Z + x y_Z \ddot{z} - x_Z y \ddot{z} + z z_Z \ddot{z}\right) \end{array} \end{equation*}

SymPy has some routines to make formulas more palatable. For instance, it can print sympy.Expr objects (expressions) in \(\LaTeX\):

from sympy import latex

print "\\begin{array}{rcl}"
print "xdd & = & %s \\\\" % latex((sols[xdd]).factor())
print "ydd & = & %s \\\\" % latex((sols[ydd]).factor())
print "Ldy & = & %s" % latex((sols[Ldy]).factor())
print "\\end{array}"

If you use IPython's QTConsole, you can even render \(\LaTeX\) formulas directly in your console. See Printing from the SymPy documentation for details.

Pages of this website are under the CC-BY 4.0 license.