SymPy#

Up until this point, we have worked with “non-symbolic” Python features. For example, if we set \(x = 1\) and then \(y = x\), the apparent relationship between \(x\) and \(y\) disappears when we change \(x\).

x = 1
y = x
x = 2
y == x

False

This is because we are working with values for \(x\) and \(y\), rather than the symbols \(x\) and \(y\) themselves - i.e. we are doing arithmetic rather than algebra. The module SymPy is Python’s Computer Algebra System (CAS). You might have used other CASs like Maple or Mathematica; SymPy has the major advantage of being free and open source. It has less features than the commercial packages, but is steadily growing.

First we will import sympy, and use init_printing to get SymPy to output maths nicely using LaTeX - this is optional.

import sympy as sym
sym.init_printing()

Creating symbols#

To create “symbols” to manipulate, use sym.symbols. Note that it takes a tuple of strings.

x, y, z = sym.symbols(("x", "y", "z"))

It’s worth noting that your variable names can be different from the names that SymPy stores for the symbols:

# This is a recipe for confusion - you better have a good reason to do it!
u, v = sym.symbols(("s", "t"))

# This is a more reasonable thing to do
x0 = sym.symbols(("x_0"))
x0
../_images/3Sympy_7_0.png

If you know that you are only interested in real values of the variables, use real=True when creating it. There are other assumptions you can put on symbols (like being positive) - see here for a guide.

s = sym.symbols("s", real=True)

Expressions and substitution#

We can use symbols to create expressions, mostly using the same syntax as for regular Python variables:

(3 * (x + y)**2 * sym.exp(x*y))/sym.log(z)
../_images/3Sympy_11_0.png

Note that we used the SymPy version of exp and log, rather than the numpy version. This is important - the numpy version simply will not work. This is because numpy.exp expects a number (or an array of numbers), rather than a symbol. There are SymPy equivalents of most of the standard mathematical functions you would use in numpy.

Let’s assign this expression to a variable:

V = (3 * (x + y)**2 * sym.exp(x*y))/sym.log(z)

V
../_images/3Sympy_14_0.png

V looks like a function of three variables, but it is not a Python function, which is why we cannot “call” it with three arguments:

V(0, 1, 2)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [8], in <cell line: 1>()
----> 1 V(0, 1, 2)

TypeError: 'Mul' object is not callable

Instead, we should “substitute” values for \(x\), \(y\), and \(z\) like so:

V.subs(x, 2)
../_images/3Sympy_18_0.png

Note that this has not changed V itself - it has simply returned a new expression where \(x\) has been replaced with \(2\).

V
../_images/3Sympy_20_0.png

To substitute multiple values, pass a list of (symbol, replacement) pairs:

V.subs([(x, 2), (y, 3), (z, 2)])
../_images/3Sympy_22_0.png

We are not restricted to substituting in numbers. Suppose that we have a parameterisation \(x = \sin(u)\), \(y = \cos(u)\) for \(u \in [0, 2\pi)\). We could make these substitutions like so:

u = sym.symbols("u")

V.subs([(x, sym.sin(u)),(y, sym.cos(u))])
../_images/3Sympy_24_0.png

Evaluation#

Note that the expression above is still being stored as a symbolic expression rather than as a float, which helps to avoid rounding errors and the other perils associated with floats. If you really want a float, use evalf, providing it with a dict of symbol : value pairs.

# The first argument to evalf controls how precise the output is.
# The default is 15 places.

V.evalf(30, subs={x: 2, y: 3, z: 2})
../_images/3Sympy_26_0.png

Exact values#

When creating expressions, we want to avoid floats when possible. In the following example, we try to create the exact expression \(\frac{8}{3}x\) but instead get a rounded float coefficient. This is because the fraction \(8/3\) is evaluated by Python before SymPy can get involved.

8/3 * x
../_images/3Sympy_28_0.png

There are several ways of avoiding this: you can create the fraction as a Rational number:

from sympy import Rational

u = Rational(8, 3) * x
u
../_images/3Sympy_30_0.png

Alternatively, you can “split up” your fraction like so:

u = 8 * x / 3
u
../_images/3Sympy_32_0.png

The same is true of roots: instead of doing

x ** (1/3)
../_images/3Sympy_34_0.png

You should do

sym.root(x, 3)
../_images/3Sympy_36_0.png

Manipulating expressions#

SymPy comes with the ability to simplify, factorise, and expand expressions. As an example, we will work with the cubic \(p = (x-2)(x + 3)(x + 7)\)

P = (x - 2) * (x + 3) * (x + 7)
P
../_images/3Sympy_39_0.png 
Q = P.expand()
Q
../_images/3Sympy_40_0.png

To go the other way, use factor:

Q.factor()
../_images/3Sympy_42_0.png

An interesting feature is that P == Q evaluates to False. This is because they are not “structurally” the same expression, even if they are the same “algebraic” expression.

P == Q

False

The same is true of P - Q == 0:

P - Q == 0

False

In order to check whether P and Q represent the same expression, the simplest thing to do is to simplify the expression P - Q and see whether it is \(0\).

(P - Q).simplify() == 0

True

Note that in general, it is not clear what is meant by a “simpler” expression, so simplify may return results you do not expect. It is also theoretically impossible for simplify to always work out whether an expression is equal to 0, though it will usually manage it for most reasonable expressions.

Solving equations#

One variable#

To find the roots of an expression, use solveset. For example, to find the roots of the polynomial \(Q(x) = x^2 - x - 6\):

Q = x**2 - x - 6

sym.solveset(Q)
../_images/3Sympy_52_0.png

If we don’t want roots (i.e. solutions to \(Q(x) = 0\)) but instead solutions to \(Q(x) = a\) for some \(a \in \mathbb{R}\), simply find the roots of \(Q - a\). For example, to get the solutions of \(Q(x) = 1\) do:

sym.solveset(Q - 1)
../_images/3Sympy_54_0.png

If we’re only interested in the real solutions, we can use domain=sym.Reals:

sym.solveset(Q + 7, domain = sym.Reals)
../_images/3Sympy_56_0.png

Note that solveset does not give us the multiplicity of roots:

sym.solveset(x**2 + 2*x + 1)
../_images/3Sympy_58_0.png

For polynomials only, we can get the multiplicities using roots:

sym.roots(x**2 + 2*x + 1)
../_images/3Sympy_60_0.png

You can also obtain the roots of an expression in terms of another variable:

# solve x^2 + yx + 1 = 0 for x
sym.solveset(x**2 + y*x + 1, x)
../_images/3Sympy_62_0.png

To access these solutions, you can convert them to a list:

sols = sym.solveset(x**2 + y*x + 1, x)

sols = list(sols)
sols[0]
../_images/3Sympy_65_0.png

Alternatively, you can iterate through them directly:

sols = sym.solveset(x**2 + y*x + 1, x)
for sol in sols:
    print(f"{sol} evaluated at y = 1: {sol.evalf(subs={y: 1})}")

-y/2 - sqrt((y - 2)*(y + 2))/2 evaluated at y = 1: -0.5 - 0.866025403784439*I
-y/2 + sqrt((y - 2)*(y + 2))/2 evaluated at y = 1: -0.5 + 0.866025403784439*I

Solving systems of equations#

Suppose we want to find stationary points of the following system of ODES

(1)#\[\begin{align} \frac{du}{dt} &= (u+1)v \\ \frac{dv}{dt} &= v(u + 2v). \end{align}\]

That is, we want to solve the system

(2)#\[\begin{align} (u+1)v &= 0 \\ v(u + 2v) &= 0. \end{align}\]

First we create our symbols. We will use real=True to only get real solutions to our system.

u, v = sym.symbols("u v", real=True)

Now our system is represented by the following expressions (which SymPy will assume we are trying to find roots of, as above)

eqs = [(u + 1) * v, v*(u + 2*v)]

And we can solve like so:

sols = sym.solve(eqs, [u, v])

sols
../_images/3Sympy_74_0.png

Note that unlike solveset, the output from solve is already in list form.

sols[1][1]
../_images/3Sympy_76_0.png

Differentiation#

One of the nice features of SymPy is its ability to perform differentiation (both partial and total). Let’s consider a moderately complicated function with a high likelihood of mistakes if we took derivatives by hand, \(V(x) = (\tan^2(x) + 1)^3\sinh(x^3)\):

V = (sym.tan(x) ** 2 + 1)**3 * sym.sinh(x**3)

To differentiate, use sym.diff, specifying which variable to differentiate with respect to (here there is of course only one option):

# Differentiate V with respect to x
sym.diff(V, x)
../_images/3Sympy_80_0.png

To take higher derivatives, just add a number to the end:

# d^2/dx^2 (V):
sym.diff(V, x, 2)
../_images/3Sympy_82_0.png

Alternatively, you can write out the variable you are differentiating with respect to repeatedly:

sym.diff(V, x, x)
../_images/3Sympy_84_0.png

Partial Derivatives#

Now let’s consider a function of three variables \(x, y, z\):

V = (3 * (x + y)**2 * sym.exp(x*y))/sym.log(z)

V
../_images/3Sympy_86_0.png

We can compute \(\frac{\partial V}{\partial x}\) in precisely the same way as above:

sym.diff(V, x)
../_images/3Sympy_88_0.png

The same would work for any of the other first partial derivatives; just replace x with the variable of interest.

If we instead wanted \(\frac{\partial^2 V}{\partial z\partial x}\), we could obtain it like so:

sym.diff(V, x, z)
../_images/3Sympy_90_0.png

In general, to compute \(\frac{\partial^N V}{\partial x_1\partial x_2 \cdots \partial x_N}\), use sym.diff(V, x_N, x_(N-1), ..., x_2, x_1). Remember that the order can matter (though very often does not for the functions we are interested in).

Unevaluated derivative#

It is possible to construct an object representing a derivative without actually computing it:

D = sym.Derivative(V, x, y)
D
../_images/3Sympy_93_0.png

This might seem like an odd thing to do, but you might want to build up an expression involving a very messy derivative without having it expand and make it difficult to see what’s going on. To evaluate such a derivative, use doit:

D.doit()
../_images/3Sympy_95_0.png

You can also use doit on any expression in which D appears.

Integration#

CASs are very well suited to performing differentiation for us, as there’s nothing particularly “clever” going on - the rules for differentiation are relatively straightforward and you could apply them by hand. Integration, on the other hand, can be very difficult. SymPy can integrate many functions that would be annoying to do manually. For example, the function \(f(x) = x^{5} e^{2x} \sin(\frac{x}{3})\) can be integrated using repeated integration by parts, but it would be a little unpleasant. SymPy can compute it easily:

f = x**3 * sym.exp(2*x) * sym.sin(x/3)
f
../_images/3Sympy_98_0.png 
sym.integrate(f, x)
../_images/3Sympy_99_0.png

Definite integrals#

To compute a definite integral, just provide integrate with the limits as well. For example, to compute \(\int_{-3}^4 f(x) dx\) we would do

# add the argument (x, -3, 4)
sym.integrate(f, x, (x, -3, 4))
../_images/3Sympy_101_0.png

To evaluate definite integrals with limits at infinity, use sym.oo (two o’s). For example, to compute \(\int_{-\infty}^0 f(x)dx\):

sym.integrate(f, x, (x, -sym.oo, 0))
../_images/3Sympy_103_0.png

Multiple integrals#

Computing multiple integrals works just like computing multiple derivatives: just give the variables and limits in the order you want the integrals to be performed. As an example, here is how to compute \(\int_{-\infty}^{0} \int_{-2\pi}^{y} x^2 e^y dx\, dy\):

g = x ** 2 * sym.exp(y)

sym.integrate(g, (x, -2*sym.pi, y), (y, -sym.oo, 0))
../_images/3Sympy_105_0.png

Unevaluated integrals#

As with derivatives, you can leave an integral in an unevaluated form. This can be useful to check that you have entered the limits correctly and are integrating in the correct order.

sym.Integral(g, (x, -2*sym.pi, y), (y, -sym.oo, 0))
../_images/3Sympy_107_0.png

To evaluate the integral, use doit:

I = sym.Integral(g, (x, -2*sym.pi, y), (y, -sym.oo, 0))

I.doit()
../_images/3Sympy_109_0.png

Special functions#

Many integrals do not have nice solutions, i.e. cannot be expressed in closed-form in terms of standard functions. Some of these integrals have solutions in terms of special functions like the error function or the gamma function:

sym.integrate(sym.exp(-x**2), x)
../_images/3Sympy_111_0.png 
sym.integrate(sym.exp(-x**3), (x, 0, 3))
../_images/3Sympy_112_0.png

Sometimes SymPy simply cannot perform the integration, in which case it just returns the integral:

sym.integrate(sym.sin(x) * sym.exp(-x**2))
../_images/3Sympy_114_0.png

Plotting#

If we want to produce plots from expressions, it is convenient to convert the expression into a form that numpy can understand. To do this, use the lambdify function:

import numpy as np

T = sym.cos(x) * (x**2 + y**2)

# First argument is the arguments of the new function (in order)
# Second argument is the expression to convert
# Third argument is the module to use (typically numpy)
U = sym.lambdify([x, y], T, "numpy")

Now the function can be plotted as usual:

import matplotlib.pyplot as plt

%matplotlib inline

xp = np.arange(-5, 5, 0.01)
yp = np.arange(-5, 5, 0.01)

X, Y = np.meshgrid(xp, yp)
Z = U(X, Y)

fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, figsize=(8, 8))

ax.plot_surface(X, Y, Z, cmap="coolwarm");
ax.set_xlabel('x',fontsize=14)
ax.set_ylabel('y',fontsize=14)
ax.set_zlabel('T(x,y)',fontsize=14);
../_images/3Sympy_119_0.png

We might want to just look at a “slice” of that surface, say through the plane \(y = 1\). We would do the same thing, but substitute in \(y = 1\) first:

T2 = T.subs(y, 1)
U2 = sym.lambdify(x, T2, "numpy")

%matplotlib inline

X = np.arange(-5, 5, 0.01)
Z = U2(X)

fig, ax = plt.subplots(figsize=(8, 8))

ax.plot(xp, Z);
ax.set_xlabel('x',fontsize=14)
ax.set_ylabel('T(x, 1)',fontsize=14);
../_images/3Sympy_122_0.png