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

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)

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

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)

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

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

V.subs([(x, 2), (y, 3), (z, 2)])

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)])

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [12], in <cell line: 3>()
      1 u = sym.symbols("u")
----> 3 V.subs([x, sym.sin(u)])

File ~/opt/anaconda3/envs/jupyterbook/lib/python3.10/site-packages/sympy/core/basic.py:933, in Basic.subs(self, *args, **kwargs)
    931 sequence = list(sequence)
    932 for i, s in enumerate(sequence):
--> 933     if isinstance(s[0], str):
    934         # when old is a string we prefer Symbol
    935         s = Symbol(s[0]), s[1]
    936     try:

TypeError: 'Symbol' object is not subscriptable

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})

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

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

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

u = 8 * x / 3
u

The same is true of roots: instead of doing

x ** (1/3)

You should do

sym.root(x, 3)

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

Q = P.expand()
Q

To go the other way, use factor:

Q.factor()

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)

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)

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

sym.solveset(Q + 7, domain = sym.Reals)

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

sym.solveset(x**2 + 2*x + 1)

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

sym.roots(x**2 + 2*x + 1)

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)

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]

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

\[ \frac{du}{dt} = (u+1)v \ \frac{dv}{dt} = v(u + 2v). \]

That is, we want to solve the system

\[ (u+1)v = 0 \ v(u + 2v) = 0. \]

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

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

sols[1][1]

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)

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

# d^2/dx^2 (V):
sym.diff(V, x, 2)

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

sym.diff(V, x, x)

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

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

sym.diff(V, x)

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)

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

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()

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

sym.integrate(f, x)

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))

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))

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))

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))

To evaluate the integral, use doit:

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

I.doit()

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)

sym.integrate(sym.exp(-x**3), (x, 0, 3))

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

sym.integrate(sym.sin(x) * sym.exp(-x**2))

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);

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);

MT3510 Introduction to Mathematical Computing

SymPy

Contents