SymPy Demo 0: Intro to SymPy and Jupyter Notebooks

Demo by Karl Johan Måstrup Kristensen and Magnus Troen. Revised 22-10-26 by shsp.

This demo will introduce you to two tools that allow the use of Python as a proper CAS tool:

• SymPy

• Jupyter Notebook

SymPy is a Python package which permits the use of mathematical symbols, equation solving and more. Jupyter Notebooks is a tool for writing Python code and text in the same document for professional presentation.

Jupyter Notebook

Follow DTU’s Python installation guide via https://pythonsupport.dtu.dk/index.html to install all relevant tools including Jupyter Notebook. A Jupyter Notebook file will have the file extension .ipynb.

These SymPy demos are written as Jupyter Notebooks. Download this page via the top icon and follow the rest of the demo in your own editor.

Code vs markdown

To bridge the gap between mathematical output from software and the mathematical reasoning associated with it, Jupyter Notebook allows for executable code and formatted text mixed in one document. A Jupyter Notebook document is built of cells: python cells for code and markdown cells for text.

In text, the markdown syntax allows for headlines with a number of #, bold face with **...**, italic with *...* and more - see a full overview in the official markdown syntax guide. Furthermore, LaTeX code will render when written with $...$ or $$...$$, so that LaTeX expressions like $c=\sqrt{a^2+b^2}$ will render beautifully as: \(c=\sqrt{a^2+b^2}\).

We will not go further into markdown and LaTeX syntax in these demos.

Execution and key combinations

The key combinations presented in these demos work in the code editor Visual Studio Code (VS Code).

• Execute a cell with the key combination Ctrl+Enter, or

• execute and jump to the next cell in one go with Shift+Enter (if there is no next cell, Jupyter Notebook will create one).

Try this out by typing and running the following code in a new python cell in a Jupyter Notebook file:

print("Wow, this output is shown below the cell!")

Wow, this output is shown below the cell!

Note that not just python cells but also markdown cells must be “executed” for markdown to render.

• Press Esc to leave editing mode in a cell without execution.

• Press Enter or double-click to reenter editing mode.

When not in editing mode, press

• m to convert a Python cell to a markdown cell,

• y to convert a markdown cell to a Python cell or

• Del to delete the cell.

Jupyter Notebook automatically prints the value of the last line in a python cell.

2+2
5+4+3+2+1
2+3+4

9

Several comma-separated values will display as a tuple.

2+3+4, 4**2, 5-6

(9, 16, -1)

SymPy

We are now about to load and work with the Python package SymPy whose documentation can be found here: https://docs.sympy.org/latest/reference/index.html. Look up in the documentation by typing ? after the object name or command name (here commented out with # as the package has not been loaded yet):

# diff?

Package import

Your Python installation already comes with several pre-installed Python packages that can be fetched (imported) when needed. Two ways to import packages are:

import cmath
import numpy as np

The contents from these packages are now available by typing the package name (or chosen abbreviation) as a prefix: cmath.polar() or np.array().

An alternative import method that allows us to use the contents without prefixes is:

from sympy import *
init_printing()

from sympy import * imports everything from the SymPy package. If there’s a risk of package conflicts then specify only the content you want: from sympy import sqrt. The command init_printing() forces SymPy to print mathematical output from Python cells using nice LaTeX formatting.

pi**3 + 3 * pi**2 + pi, sqrt(2)

Exact numbers

Python stores values numerically, meaning as a finite number of \(1\)’s and \(0\)’s in memory. This can lead to unexpected results. For example, we can all agree that

(1)\[\begin{equation} 0.3-(0.2+0.1) = 0. \end{equation}\]

But Python disagrees:

0.3-(0.2+0.1)

The output from Python is not exactly \(0\). Consequently, if we ask Python whether they are equal, using the logical operator ==, we get

0 == 0.3-(0.2+0.1)

False

Issues of this kind are avoided by letting SymPy handle the numbers as exact rationals. If our numbers are expressed as fractions, syntax such as S(number), Rational(numerator, denominator) or S('fraction') can be used. SymPy will then keep the numbers on rational form (which is exact in contrast to the approximated numerical form).

Rational(3,10)-(S(2)/10 + S('1/10'))

The S object is called a singleton. See more about singletons in the SymPy documentation.

Mathematical symbols and basic operations

Python does not allow for symbolic mathematics. We can assign values to symbols with y=5, but y=2*x+3 will give an error since x is not defined.

Symbols and arithmetic

This is where SymPy flexes its muscles! With SymPy we can define our variable as a symbol after which it can be used freely.

x = symbols('x')

y=2*x+3
y

When solving equations, reducing expressions and the like, all involved symbols must be defined, which can be done in one go. Note how certain symbols that are called with longer words can be defined as a simpler symbol.

x,y,z,l,p,t = symbols('x y z lambda phi theta')
x,y,z,l,p,t

Many mathematical operations come with the SymPy package, such as the square root sqrt(), the natural exponential function exp(x) and the factorial function factorial(). Note that unusual notation sometimes appears in Python, as is the case for powers \(x^y\), which you obtain with x**y, not with x^y.

sqrt(x), exp(x), factorial(x), x**y, sqrt(16), exp(2), factorial(4), 2**3

Use .evalf() or N(...) to force an exact output down to a numerical, decimal form, where .evalf(no-of-significant-digits) controls how many significant figures to show.

exp(2).evalf(), exp(2).evalf(5), N(exp(2))

Euler’s number is furthermore defined as E, so you can also write:

E**x

Assumptions

By default, defined symbols are treated as complex numbers (that is, without any assumptions). If your symbol represents a number that is real, non-zero or the like, then you can set that as an assumption where the symbol is defined.

# a is a real variable
a = symbols('a', real = True)

# b can't be 0
b = symbols('b', nonzero = True)

# c is a real variable and can't be 0
c = symbols('c', real = True, nonzero = True)

Assumptions let SymPy execute quicker and simplify further, and sometimes they are even essential for getting a useful output, as in the following example.

Example on symbols with and without assumptions

Let us define the expression

\[ \displaystyle\frac{\sqrt{x_1^2 + x_2^2}}{\sqrt{|x_1^2 + x_2^2|}},\quad x_1,x_2 \in \mathbb{R}. \]

x1,x2 = symbols('x1,x2')

expression_x = sqrt(x1**2 + x2**2)/sqrt(abs(x1**2 + x2**2))
expression_x

It seems fine at first sight, but shouldn’t a simplification be possible? You try with .simplify() but see no difference:

expression_x.simplify()

Then you realise that you forgot the assumption of real variables! You adjust the code:

# Symbols with constraints
x1,x2 = symbols('x1,x2', real = True)

expression_y = sqrt(x1**2 + x2**2)/sqrt(abs(x1**2 + x2**2))
expression_y

Wow, a significant simplification. Don’t forget to set assumptions when you define symbols.

Equations

= assigns values to symbols. == performs logical comparison. With the equal sign reserved for this, how, then, can we write equations?

We do that with SymPy’s command Eq(left-hand-side,right-hand-side) in which you comma-separate the left-hand side and right-hand side. For example, the equation \(x^2=\frac13 y\) is written Eq(x**2, S('1/3')*y).

x,y = symbols('x,y')
equation = Eq(x**2, y*S('1/3'))
equation

Note: If you see non-formatted outputs now that we are working with larger expressions, you most likely have forgotten to run init_printing().

Solving equations

For equation solving SymPy offers the command solveset(equation, variable).

sols = solveset(equation, x)
sols

solveset() outputs all solutions in a set. As a set has no ordering, we have to convert the set to a list before we can access specific solutions.

list_of_solutions = list(sols)

list_of_solutions[0], list_of_solutions[1]

Note that Python uses 0-indexing (meaning, the first one in the list has index 0 and is fetched with ...[0], the second one has index 1 etc.)

Alternatively, you can loop through the contents of the set:

for solution in sols:
    display(solution)

Other solver commands exist such as solve(). The benefit of solveset() is that it gives all solutions, even when there are infinitely many. For example, solve the equation \(\sin(x) = 0\).

equation = Eq(sin(x), 0)
solveset(equation, x)

solve() only yields a few of these solutions.

solve(equation, x)

Constraints on solutions

solveset() by default provides all solutions in \(\mathbb{C}\). We can constrain this by providing the wanted domain: solveset(equation, variable, domain). As an example, the equation

x = symbols('x')
equation = Eq(exp(x), 1)
equation

has infinitely many solutions when \(x\in \mathbb{C}\),

solveset(equation, x)

but it only has a single solution when \(x\in\mathbb{R}\),

solveset(equation, x, S.Reals)

Common domains can be stated as follows:

• S.Reals

• S.Complexes

• Interval(0, 1)

Mathematical Functions

Using a Standard Python Function

Python’s native function syntax can be utilised for writing mathematical functions.

# Standard Python function
def f(x):
    return x**2

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

With this definition, we can obtain function values with usual mathematical notation:

f(2), f(10)

While Python functions can handle SymPy objects, they can be problematic along with SymPy functions, such as when differentiating with .diff().

Using Substitutions

An alternative is therefore to define mathematical functions as simple expressions.

x = symbols('x')

# Mathematical function
fx = x**2
fx

We now cannot do fx(2) to find function values, but for that purpose we can use .subs(variable, value) to substitute in a value for \(x\):

# Substituting x = 2
fx.subs(x, 2)

.subs() can also handle dictionaries with the syntax {variable1:value1, variable2:value2, ...} or dict(variable1=value1, variable2=value2, ...) if we wish to substitute in multiple values in one go.

x,y = symbols('x,y')

# Function with two variables
fxy = x*y

# Substituting in x=1, y=2
fxy.subs({x: 1, y: 2}), fxy.subs(dict(x = 1, y = 2))

Mathematical functions defined in this manner can be differentiated with diff(function, variable, order).

diff(fx, x, 1), diff(fx, x, 2) # the first and second derivatives

diff(fxy, y,1), diff(fxy,x,1) # partial derivatives

Complex numbers

In SymPy, the imaginary unit is symbolised I.

# With S(1)/2, 1/2 is displayed as a fraction rather than as 0.5
z = S(1)/2 + I * sqrt(3)/2
z

SymPy provides several tools for work with complex numbers, among others the following - note how the modulus is obtained with abs().

z.as_real_imag(), abs(z), arg(z), conjugate(z)

Remember to try simplify() when an expression does not simplify fully on its own.

z/(5*z)

simplify(z/(5*z))

Switching between complex forms

Converting a complex number from rectangular to polar form is most easily done manually.

z_polar = abs(z)*exp(arg(z)*I)
z_polar

In some cases, but not always, .rewrite(exp) can be used:

t = symbols('t')
w = cos(pi/3 * t) + I*sin(pi/3 * t)
w

w.rewrite(exp)

Rewriting to rectangular form can be done with .rewrite(sin) or .rewrite(cos).

z_rect = z_polar.rewrite(sin)
z_rect

Plots, a brief look

A complete overview of SymPy’s plot functionalities can be found in the SymPy documentation: https://docs.sympy.org/latest/modules/plotting.html.

Note that SymPy draws its plots using the Matplotlib package. SymPy imports it itself, so you don’t have to do that manually. But if Matplotlib is not installed and thus can’t be imported, then your plots will appear as crude text-based sketches instead of figures. Then you must install the Matplotlib package, e.g. by running pip install matplotlib in a terminal.

Let’s define some symbols for the examples below.

x,y,t = symbols('x,y,t')

Plot function of one variable

The plot() command plots functions of one variable. A vast amount of parameters can be set to customise the plot.

func = 3*x**2
plot(func, (x, -5, 10), ylim = [0,100], xlabel = 'x', ylabel = '$3x^2$');

Note how the domain you wish to plot over is set as a parenthesis (variable, from, to) right after the function you are plotting.

Also note that the plot() line is ended with a ;. Without it, a textline with object information would be shown under the plot, which rarely is useful.

The plot() command allows for plotting several expressions in one go.

expression1 = exp(-x/2)*cos(3*x)
expression2 = exp(-x/2)
expression3 = -exp(-x/2)

plot(expression1, expression2, expression3, (x, 0, 10),
     legend=True,
     title='Damped oscillation',
     xlabel='x',
     ylabel='u(x)'
     );

Note how line breaks inside the otherwise rather long line of code give a better overview.

Plot equations and inequalities

The plot_implicit() command plots equations and inequalities. First, an equation:

eq = Eq((x - 1)**2 + (y+1)**2, 10)
plot_implicit(eq, (x,-5, 5), (y,-5, 5), aspect="equal");

And here an inequality:

plot_implicit(x**2 + y**2 <= 4, (x,-3,3), (y,-3,3), aspect="equal");

Plot parametric expressions

The plot_parametric() command plots parametric representations in two dimensions.

parametricrep = Matrix([t*cos(t), sin(t)])
plot_parametric(*parametricrep, (t,0,2*pi));

Plot in 3D

Plotting in 3D can be done with plot3d(), but this command is not included in the SymPy import. It must instead be imported from the sympy.plotting subpackage.

from sympy.plotting import plot3d
plot3d(x**2 - y**2, (x,-2,2), (y,-2,2));

sympy.plotting is a subpackage within the sympy package, so you already have it installed along with SymPy. The from sympy import * line imports SymPy’s most common commands, including some from subpackages (such as plot()), but not everything. The rest must be imported from the relevant subpackage, a separation that has been made to avoid conflicts and clutter. Fortunately, no content in sympy.plotting is in conflict with any top-level content in sympy, so in the future you can safely import everything with * from both at the beginning of your documents.

from sympy import *
from sympy.plotting import *
init_printing()