## Chapter 05: Symbolic mathematics with Sympy¶

This lab sheet introduces a specific mathematical library called Sympy which allows us to carry out symbolic mathematics.

In the previous week we saw a variety of different libraries:

In [1]:
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.

### 1. Exact calculations¶

A video describing the concept.

A video demo.

Using Python we can calculate the square root and trigonometric values of numbers (we do this by importing the math library)::

In [2]:
import math
math.sqrt(20)

Out[2]:
4.47213595499958
In [3]:
math.cos(math.pi / 4)

Out[3]:
0.7071067811865476

These are fine for numerical work but not when it comes to carrying out exact mathematical calculations, where for example we know that:

$$\cos(\pi / 4) = \sqrt{2} / 2$$

This is where Sympy is useful: it can carry out exact mathematical calculations:

In [4]:
import sympy as sym
sym.sqrt(20)

Out[4]:
2*sqrt(5)
In [5]:
sym.cos(sym.pi / 4)

Out[5]:
sqrt(2)/2

In [6]:
sym.I ** 2

Out[6]:
-1
In [7]:
sym.sqrt(-20)

Out[7]:
2*sqrt(5)*I

Sympy also has numerous functions to manipulate natural numbers:

In [8]:
N = 45 * 63
sym.isprime(N)

Out[8]:
False
In [9]:
sym.primefactors(N)

Out[9]:
[3, 5, 7]
In [10]:
all(sym.isprime(p) for p in sym.primefactors(N))  # All prime factors are prime

Out[10]:
True
In [11]:
sym.factorint(N)

Out[11]:
{3: 4, 5: 1, 7: 1}
In [12]:
N == 3 ** 4 * 5 * 7  # Checking the output of factorint

Out[12]:
True

Repeat the above example with a different value of N.

### 2. Symbolic expressions¶

A video describing the concept.

A Video demo.

Using Sympy it is possible to carry out symbolic computations. To do this we need to creates instances of the Sympy Symbols class.

In [13]:
x = sym.Symbol("x")
x

Out[13]:
x
In [14]:
type(x)

Out[14]:
sympy.core.symbol.Symbol

We can then manipulate this abstract symbolic object without giving it a specific numerical value:

In [15]:
x + x

Out[15]:
2*x
In [16]:
x - x

Out[16]:
0
In [17]:
x ** 2

Out[17]:
x**2

Sympy has a helpful symbols (with a small s) function that lets us create multiple sympy.Symbol objects at a time:

In [18]:
y, z = sym.symbols("y, z")
y, z

Out[18]:
(y, z)

Symbolic expressions can be manipulated using Sympy's:

• factor
• expand

Here we confirm some well known formula:

In [19]:
expr = x ** 2 + 2 * x * y + y ** 2
expr

Out[19]:
x**2 + 2*x*y + y**2
In [20]:
expr.factor()

Out[20]:
(x + y)**2
In [21]:
expr = (x - y) * (x + y)
expr

Out[21]:
(x - y)*(x + y)
In [22]:
sym.expand(expr)  # Note we could also use expr.expand

Out[22]:
x**2 - y**2

Sympy also has a simplify command that can be powerful. Experiment with all of these as well as more complex expressions.

### 3. Symbolic equations¶

A video describing the concept.

A video demo.

We can use Sympy to solve symbolic equations. Let us solve the following symbolic equation:

$$x ^ 2 + 3 x - 2 = 0$$

We do this using the solveset function:

In [23]:
sym.solveset(x ** 2 + 3 * x - 2, x)

Out[23]:
{-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:

$$x^2 + 3x - 2=y$$

1. Modify the equation so that it corresponds to an equation with zero right hand side:

In [24]:
sym.solveset(x ** 2 + 3 * x - 2 - y, x)

Out[24]:
{-sqrt(4*y + 17)/2 - 3/2, sqrt(4*y + 17)/2 - 3/2}

2. Create an Eq object:

In [25]:
eqn = sym.Eq(x ** 2 + 3 * x - 2, y)
sym.solveset(eqn, x)

Out[25]:
{-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):

$$x^2 = -9$$

In [26]:
sym.solveset(x ** 2 + 9, x)

Out[26]:
{-3*I, 3*I}

However if we restrict ourselves to the Reals this is no longer the case:

In [27]:
sym.solveset(x ** 2 + 9, x, domain=sym.S.Reals)

Out[27]:
EmptySet()

### 4. Symbolic calculus¶

A video describing the concept.

A video demo.

It is possible to carry out various symbolic calculus related operations using Sympy:

Let us consider the function:

$$f(x) = 1 / x$$

We will do this by defining a standard Python function:

In [28]:
def f(x):
return 1 / x


and passing it our symbolic variable:

In [29]:
f(x)

Out[29]:
1/x

We can compute the limits at $\pm\infty$

In [30]:
sym.limit(f(x), x, sym.oo)

Out[30]:
0
In [31]:
sym.limit(f(x), x, -sym.oo)

Out[31]:
0
In [32]:
sym.limit(f(x), x, +sym.oo)

Out[32]:
0

We can also compute the limit at $0$ however we must be careful here (you will recall from basic calculus that the limit depends on the direction):

In [33]:
sym.limit(f(x), x, 0)  # The default direction of approach is positive.

Out[33]:
oo
In [34]:
sym.limit(f(x), x, 0, dir="+")

Out[34]:
oo
In [35]:
sym.limit(f(x), x, 0, dir="-")

Out[35]:
-oo

We can use Sympy to carry out differentiation:

In [36]:
sym.diff(f(x), x)

Out[36]:
-1/x**2
In [37]:
sym.diff(sym.cos(x), x)

Out[37]:
-sin(x)

We can carry out various orders of differentiation. These all give the second order derivative of $\cos(x)$:

In [38]:
sym.diff(sym.diff(sym.cos(x), x), x)

Out[38]:
-cos(x)
In [39]:
sym.diff(sym.cos(x), x, x)

Out[39]:
-cos(x)
In [40]:
sym.diff(sym.cos(x), x, 2)

Out[40]:
-cos(x)

As well as differentiation it is possible to carry out integration.

We can do both definite and indefinite integrals:

In [41]:
sym.integrate(f(x), x)  # An indefinite integral

Out[41]:
log(x)
In [42]:
sym.integrate(f(x), (x, sym.exp(1), sym.exp(5)))  # A definite integral

Out[42]:
4

### 5. Differential equations¶

A video describing the concept.

A video demo.

We can use SymPy to solve differential equations. For example:

$$\frac{dy}{dx} = y$$

In [43]:
y = sym.Function('y')
x = sym.symbols('x')
sol = sym.dsolve(sym.Derivative(y(x), x) - y(x), y(x))
sol

Out[43]:
Eq(y(x), C1*exp(x))

Let us verify that the solution is correct:

In [44]:
sym.diff(sol.rhs, x) == sol.rhs

Out[44]:
True

We can also solve higher order differential equations. For example, the following can be used to model the position of a mass on a spring:

$$m\frac{d^2x}{dt^2} + c\frac{dx}{dt} + kx = 0$$

In [45]:
m, c, k, t = sym.symbols('m, c, k, t')
x = sym.Function("x")
sym.dsolve(m * sym.Derivative(x(t), t, 2) + c * sym.Derivative(x(t), t) + k * x(t), x(t))

Out[45]:
Eq(x(t), C1*exp(t*(-c - sqrt(c**2 - 4*k*m))/(2*m)) + C2*exp(t*(-c + sqrt(c**2 - 4*k*m))/(2*m)))

We can solve systems of differential equations like the following:

\begin{aligned} \frac{dx}{dt} & = 1-y\\ \frac{dy}{dt} & = 1-x\\ \end{aligned}

In [46]:
eq1 = sym.Derivative(x(t), t) - 1 + y(t)
eq2 = sym.Derivative(y(t), t) - 1 + x(t)
sym.dsolve((eq1, eq2))

Out[46]:
[Eq(x(t), -C1*exp(-t) - C2*exp(t) + 1), Eq(y(t), -C1*exp(-t) + C2*exp(t) + 1)]

The solution is given as:

\begin{align} x(t) & =-C_1e^{-t}-C_2e^{t} + 1\\ y(t) & =-C_1e^{-t}-C_2e^{t} + 1\\ \end{align}

### 6. Displaying output using $\LaTeX$¶

A video describing the concept.

A video demo.

We can make use of $\LaTeX$ to display the output of Sympy in a human friendly way:

In [47]:
sym.init_printing()
m, c, k, t = sym.symbols('m, c, k, t')
x = sym.Function("x")
sym.dsolve(m * sym.Derivative(x(t), t, 2) + c * sym.Derivative(x(t), t) + k * x(t), x(t))

Out[47]:
$$x{\left (t \right )} = C_{1} e^{\frac{t}{2 m} \left(- c - \sqrt{c^{2} - 4 k m}\right)} + C_{2} e^{\frac{t}{2 m} \left(- c + \sqrt{c^{2} - 4 k m}\right)}$$

On some occasions it might be helpful to be able to turn this off (for example if we wanted to copy and paste the output):

In [48]:
sym.init_printing(False)
m, c, k, t = sym.symbols('m, c, k, t')
x = sym.Function("x")
sym.dsolve(m * sym.Derivative(x(t), t, 2) + c * sym.Derivative(x(t), t) + k * x(t), x(t))

Out[48]:
Eq(x(t), C1*exp(t*(-c - sqrt(c**2 - 4*k*m))/(2*m)) + C2*exp(t*(-c + sqrt(c**2 - 4*k*m))/(2*m)))

## Exercises¶

Here are a number of exercises that are possible to carry out using Sympy:

• Using exact arithmetic;
• Algebraic manipulation of symbolic variables;
• Limits, differentiation and integration;
• Solving differential equations.

### Exercise 1.¶

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

### Exercise 2.¶

Use Sympy's simplify method (and other things) to verify the follow trigonometric identities:

1. $\sin^2(\theta) + \cos^2(\theta) = 1$
2. $2\cos(\theta) \sin(\theta) = \sin(2\theta)$
3. $(1 - \cos(\theta)) / 2 = \sin^2(\theta / 2)$
4. $\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)

### Exercise 3.¶

The point of this question is to investigate the definition of a derivative:

$$\frac{df}{dx}=\lim_{h\to 0}\frac{f(x+h)-f(x)}{h}$$

1. Consider $f(x) = x^3 + 3x - 20$;
2. Compute $\frac{f(x+h)-f(x)}{h}$;
3. Compute the above limit as $h\to 0$ and verify that this is the derivative of $f$.

### Exercise 4.¶

Find the general solutions to the following 4 differential equations:

1. $\frac{dy}{dx}-6y=3e^x$
2. $\frac{dy}{dx}+\frac{x(2x-3)}{x^2+1}=\sin(x)$
3. $\frac{d^2y}{dx^2}-y=\sin(5x)$
4. $\frac{d^2y}{dx^2}+2\frac{dy}{dx}+2x=\cosh(x)$

### Exercise 5.¶

A battle between two armies can be modelled with the following set of differential equations:

$$\begin{cases} \frac{dx}{dt} = - y\\ \frac{dy}{dt} = -5x \end{cases}$$

Obtain the solution to this system of equations.

## Further resources¶

Previous

Next

Source code: @drvinceknight Powered by: Python Jupyter Mathjax Github pages Skeleton css