SymPy Examples

sources: http://docs.sympy.org/latest/tutorial/ http://nbviewer.ipython.org/github/ipython/ipython/blob/master/examples/notebooks/SymPy Examples.ipynb

SymPy provides support for symbolic math to python, similar to what you would do with Mathematica or Maple. The major difference is that it acts just like any other python module, so you can use the symbolic math together in your own python projects with the rest of python functionality.

The following import and function (init_session()) sets up a nice environment for us when working in Jupyter

from sympy import init_session
init_session(use_latex="mathjax")

IPython console for SymPy 1.13.1 (Python 3.11.9-64-bit) (ground types: python)

These commands were executed:
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at https://docs.sympy.org/1.13.1/

import math

SymPy types and basic symbolic manipulation

Sympy defines its own types, you can convert them to python types, but you don’t always want to (and will probably lose accuracy when you do).

print(math.sqrt(2))

1.4142135623730951

print(sqrt(2))

sqrt(2)

print(sqrt(8))

2*sqrt(2)

#help(sqrt(8))

We can do symbolic math not just on numbers, but we can tell SymPy what to treat as a symbol, using symbols()

from sympy import symbols
x, y, z = symbols("x y z")

expr = x + 2*y
expr

\[\displaystyle x + 2 y\]

expr - 1

\[\displaystyle x + 2 y - 1\]

expr - y

\[\displaystyle x + y\]

f = x*expr
f

\[\displaystyle x \left(x + 2 y\right)\]

g = expand(f)
g

\[\displaystyle x^{2} + 2 x y\]

factor(g)

\[\displaystyle x \left(x + 2 y\right)\]

substitution

SymPy provides methods to substitute values for symbols in symbolic expressions. Note, the follow likely does not do what you expect:

expr = sin(z*2*pi)
z = 0
expr

\[\displaystyle \sin{\left(2 \pi z \right)}\]

We’ve now redefined z to be a python type

type(z)

int

to do substitution, we use the subs() method

expr = sin(x*2*pi)
expr

\[\displaystyle \sin{\left(2 \pi x \right)}\]

a = expr.subs(x, 0.125)
a

\[\displaystyle \frac{\sqrt{2}}{2}\]

Note that this is not a floating point number – it is still a SymPy object. To make it floating point, we can use evalf()

b = a.evalf()
print(b, type(b))

0.707106781186548 <class 'sympy.core.numbers.Float'>

This is still a SymPy object, because SymPy can do arbitrary precision

a.evalf(50)

\[\displaystyle 0.70710678118654752440084436210484903928483593768847\]

want regular python types?

c = float(b)
print(c, type(c))

0.7071067811865476 <class 'float'>

Python and SymPy

x, y, z, t = symbols('x y z t')

SymPy symbols are just objects and when you do operations on two sympy objects the result is a sympy object.

When you combine a sympy and python object, the result is also a sympy object.

But we need to be careful when doing fractions. For instance doing x + 1/3 will first compute 1/3 in python (giving 0.333...) and then add it to the sympy x symbol. The Rational() function makes this all happen in sympy

f = expr + Rational(1,3)
f

\[\displaystyle \sin{\left(2 \pi x \right)} + \frac{1}{3}\]

expr + 1/3

\[\displaystyle \sin{\left(2 \pi x \right)} + 0.333333333333333\]

equality

= is still the assignment operator of python (it does not mean symbolic equality), and == is still the logical test (exact structural equality). There is a separate object, Eq() to specify symbolic equality.

And testing for algebraic equality is not always accomplished using ==, since that tests for structural equality.

x + 1 == 4

False

Eq(x + 1, 4)

\[\displaystyle x + 1 = 4\]

a = (x + 1)**2
b = x**2 + 2*x + 1    # these are algebraically equal

a == b

False

We can use simplify() to test for algebraic equality

simplify(a - b)

\[\displaystyle 0\]

a = cos(x) + I*sin(x)
a

\[\displaystyle i \sin{\left(x \right)} + \cos{\left(x \right)}\]

simplify(a)

\[\displaystyle e^{i x}\]

More substitution

note that substitution returns a new expression: SymPy expressions are immutable

expr = cos(x)
expr.subs(x, 0)

\[\displaystyle 1\]

expr

\[\displaystyle \cos{\left(x \right)}\]

x

\[\displaystyle x\]

multiple substitutions, pass a list of tuples

expr = x**3 + 4*x*y - z
expr

\[\displaystyle x^{3} + 4 x y - z\]

expr.subs([(x, 2), (y, 4), (z, 0)])

\[\displaystyle 40\]

simplifying

There is not unique definition of what the simplest form of an expression is.

simplify() tries lots of methods for simplification

simplify(sin(x)**2 + cos(x)**2)

\[\displaystyle 1\]

simplify( (x**3 + x**2 - x - 1)/(x**2 + 2*x + 1) )

\[\displaystyle x - 1\]

simplify(gamma(x)/gamma(x - 2))

\[\displaystyle \left(x - 2\right) \left(x - 1\right)\]

but sometimes it doesn’t have your idea of what the simplest form is

simplify(x**2 + 2*x + 1)

\[\displaystyle x^{2} + 2 x + 1\]

instead factor may be what you want

factor(x**2 + 2*x + 1)

\[\displaystyle \left(x + 1\right)^{2}\]

polynomial simplification

expand((x + 1)**2)

\[\displaystyle x^{2} + 2 x + 1\]

expand((x + 2)*(x - 3))

\[\displaystyle x^{2} - x - 6\]

expand( (x + 1)*(x - 2) - (x - 1)*x)

\[\displaystyle -2\]

factor(x**2*z + 4*x*y*z + 4*y**2*z)

\[\displaystyle z \left(x + 2 y\right)^{2}\]

factor_list(x**2*z + 4*x*y*z + 4*y**2*z)

\[\displaystyle \left( 1, \ \left[ \left( z, \ 1\right), \ \left( x + 2 y, \ 2\right)\right]\right)\]

collect collects common powers

expr = x*y + x - 3 + 2*x**2 - z*x**2 + x**3
expr

\[\displaystyle x^{3} - x^{2} z + 2 x^{2} + x y + x - 3\]

collected_expr = collect(expr, x)
collected_expr

\[\displaystyle x^{3} + x^{2} \left(2 - z\right) + x \left(y + 1\right) - 3\]

cancel cancels

a = (x**2 + 2*x + 1)/(x**2 + x)
a

\[\displaystyle \frac{x^{2} + 2 x + 1}{x^{2} + x}\]

cancel(a)

\[\displaystyle \frac{x + 1}{x}\]

trigsimp simplifies trigonometric identities

trigsimp(sin(x)**4 - 2*cos(x)**2*sin(x)**2 + cos(x)**4)

\[\displaystyle \frac{\cos{\left(4 x \right)}}{2} + \frac{1}{2}\]

trigsimp(sin(x)*tan(x)/sec(x))

\[\displaystyle \sin^{2}{\left(x \right)}\]

the tutorial discusses some of the nuances of simplification of powers and special functions

Calculus

Calculus operations are simple in SymPy

derivatives

diff(cos(x), x)

\[\displaystyle - \sin{\left(x \right)}\]

diff(exp(x**2), x)

\[\displaystyle 2 x e^{x^{2}}\]

third derivative

diff(x**4, x, 3)

\[\displaystyle 24 x\]

differentiate different variables

expr = exp(x*y*z)
diff(expr, x, y, z)

\[\displaystyle \left(x^{2} y^{2} z^{2} + 3 x y z + 1\right) e^{x y z}\]

unevaluated derivatives can be useful for building up ODEs and PDEs

deriv = Derivative(expr, x, y, z)
deriv

\[\displaystyle \frac{\partial^{3}}{\partial z\partial y\partial x} e^{x y z}\]

deriv.doit()

\[\displaystyle \left(x^{2} y^{2} z^{2} + 3 x y z + 1\right) e^{x y z}\]

integrals

definite and indefinite integrals are supported

integrate(cos(x), x)

\[\displaystyle \sin{\left(x \right)}\]

definite integral – note the construction of the infinity

integrate(exp(-x), (x, 0, oo))

\[\displaystyle 1\]

double integral

integrate(exp(-x**2 - y**2), (x, -oo, oo), (y, -oo, oo))

\[\displaystyle \pi\]

if it is unable to do the integral, it returns an Integral object

expr = integrate(x**x, x)
print(expr)
expr

Integral(x**x, x)

\[\displaystyle \int x^{x}\, dx\]

a = x / sqrt(x**4 + 10*x**2 - 96*x - 71)   # example from Wikipedia Risch algorithm page)
a

\[\displaystyle \frac{x}{\sqrt{x^{4} + 10 x^{2} - 96 x - 71}}\]

integrate(a, x)     # this has a known solution, but SymPy fails to find it

\[\displaystyle \int \frac{x}{\sqrt{x^{4} + 10 x^{2} - 96 x - 71}}\, dx\]

limits

limit(sin(x)/x, x, 0)

\[\displaystyle 1\]

series expansions

expr = exp(sin(x))
a = expr.series(x, 0, 10)
a

\[\displaystyle 1 + x + \frac{x^{2}}{2} - \frac{x^{4}}{8} - \frac{x^{5}}{15} - \frac{x^{6}}{240} + \frac{x^{7}}{90} + \frac{31 x^{8}}{5760} + \frac{x^{9}}{5670} + O\left(x^{10}\right)\]

c = log(x).series(x, x0=1, n=6)
c

\[\displaystyle -1 - \frac{\left(x - 1\right)^{2}}{2} + \frac{\left(x - 1\right)^{3}}{3} - \frac{\left(x - 1\right)^{4}}{4} + \frac{\left(x - 1\right)^{5}}{5} + x + O\left(\left(x - 1\right)^{6}; x\rightarrow 1\right)\]

simplify(c.removeO())

\[\displaystyle \frac{x^{5}}{5} - \frac{5 x^{4}}{4} + \frac{10 x^{3}}{3} - 5 x^{2} + 5 x - \frac{137}{60}\]

solvers

solveset() is the main interface to solvers in SymPy. Note that it used to be solve(), but this has been replaced (see http://docs.sympy.org/latest/modules/solvers/solveset.html)

If no Eq() is done, then it is assumed to be equal to 0

solveset(x**2 - x, x)

\[\displaystyle \left\{0, 1\right\}\]

you can restrict the domain of the solution (e.g. to reals). Recall that Z is the set of integers

solveset(sin(x) - 1, x, domain=S.Reals)

\[\displaystyle \left\{2 n \pi + \frac{\pi}{2}\; \middle|\; n \in \mathbb{Z}\right\}\]

linear systems

linsolve() is the interface to linear systems

linsolve([x - y + 2, x + y - 3], [x, y])

\[\displaystyle \left\{\left( \frac{1}{2}, \ \frac{5}{2}\right)\right\}\]

linsolve([x + y + z - 1, x + y + 2*z - 3 ], (x, y, z))

\[\displaystyle \left\{\left( - y - 1, \ y, \ 2\right)\right\}\]

roots will report if a solution is multiple by listing it multiple times

roots(x**3 - 6*x**2 + 9*x, x)

\[\displaystyle \left\{ 0 : 1, \ 3 : 2\right\}\]

0 is 1 root, and 3 is 2 more roots

Differential equations

you need an undefined function (f and g already are by our init_session() above, but we’ve probably reset these

f, g = symbols('f g', cls=Function)

f(x)

\[\displaystyle f{\left(x \right)}\]

f(x).diff(x)

\[\displaystyle \frac{d}{d x} f{\left(x \right)}\]

diffeq = Eq(f(x).diff(x, 2) - 2*f(x).diff(x) + f(x), sin(x))

diffeq

\[\displaystyle f{\left(x \right)} - 2 \frac{d}{d x} f{\left(x \right)} + \frac{d^{2}}{d x^{2}} f{\left(x \right)} = \sin{\left(x \right)}\]

dsolve(diffeq, f(x))

\[\displaystyle f{\left(x \right)} = \left(C_{1} + C_{2} x\right) e^{x} + \frac{\cos{\left(x \right)}}{2}\]

Matrices

consider the Euler equations:

\[q_t + A(q) q_x = 0\]

where

\[\begin{split}q = \left ( \begin{array}{c} \rho \\ u \\ p \end{array} \right ) \qquad A(q) = \left ( \begin{array}{ccc} u & \rho & 0 \\ 0 & u & 1/\rho \\ 0 & c^2 \rho & u \end{array} \right ) \end{split}\]

from sympy.abc import rho
rho, u, c = symbols('rho u c')
A = Matrix([[u, rho, 0], [0, u, rho**-1], [0, c**2 * rho, u]])
A

\[\begin{split}\displaystyle \left[\begin{matrix}u & \rho & 0\\0 & u & \frac{1}{\rho}\\0 & c^{2} \rho & u\end{matrix}\right]\end{split}\]

A.row(0)

\[\displaystyle \left[\begin{matrix}u & \rho & 0\end{matrix}\right]\]

The eigenvalues of the system are the speeds at which information propagates

A.eigenvals()

\[\displaystyle \left\{ u : 1, \ - c + u : 1, \ c + u : 1\right\}\]

You can diagonalize it, such that $\( A = PDP^{-1}\)$

P, D = A.diagonalize()

\(D\) will be a matrix of the eigenvalues

D

\[\begin{split}\displaystyle \left[\begin{matrix}u & 0 & 0\\0 & - c + u & 0\\0 & 0 & c + u\end{matrix}\right]\end{split}\]

\(P\) will be the matrix of right eigenvectors

P

\[\begin{split}\displaystyle \left[\begin{matrix}1 & \frac{1}{c^{2}} & \frac{1}{c^{2}}\\0 & - \frac{1}{c \rho} & \frac{1}{c \rho}\\0 & 1 & 1\end{matrix}\right]\end{split}\]

Inverse

A**-1

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{1}{u} & - \frac{\rho}{- c^{2} + u^{2}} & \frac{1}{- c^{2} u + u^{3}}\\0 & \frac{u}{- c^{2} + u^{2}} & - \frac{1}{- c^{2} \rho + \rho u^{2}}\\0 & - \frac{c^{2} \rho}{- c^{2} + u^{2}} & \frac{u}{- c^{2} + u^{2}}\end{matrix}\right]\end{split}\]

Units

Sympy can attach units to numbers and propagate them through

from sympy.physics.units import newton, kilogram, meter, second, convert_to

F = 1 * kilogram * 9.81 * meter / second**2
convert_to(F, newton)

\[\displaystyle 9.81 \text{N}\]