Lab Sheet 06: Symbolic Mathematics
In the previous week we saw a variety of different libraries:
>>> import random
>>> import math
and there are many more. These libraries are part of the ‘‘standard library’’. This means that they come with the base version of Python. There are a variety of other libraries that exist and are developed independently. Some of these come as standard with Anaconda.
This lab sheet will introduce one such library: SymPy which allows us to do symbolic mathematics.
Building blocks
These questions aim to show you the basic building blocks of programming

Exact valued computations.
A video describing the concept.
Using python we can calculate the square root and trigonometric values of numbers (we do this by importing the
math
library)::>>> import math >>> math.sqrt(20) 4.472... >>> math.cos(math.pi / 4) 0.707...
These are fine for numerical work but not when it comes to carrying out mathematical calculations, where for example we know that:
We can use a python library called SymPy to import various new commands that allow for exact calculations.
>>> import sympy as sym >>> sym.sqrt(20) 2*sqrt(5) >>> sym.cos(sym.pi / 4) sqrt(2)/2
We also have access to complex numbers:
>>> sym.I ** 2 1 >>> sym.sqrt(20) 2*sqrt(5)*I
Whereas this would not have worked in the standard library:
>>> math.sqrt(20) Traceback (most recent call last): ... ValueError: math domain error
Note that if we do want to get the inexact numeric values we can:
>>> sym.sqrt(2).evalf() 1.41...
Experiment with other mathematical functions:
tan
,sin
,exp
. 
Manipulating natural numbers.
A video describing the concept.
As well as exact numerics as described above, SymPy can also be used to explore the factors and primality of a number.
>>> N = 45 * 63 >>> sym.isprime(N) False >>> sym.primefactors(N) [3, 5, 7] >>> all(sym.isprime(p) for p in sym.primefactors(N)) # All prime factors are prime True >>> sym.factorint(N) {3: 4, 5: 1, 7: 1} >>> N == 3 ** 4 * 5 * 7 # Checking the output of `factorint` True
Repeat the above example with a different value of
N
. 
Creating a symbolic variable.
A video describing the concept.
Using SymPy it is possible to carry out symbolic computations. To do this we need to creates instances of the SymPy Symbols class. We do this using the
symbols
function:>>> x = sym.symbols('x') >>> x x >>> type(x) <class 'sympy.core.symbol.Symbol'>
We can also create multiple symbols at a time:
>>> y, z = sym.symbols('y, z') >>> y, z (y, z)
A symbolic variable no longer needs to have a numeric value to be manipulated:
>>> x + x 2*x
Experiment with other symbolic manipulations.

Manipulating symbolic expressions.
A video describing the concept.
We can substitute values in to our symbolic expressions:
>>> expr = x + y + z >>> 2 * expr 2*x + 2*y + 2*z >>> 2 * expr.subs({z:x, y:2}) # Replacing z with x and y with 2 4*x + 4
We can also
factor
,expand
and/orsimplify
them:>>> expr2 = 2 * expr.subs({z:x, y:2}) # Repeating the previous >>> expr2 4*x + 4 >>> expr2.factor() 4*(x + 1)
Here is another example:
>>> alpha, beta = sym.symbols('alpha, beta') >>> ((alpha + beta) ** 2).expand() alpha**2 + 2*alpha*beta + beta**2
Experiment with other expressions and
expand
/simplify
. 
Solving equations.
A video describing the concept.
Sympy has the ability to solve equations. For example let us solve the following quadratic:
We do this using the
solveset
function>>> sym.solveset(x ** 2 + 3 * x  2, x) {3/2 + sqrt(17)/2, sqrt(17)/2  3/2}
If our equation had a non zero right hand side we can use one of two approaches:

Modify the equation so that it corresponds to an equation with zero right hand side:
>>> sym.solveset(x ** 2 + 3 * x  2  y, x) {sqrt(4*y + 17)/2  3/2, sqrt(4*y + 17)/2  3/2}

Create an
Eq
object:>>> eqn = sym.Eq(x ** 2 + 3 * x  2, y) >>> sym.solveset(eqn, x) {sqrt(4*y + 17)/2  3/2, sqrt(4*y + 17)/2  3/2}
We can also specify a domain. For example the following equation has two solutions (it’s a quadratic):
>>> sym.solveset(x ** 2 + 9, x) {3*I, 3*I}
However if we restrict ourselves to the Reals this is no longer the case:
>>> sym.solveset(x ** 2 + 9, x, domain=sym.S.Reals) EmptySet()
Experiment with solving different equations.


Worked example
A video describing the concept.
In this example we are going to investigate the claim that the following polynomial expression can be used to generate primes:
Here are a couple of values:
>>> def polynomial(n): ... return n ** 2  79 * n + 1601 >>> sym.isprime(polynomial(0)) True >>> sym.isprime(polynomial(1)) True >>> sym.isprime(polynomial(2)) True >>> sym.isprime(polynomial(3)) True
Let us first identify for which numbers this is true. We will do this by using a while loop to keep generating the output of this polynomial until it is no longer prime:
>>> n = 0 >>> while sym.isprime(polynomial(n)): # Keep checking values until a composite ... n += 1 >>> n 80
We see that for natural numbers \(n\) from 0 to 80, the quadratic gives primes in our polynomial.
>>> p79 = polynomial(79) >>> sym.isprime(p79), p79 (True, 1601) >>> sym.factorint(p79) {1601: 1} >>> p80 = polynomial(80) >>> sym.isprime(p80), p80 (False, 1681) >>> sym.factorint(p80) {41: 2}
We can list all these primes:
>>> for n in range(80): ... p = polynomial(n) ... print(p, sym.isprime(p)) 1601 True 1523 True ... 1447 True 1523 True 1601 True
Further work
These questions aim to push a bit further.

What is the lowest value of \(n\in\mathbb{Z}\) for which the following polynomials are not prime:

(Optional) Use SymPy to verify that:
Use this and the work done in question 7 to try and explain why \(n^279n+1601\) gives primes for \(0\leq n \leq 79\).

(Optional) Solve the following equations (for \(x\)):
 \(x^2=−1\)
 \(x^2−53x+2a=0\)
 \(xy=\frac{y}{x}\)
Note that not all equations can be solved exactly. There is a whole area of mathematics that studies how to solve equations numerically. One such implementation is in a library called
scipy
. Investigate thenewton
method in scipy (which implements the NewtonRaphson method for finding roots of equations) and use it to solve the following equation (try and solve it using SymPy first):Here is a similar example (finding the solution to \(x^2=9\), which we could of course do exactly):
>>> from scipy.optimize import newton >>> def eqn(x): ... return x ** 2  9 >>> newton(eqn, 0) # How could we obtain the other root? 3.000...

(Optional) Trigonometry
Use SymPy’s
simplify
method (and other things) to verify the follow trigonometric identities: \(\sin^2(\theta) + \cos^2(\theta) = 1\)
 \(2\cos(\theta) \sin(\theta) = \sin(2\theta)\)
 \((1  \cos(\theta)) / 2 = \sin^2(\theta / 2)\)
 \(\cos(n\pi)=(1) ^ n\) (for \(n\in\mathbb{Z}\) (Hint: you will need to
look in to options that can be passed to
symbols
for this).

(Optional) Use SymPy to write the first \(10^6\) prime numbers to file. Compare this file to
primes.csv
(download) (not by hand!) and check that it is the same.