Chapter 26 Symbolic Maths with SymPy

The basic SymPy commands were covered in lectures; here we will review some additional topics that may be useful to you.

26.1 Further calculus with SymPy

When we discussed calculus we considered differentiation and integration. A precursor to both is the notion of a limit, and SymPy has methods for determining these. Before discussing this further, it is useful to introduce the infinity object in SymPy. For once we will use the from...import syntax, so that we do not have to preface our function with sp. The infinity object is given by oo and so we will enter

from sympy import oo

Sympy correctly carries out arithmetic involving infinity, so 1 + oo will equal oo and so on.

As an example of how to evaluate a limit, suppose that we have defined a symbol x and consider

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

which calculates \[\lim_{x\longrightarrow 0}\left(\frac{\sin(x)}{x}\right)=1\] When we discussed derivatives we saw that the same syntax could be used both for normal and partial derivatives. If we want to take higher derivatives there are two ways to do this, one of which is better suited to multivariable functions. If we enter

d = sp.Derivative(3*x**2+4*x+5, x, 2)

then this will create a derivative object corresponding to the second derivative of the given polynomial (with respect to \(x\)).

Alternatively, we can include a sequence of additional arguments in our function consisting of the variables that we want to differentiate with respect to at each stage. So we could achieve the same result as above by entering

d = sp.Derivative(3*x**2+4*x+5, x, x)

This has the advantage that we can now take higher partial derivatives. For example

d = sp.Derivative(3*x**2*y**2+4*x+3*y+5, x, y)

will calculate the mixed partial derivative of the given polynomial where we first differentiate with respect to \(x\) and then with respect to \(y\).

SymPy can also work with sequences and series. We can define a sequence as in the following example (where we assume that we have previously defined n to be a symbol)

a_n = 1/(n**2)

and then

sp.limit(a_n, n, oo)

calculates \[\lim_{n\longrightarrow \infty}a_n=0\] and

sp.summation(a_n, [n, 1, oo])

calculates \[\sum_{n=1}^\infty a_n=\frac{\pi^2}{6}.\] We can also just work with the formula for the \(n\)th term directly, for example

sp.limit((1+1/n)**n, n, oo)

will evaluate \[\lim_{n\longrightarrow \infty}\left(1+\frac{1}{n}\right)^n\] and return the answer \(e\) (ie the base of the natural logarithms).

26.2 Further algebra with SymPy

A matrix in SymPy is defined in terms of a list of lists. For example

a = sp.Matrix([[1,2,3],[4,5,6], [7,8,9]])

corresponds to the matrix \[\left(\begin{array}{ccc} 1 & 2 & 3\\ 4 & 5 & 6\\ 7 & 8 & 9\end{array}\right).\] To access individual elements we do not need to consider them as we did in standard Python in terms of elements of elements, but instead can just use the usual matrix notation — but starting our count of rows and columns from zero. So a[1,2] refers to the element in row 2 and column 3 of our matrix a, ie the number 6.

Special cases of matrices are row and column vectors, and so SymPy does not have a separate notation for these. Given two row (or two column) vectors u and v of the same dimension we can consider the dot product, u.dot(v), or (for vectors in three dimensions) the cross product u.cross(v).

Many of the operations we learnt in Linear Algebra have analogues in SymPy. For example, we can determine the reduced row echelon form of a matrix a by entering a.rref()[0]. We need to include the [0] here as the rref method actually returns two pieces of information as a tuple, and the desired matrix is in the first part of this.

The determinant of a square matrix a is given by a.det() and the inverse (if it exists) by a.inverse(). If a has no inverse this will result in an error.

The method a.eigenvals() will return the eigenvalues of a and their multiplicities. For example if

a = sp.Matrix([[2,3],[3,2]])

then a.eigenvals() will return

{5: 1, -1: 1}

which gives the two eigenvalues 5 and -1, each with multiplicity 1. Similarly a.eigenvects() will return the eigenvalues together with their corresponding multiplicities and respective bases of eigenvectors.

Finally, if a is a diagonalisable matrix then we can diagonalise it using the command

p, d = a.diagonalize()

This will produce two new matrices p and d such that \[a=pdp^{-1}.\]