Symbolic Maths with SymPy

Algebra, calculus, matrices and exact arithmetic — done symbolically, with the working shown, on the device in front of you.

Most calculators give you a number. sympy gives you the expression: it factorises, differentiates, integrates, solves for x, and hands back pi/2 rather than 1.5707963. It is a full computer algebra system, it is already here, and almost nobody knows it.

This page is a lesson. Type each block, run it, look at what comes back, then change a number and run it again. Every output shown below is the real output from the version that runs on your device.

Where to type it

Two places, both with a Run button and an output panel underneath:

• A CodeBook cell — press Run Cell Run Cell or ⇧⏎. This is the better one for a lesson, because variables stay alive from one Python cell to the next: define x in the first cell and every cell after it can use x. The blocks below are written to be run in order, one per cell.
• A Python code document — one file, one Run button. Each run starts fresh, so paste a whole block at a time.

It runs on the device you are holding — phone, tablet or computer. The one exception is the very first run: the library is fetched then (about seventeen megabytes), so that first run needs a connection and takes a moment. Every run after that on the same device starts straight away.

You never install anything by hand. Writing import sympy is enough — the library is fetched because your code asked for it.

Five minutes in

Start here. Two symbols, and three things you cannot do with a calculator.

from sympy import symbols, expand, factor, simplify

x, y = symbols('x y')

print(expand((x + y)**3))
print(factor(x**2 - 4))
print(simplify((x**2 - 1) / (x - 1)))

x**3 + 3*x**2*y + 3*x*y**2 + y**3
(x - 2)*(x + 2)
x + 1

Three things just happened that are worth naming.

symbols('x y') makes maths, not data. x is not a number waiting for a value — it is the letter x, and it stays a letter. x + x is 2*x. Everything on this page rests on that one idea.

Expressions are things you can hold. (x + y)**3 is an object you can pass around, store in a variable, and hand to expand or factor or diff later. Nothing is evaluated until you ask.

Nothing was approximated. (x**2 - 1)/(x - 1) became exactly x + 1, by cancelling — not by trying numbers and rounding.

Now change something. Make it (x + y)**5 and run it again. Try factor(x**3 - 8). Try simplify on a fraction of your own. Nothing here can break.

Making it readable

print gives you a single line — correct, copyable, and hard on the eyes as soon as anything gets nested:

from sympy import symbols, Integral, sqrt

x = symbols('x')
print(Integral(sqrt(1/x), x))

Integral(sqrt(1/x), x)

Use pprint instead and you get the thing as it would be written by hand:

from sympy import symbols, pprint, Integral, sqrt

x = symbols('x')
pprint(Integral(sqrt(1/x), x))

  /
 |
 |     ___
 |    / 1
 |   /  -  dx
 | \/   x
 |
/

That is the expression drawn in two dimensions, out of plain keyboard characters. Ask for the rounded ones and the same thing gets a good deal easier on the eye:

from sympy import symbols, pprint, Integral, sqrt

x = symbols('x')
pprint(Integral(sqrt(1/x), x), use_unicode=True)

⌠
⎮     ___
⎮    ╱ 1
⎮   ╱  ─  dx
⎮ ╲╱   x
⌡

Either form arrives in the output panel intact: the panel keeps every space and every line break, and draws them in a fixed-width font, so the columns line up the way they were built. Fractions stack properly too:

from sympy import symbols, pprint

x = symbols('x')
pprint((x**2 + 2*x + 1)/(x + 1), use_unicode=True)

 2
x  + 2⋅x + 1
────────────
   x + 1

Three things to know about this, because the usual advice is only half right here.

Every tutorial elsewhere opens with init_printing(), and it is worth calling: init_printing(use_unicode=True) makes the rounded characters the default, so every pprint after it can leave the argument off. In a CodeBook you need it once, in an early cell, and it holds for every cell after. What it does not do is change print — that still gives the one-line form whatever you set.

pprint draws straight to the output. If you want the drawing as text — to put in a variable, or return it — pretty(expr, use_unicode=True) gives you exactly the same picture as an ordinary string.

Wide results wrap. The panel wraps long lines rather than running off the side, and a wrapped picture loses its alignment. Most things are narrow enough not to care; when one isn't, drag the top edge of the output panel to make it taller, or fall back to print for that one line.

Solving equations

solve takes an expression and the letter to solve for, and assumes the expression equals zero.

from sympy import symbols, solve, Eq

x, y = symbols('x y')

print(solve(x**2 - 5*x + 6, x))
print(solve(Eq(x**2, 2), x))
print(solve([Eq(x + y, 10), Eq(2*x - y, 2)], [x, y]))

[2, 3]
[-sqrt(2), sqrt(2)]
{x: 4, y: 6}

Note the second answer: -sqrt(2) and sqrt(2), not -1.414…. It stays exact until you ask for digits.

Use Eq(left, right) when you want to write an equation with both sides, and a plain expression when it is already equal to zero. Give solve a list of each and a list of the unknowns to solve a system.

The part that makes this different from a calculator is that the coefficients can be letters too:

from sympy import symbols, solve

a, b, c, x = symbols('a b c x')
print(solve(a*x**2 + b*x + c, x))

[(-b - sqrt(-4*a*c + b**2))/(2*a), (-b + sqrt(-4*a*c + b**2))/(2*a)]

That is the quadratic formula, derived rather than remembered. Feed it a cubic and it will derive that one too.

When there is no formula. Some equations have no answer expressible in symbols. solve says so honestly rather than guessing — a fifth-degree polynomial comes back as CRootOf(...) placeholders, which are exact references to the real roots. When you want numbers instead, ask for numbers:

from sympy import symbols, solve, nsolve, cos, N

x = symbols('x')
print([N(r, 10) for r in solve(x**5 - x - 1, x)][:2])
print(nsolve(cos(x) - x, x, 1))

[1.167303978, -0.7648844336 - 0.3524715460*I]
0.739085133215161

nsolve is the one to reach for when an equation mixes ordinary functions in a way nothing can untangle — cos(x) = x has a perfectly good answer and no symbolic form. Give it a starting guess.

Calculus

This is where the working being shown matters most.

Differentiating

from sympy import symbols, diff, sin, exp

x = symbols('x')

print(diff(x**3 + 2*x, x))
print(diff(sin(x)*exp(x), x))
print(diff(x**4, x, 2))

3*x**2 + 2
exp(x)*sin(x) + exp(x)*cos(x)
12*x**2

The third argument is how many times: diff(x**4, x, 2) is the second derivative.

Integrating

Same shape. One argument integrates indefinitely; a tuple (x, from, to) gives you a definite integral.

from sympy import symbols, integrate, sin, exp, oo, pi

x = symbols('x')

print(integrate(2*x, x))
print(integrate(sin(x), (x, 0, pi)))
print(integrate(exp(-x**2), (x, -oo, oo)))

x**2
2
sqrt(pi)

oo — two lower-case letter o's — is infinity. That last line is the Gaussian integral, and the answer is exactly the square root of pi.

Limits and series

from sympy import symbols, limit, series, sin, cos, oo

x = symbols('x')

print(limit(sin(x)/x, x, 0))
print(limit((1 + 1/x)**x, x, oo))
print(series(cos(x), x, 0, 8))

1
E
1 - x**2/2 + x**4/24 - x**6/720 + O(x**8)

E is Euler's number, exactly. The O(x**8) on the end of a series is the honest part: it marks where the expansion was cut off.

A worked problem

Throw a ball upward at 15 metres per second from a height of 20 metres. When does it stop rising, how high does it get, and when does it land?

from sympy import symbols, Eq, solve, diff

t = symbols('t')

h = 20 + 15*t - 4.9*t**2
v = diff(h, t)

print('velocity:', v)
print('time at the top:', solve(Eq(v, 0), t))
print('height at the top:', h.subs(t, solve(Eq(v, 0), t)[0]))
print('back at ground level:', solve(Eq(h, 0), t))

velocity: 15 - 9.8*t
time at the top: [1.53061224489796]
height at the top: 31.4795918367347
back at ground level: [-1.00402905068862, 4.06525354048453]

Nothing was plotted and nothing was guessed. The velocity is the derivative of the height; the top is where the velocity is zero; .subs puts a value back into an expression. The negative landing time is real and worth keeping — it is when the ball would have left the ground had it been thrown from there.

Differential equations

dsolve solves for a function rather than a number. Declare it with Function, write the equation with .diff(x), and read the answer with its constants of integration in place.

from sympy import symbols, Function, dsolve, Eq

x = symbols('x')
f = Function('f')

print(dsolve(Eq(f(x).diff(x), f(x)), f(x)))
print(dsolve(Eq(f(x).diff(x, 2) + f(x), 0), f(x)))

Eq(f(x), C1*exp(x))
Eq(f(x), C1*sin(x) + C2*cos(x))

The first is the equation that says a thing grows at a rate equal to its own size, and the answer is exponential growth. The second is the equation of anything that oscillates, and the answer is sine and cosine. C1 and C2 are the unknown constants.

Pin them down by giving the starting conditions:

from sympy import symbols, Function, dsolve, Eq

x = symbols('x')
f = Function('f')

print(dsolve(Eq(f(x).diff(x), -2*f(x)), f(x), ics={f(0): 5}))

Eq(f(x), 5*exp(-2*x))

That is a decay curve starting at 5, derived from the statement that it decays at twice its size.

Matrices and eigenvalues

Matrix takes a list of rows. Arithmetic, determinants and inverses are all exact.

from sympy import Matrix, pprint

M = Matrix([[1, 2], [3, 4]])

pprint(M, use_unicode=True)
print(M.det())
pprint(M.inv(), use_unicode=True)
pprint(M * M, use_unicode=True)

⎡1  2⎤
⎢    ⎥
⎣3  4⎦
-2
⎡-2    1  ⎤
⎢         ⎥
⎣3/2  -1/2⎦
⎡7   10⎤
⎢      ⎥
⎣15  22⎦

Note 3/2 in the inverse — a fraction, not 1.5.

Solving a linear system is one call:

from sympy import Matrix, pprint

A = Matrix([[2, 1, -1], [-3, -1, 2], [-2, 1, 2]])
rhs = Matrix([8, -11, -3])

pprint(A.solve(rhs), use_unicode=True)
print('rank', A.rank(), 'determinant', A.det())

⎡2 ⎤
⎢  ⎥
⎢3 ⎥
⎢  ⎥
⎣-1⎦
rank 3 determinant -1

And eigenvalues come back exact, as a dictionary of value to how many times it occurs:

from sympy import Matrix

M = Matrix([[2, 1], [1, 2]])

print(M.eigenvals())
print(M.eigenvects())

{3: 1, 1: 1}
[(1, 1, [Matrix([
[-1],
[ 1]])]), (3, 1, [Matrix([
[1],
[1]])])]

Here is a real edge, and it will surprise you. Those numbers were clean because the matrix was chosen to be. Eigenvalues of an everyday 3×3 are exact and enormous — nested cube roots of complex numbers, running to several lines, all of it correct and none of it readable. Exactness is not always what you want. When it isn't, ask for numbers:

from sympy import Matrix

M = Matrix([[4, 1, 2], [1, 3, 0], [2, 0, 5]])

print([complex(v).real for v in M.eigenvals()])

[6.669079088282288, 1.8548973087995775, 3.476023602918134]

That is the general lesson of this page in one line: stay exact while it is useful, and convert when it stops being useful.

Exact numbers, and as many digits as you like

Ordinary arithmetic on a computer is approximate, and quietly so:

from sympy import Rational, sqrt, S

print(0.1 + 0.2)
print(Rational(1, 10) + Rational(2, 10))
print(sqrt(8))
print(S(1)/3)

0.30000000000000004
3/10
2*sqrt(2)
1/3

The first line is not a bug in this library — that is what a decimal fraction does in ordinary floating-point arithmetic, everywhere. Rational and S step around it by keeping the numbers as fractions.

When you do want digits, say how many:

from sympy import N, pi, sqrt

print(N(pi, 50))
print(N(sqrt(2), 30))
print(N(pi**pi, 40))

3.1415926535897932384626433832795028841971693993751
1.41421356237309504880168872421
36.46215960720791177099082602269212366637

There is no practical ceiling on that second argument — a couple of thousand digits of pi comes back in about a millisecond. And exact arithmetic on whole numbers has no size limit at all:

from sympy import factorial, factorint, prime

print(factorial(30))
print(factorint(600851475143))
print(prime(1000))

265252859812191058636308480000000
{71: 1, 839: 1, 1471: 1, 6857: 1}
7919

factorint gives the prime factorisation as a dictionary of prime to power; prime(1000) is the thousandth prime.

Typesetting the answer

latex() converts any expression into LaTeX — the notation used to typeset mathematics in books and papers.

from sympy import symbols, latex, Integral, sqrt, Matrix

x = symbols('x')

print(latex(Integral(sqrt(1/x), x)))
print(latex((x**2 + 1)/(x - 1)))
print(latex(Matrix([[1, 2], [3, 4]])))

\int \sqrt{\frac{1}{x}}\, dx
\frac{x^{2} + 1}{x - 1}
\left[\begin{matrix}1 & 2\\3 & 4\end{matrix}\right]

Those three lines are what typeset mathematics looks like as text. Pasted into a text document between $$ marks, they come out like this:

What it is for: you did the algebra here, and you want the result to appear properly set in something you are writing. Copy that output into a text document, wrap it in $$ at each end, and it renders as above rather than as code. Chat replies render it the same way. The LaTeX reference covers the notation and where it renders.

This is the step that closes the loop: derive it with sympy, check it with pprint, publish it with latex(), without ever retyping the expression and risking a slip.

Where the edges are

Said plainly, so you do not have to find them yourself.

Try it

Nothing here needs anything set up. Pick one.

1. Check your own homework. Differentiate something by hand, then print(diff(your_expression, x)) and compare.
2. Derive a formula you were made to memorise. solve(a*x**2 + b*x + c, x) is one. The sum of a geometric series is another.
3. Watch a decay curve fall out of its own definition. dsolve(Eq(f(x).diff(x), -k*f(x)), f(x), ics={f(0): 100}), with k a symbol.
4. Get a hundred digits of something. N(sqrt(2), 100), or N(pi**2/6, 100).
5. Typeset a result. Derive something, run it through latex(), and paste it into a text document between $$ marks.

Related

• Languages, and Where They Run — what Python can do on this device, on a connected computer, and in an unattended run.
• CodeBook Tour — every control in the notebook, including where a cell's code runs.
• Code IDE Tour — the code document, its Run button and its output panel.
• LaTeX Reference — the notation latex() produces, and the places it renders.
• Coming from Jupyter — how notebook habits carry across.