SymPy

November 23, 2023

1 SymPy
SymPy is an open-source Python library for symbolic mathematics. It aims to become
a full-featured computer algebra system (CAS) while keeping the code as simple as
possible in order to be comprehensible and easily extensible. SymPy allows users
to perform various symbolic computations, including algebraic manipulations, solving
equations, calculus, linear algebra, and more.
Key features of SymPy include:
• Symbolic Variables: SymPy allows you to work with symbolic variables, which means
you can manipulate mathematical expressions using variables that are not assigned specific
numerical values.
• Algebraic Manipulations: You can simplify, expand, factorize, and manipulate algebraic
expressions symbolically.
• Equation Solving: SymPy provides tools for solving equations and systems of equations
symbolically. It can find solutions for a wide range of algebraic and transcendental equations.
• Calculus: SymPy supports symbolic calculus operations, such as finding derivatives, inte-
grals, limits, and series expansions.
• Linear Algebra: The library includes functionalities for symbolic linear algebra operations,
allowing you to work with matrices and perform operations like matrix multiplication, inver-
sion, and determinant computation symbolically.
• Trigonometry: SymPy can handle trigonometric expressions symbolically, simplifying and
manipulating them using trigonometric identities.
• Differential Equations: SymPy has tools for solving ordinary differential equations (ODEs)
symbolically.
• Printing and Display: SymPy can generate LaTeX code for mathematical expressions,
making it easy to integrate with documents and presentations.
SymPy is designed to be used both as a standalone library and in conjunction with other scientific
computing libraries, such as NumPy and SciPy. It is a valuable tool for researchers, students,
and professionals in fields such as mathematics, physics, engineering, and computer science where
symbolic mathematics is frequently used. The library is actively developed and has a growing
community of users and contributors.

1
[1]: import sympy as sp

# Define symbolic variables

x, y, z = [Link]('x y z')

[2]: # Define expressions

expr1 = x**2 + 2*x + 1
expr2 = x**2 - 1

# Simplify expressions
simplified_expr1 = [Link](expr1)
expanded_expr2 = [Link](expr2)

print("Simplified expression 1:", simplified_expr1)

print("Expanded expression 2:", expanded_expr2)

Simplified expression 1: x**2 + 2*x + 1

Expanded expression 2: x**2 - 1

[3]: # Define an equation

equation = [Link](x**3 - 4*x**2-x-1, 0)

# Solve the equation

solutions = [Link](equation, x)

print("Solutions:", solutions)

Solutions: [4/3 + (-1/2 - sqrt(3)*I/2)*(sqrt(1005)/18 + 191/54)**(1/3) +

19/(9*(-1/2 - sqrt(3)*I/2)*(sqrt(1005)/18 + 191/54)**(1/3)), 4/3 + 19/(9*(-1/2 +
sqrt(3)*I/2)*(sqrt(1005)/18 + 191/54)**(1/3)) + (-1/2 +
sqrt(3)*I/2)*(sqrt(1005)/18 + 191/54)**(1/3), 19/(9*(sqrt(1005)/18 +
191/54)**(1/3)) + 4/3 + (sqrt(1005)/18 + 191/54)**(1/3)]

1.1 Solve a System of Equations Algebraically

Use SymPy to algebraically solve a system of equations, whether linear or nonlinear. For example,
solving
√ 𝑥2 + 𝑦 = 2𝑧, 𝑦 = −4𝑧 for x and y (assuming z is a constant or parameter) yields {(𝑥 =
6𝑧, 𝑦 = −4𝑧)}
Whether your equations are linear or nonlinear, you can use solve():

[4]: from sympy import solve

from [Link] import x, y, z
solutions = solve([x**2 + y - 2*z, y + 4*z], [x, y], dict=True)

The dict=True parameter indicates that the result should be returned as a list of
dictionaries.
You can then extract solutions by indexing (specifying in brackets) the solution number, and then
the symbol. For example solutions[0][x] gives the result for x in the first solution:

2
[5]: solutions

[5]: [{x: -sqrt(6)*sqrt(z), y: -4*z}, {x: sqrt(6)*sqrt(z), y: -4*z}]

[6]: solutions[0][x]
[6]: √ √
− 6 𝑧

[7]: solutions[1][y]
[7]:
−4𝑧
Here’s an example of finding the derivative and integral of a function:

[8]: # Define a function

f = x**2 + 3*x + 2

# Find the derivative

derivative_f = [Link](f, x)

# Find the integral

integral_f = [Link](f, x)

print("Derivative of f:", derivative_f)

print("Integral of f:", integral_f)

Derivative of f: 2*x + 3
Integral of f: x**3/3 + 3*x**2/2 + 2*x

[9]: [Link](x/(x**2+2*x+1), x)
[9]: 1
log (𝑥 + 1) +
𝑥+1
[10]: import sympy as sp

# Define the symbolic variable

x = [Link]('x')

# Define the integrand

integrand = [Link](x) / x

# Evaluate the definite integral from 0 to infinity

result = [Link](integrand, (x, 0, [Link]))

print("Result of the integral:", result)

Result of the integral: pi/2

SymPy can handle elliptic integrals symbolically. Elliptic integrals are special functions that arise
in various areas of mathematics and physics. Here’s an example of evaluating the complete elliptic
integral of the first kind, denoted as K(m)K(m):

3
[11]: import sympy as sp

# Define the symbolic variable

m = [Link]('m')

# Define the complete elliptic integral of the first kind

elliptic_integral = sp.elliptic_k(m)

print("Complete elliptic integral of the first kind:", elliptic_integral)

Complete elliptic integral of the first kind: elliptic_k(m)

In this example, sp.elliptic_k(m) represents the complete elliptic integral of the first kind with the
parameter mm. SymPy can provide symbolic representations for these elliptic integrals.

[12]: from sympy import symbols, reduce_inequalities, pi

x = symbols('x')
reduce_inequalities([x >= 0, x**2 <= pi], x)
[12]: √
0≤𝑥∧𝑥≤ 𝜋

[13]: reduce_inequalities([x < 0, x > pi], x)

[13]:
False

[14]: from sympy import reduce_inequalities, cos

from [Link] import x, y
from [Link] import periodicity
reduce_inequalities([2*cos(x) < 1, x > 0], x)
[14]: 𝜋 5𝜋
0<𝑥∧ <𝑥∧𝑥<∞∧𝑥<
3 3
[15]: periodicity(2*cos(x), x)
[15]:
2𝜋

1.2 Number Theory

[16]: from sympy import sieve, prime
# All primes less than 19:

print([i for i in [Link](19)])

[2, 3, 5, 7, 11, 13, 17]

[2, 3, 5, 7, 11, 13, 17]

[16]: [2, 3, 5, 7, 11, 13, 17]

[17]: from sympy import prime

prime(10)

4
[17]: 29

[18]: prime(1)

[18]: 2

[19]: prime(100000)

[19]: 1299709

[20]: from sympy import randprime, isprime

randprime(1, 30)

[20]: 19

[Link](n, nth=True)
Returns the product of the first n primes (default) or the primes less than or equal to n (when
nth=False).

[21]: from [Link] import primorial, primerange

from sympy import factorint, Mul, primefactors, sqrt
primorial(4) # the first 4 primes are 2, 3, 5, 7

[21]: 210

[22]: primorial(4, nth=False) # primes <= 4 are 2 and 3

[22]: 6

[23]: primorial(1)

[23]: 2

[24]: from [Link] import primefactors

primefactors(23456)

[24]: [2, 733]

Use SymPy to solve an ordinary differential equation (ODE) algebraically. For example, solving
• 𝑦″ (𝑥) + 9𝑦(𝑥) = 0 yields 𝑦(𝑥) = 𝐶1 𝑆𝑖𝑛(3𝑥) + 𝐶2 𝐶𝑜𝑠(3𝑥)

[25]: from sympy import Function, dsolve, Derivative, checkodesol

from [Link] import x
y = Function('y')
# Solve the ODE

result = dsolve(Derivative(y(x), x, x) + 9*y(x), y(x))

5
 result
[25]:
𝑦(𝑥) = 𝐶1 sin (3𝑥) + 𝐶2 cos (3𝑥)
Check that the solution is correct

[26]: checkodesol(Derivative(y(x), x, x) + 9*y(x), result)

[26]: (True, 0)

Let’s consider an example of solving a non-homogeneous ordinary differential equa-

tion (ODE) using SymPy. Suppose you want to solve the following first-order non-
homogeneous ODE:
𝑑2 𝑦
𝑑𝑥2 + 2𝑦 = 𝑐𝑜𝑠(𝑥)
Here’s the SymPy code to solve this non-homogeneous ODE:

[27]: from sympy import Function, dsolve, Eq, Derivative, symbols, cos

# Define the variable and the function

x = symbols('x')
y = Function('y')

# Define the non-homogeneous ODE

ode = Eq(Derivative(y(x), x, x) + 2 * y(x), cos(x))

# Solve the non-homogeneous ODE

solution = dsolve(ode)

# Print the solution

print("Solution:")
print(solution)

Solution:
Eq(y(x), C1*sin(sqrt(2)*x) + C2*cos(sqrt(2)*x) + cos(x))
The output of checkodesol() is a tuple where the first item, a boolean, tells whether substituting
the solution into the ODE results in 0, indicating the solution is correct.
To solve a non-homogeneous PDE with SymPy, you may need to specify additional
information or make simplifications to obtain a solution. Here’s a modified approach
using pdsolve (partial differential equation solver for PDEs) that might be helpful:
Let’s consider an example of a non-homogeneous first-order partial differential equation (PDE)
using SymPy. Suppose you want to solve the following non-homogeneous PDE:
𝜕𝑢
𝜕𝑡 + 𝑐 𝜕𝑢 𝑡
𝜕𝑥 = 𝑒 𝑥

[28]: from sympy import Function, Eq, Derivative, symbols, pdsolve, exp

# Define the variables

6
 t, x, c = symbols('t x c')

# Define the dependent variable u as a function of t and x

u = Function('u')(t, x)

# Define the non-homogeneous term f(t, x)

f = exp(t) * x

# Define the non-homogeneous PDE

pde = Eq(Derivative(u, t) + c * Derivative(u, x), f)

# Solve the non-homogeneous PDE

solution = pdsolve(pde)

# Print the solution

print("Solution:")
print(solution)

Solution:
Eq(u(t, x), -c*exp(t) + x*exp(t) + F(c*t - x))

1.3 Latex in SymPy

SymPy has built-in support for generating LaTeX code for mathematical expressions. This feature
is helpful if you want to include mathematical expressions in LaTeX documents or presentations.
Here’s an example of how to use LaTeX output in SymPy:

[29]: import sympy as sp

# Define symbolic variables

x, y, z = [Link]('x y z')

# Define a mathematical expression

expr = x**2 + 2*x + 1

# Generate LaTeX code for the expression

latex_code = [Link](expr)

# Print the LaTeX code

print("LaTeX code:", latex_code)

LaTeX code: x^{2} + 2 x + 1

[30]: # Define a more complex expression

expr = [Link](x) + [Link](y) + [Link](z**2, z)

# Generate LaTeX code for the complex expression

latex_code = [Link](expr)

7
 # Print the LaTeX code for the expression
print("LaTeX code:", latex_code)

LaTeX code: \sqrt{y} + \frac{z^{3}}{3} + \sin{\left(x \right)}

[31]: sp.init_printing() # Enables pretty printing with LaTeX in Jupyter

# Display the expression in the notebook

display(expr)

√ 𝑧3
𝑦+ + sin (𝑥)
3

1.4 Limits
[Link](e, z, z0, dir=‘+’) Computes the limit of e(z) at the point z0.

[32]: from sympy import limit, sin, oo

import sympy as sp
# from [Link] import x
x = [Link]('x')

limit(sin(x)/x, x, 0)
[32]:
1

[33]: limit(1/x, x, 0) # default dir='+'

[33]:
∞

[34]: limit(1/x, x, 0,dir='-')

[34]:
−∞

[35]: limit(1/x, x, oo)

[35]:
0

[36]: from sympy import Symbol, cos, series

x = Symbol('x')
series(cos(x),x,0,20)
[36]: 𝑥2 𝑥4 𝑥6 𝑥8 𝑥10 𝑥12 𝑥14 𝑥16
1 − + − + − + − + −
2 24 720 40320 3628800 479001600 87178291200 20922789888000
𝑥18
+ 𝑂 (𝑥20 )
6402373705728000
[37]: from sympy import tan
f = tan(x)
series(f, x, 2, 6, "+")

8
[37]: 2 3
tan (2) + (1 + tan2 (2)) (𝑥 − 2) + (𝑥 − 2) (tan3 (2) + tan (2)) + (𝑥 − 2) ⋅
2 3
1 4 tan (2) 4 5 tan (2) 2 tan (2) 5
( + + tan4 (2)) + (𝑥 − 2) (tan5 (2) + + ) + (𝑥 − 2) ⋅
3 3 3 3
2 17 tan2 (2) 6
( + + 2 tan4 (2) + tan6 (2)) + 𝑂 ((𝑥 − 2) ; 𝑥 → 2)
15 15
• Explanation
• The series function is used to calculate the Maclaurin series expansion.
• f: The function for which the series is calculated (tan(x) in this case).
• x: The variable with respect to which the series is expanded.
• 2: The point around which the expansion is done.
• 6: The number of terms in the series to compute.
• “+”: The direction of the expansion (positive, meaning from the right).

[38]: from sympy import fourier_series, pi

from [Link] import x
s = fourier_series(x**2, (x, -pi, pi))
[Link](2).truncate()
[38]: 2𝜋2
−8 cos (𝑥) + 2 cos (2𝑥) +
3
scale(s)
Scale the function by a term independent of x.

[39]: s
[39]: 𝜋2
−4 cos (𝑥) + cos (2𝑥) + +…
3

1.5 Permutation Groups

PermutationGroup([p1, p2, …, pn]) returns the permutation group generated by the list of permu-
tations. This group can be supplied to Polyhedron if one desires to decorate the elements to which
the indices of the permutation refer.

[40]: from [Link] import Permutation, PermutationGroup

from [Link] import Polyhedron

The permutations corresponding to motion of the front, right and bottom face of a Rubik’s cube
are defined:

[41]: F = Permutation(2, 19, 21, 8)(3, 17, 20, 10)(4, 6, 7, 5)

R = Permutation(1, 5, 21, 14)(3, 7, 23, 12)(8, 10, 11, 9)
D = Permutation(6, 18, 14, 10)(7, 19, 15, 11)(20, 22, 23, 21)

[42]: G = PermutationGroup(F, R, D)
[Link]()
[42]:
3674160

9
[43]: P1 = Permutation(1,2)
P2 = Permutation(3,4)
G = PermutationGroup(P1,P2)
[Link]()
[43]:
4

[44]: G._elements
[44]:
[(4) , (1 2) (4) , (3 4) , (1 2) (3 4)]

[45]: from [Link].named_groups import (SymmetricGroup,

CyclicGroup)
S = SymmetricGroup(6)
C = CyclicGroup(6)
H = [Link](C)
H.is_subgroup(C)

[45]: True

[46]: from [Link].named_groups import (SymmetricGroup,

AlternatingGroup)
S = SymmetricGroup(5)
A = AlternatingGroup(5)
G = [Link](S, A)
G.is_subgroup(A)

[46]: True

[47]: from [Link] import Permutation, SymmetricGroup

S3 = SymmetricGroup(3)
S3.conjugacy_class(Permutation(0, 1, 2))
[47]:
{(0 1 2) , (0 2 1)}

[48]: from [Link].named_groups import SymmetricGroup

from [Link].named_groups import CyclicGroup
S = SymmetricGroup(12)
G = S.sylow_subgroup(2)
C = G.composition_series()
[[Link]() for H in C]
G = S.sylow_subgroup(3)
C = G.composition_series()
[[Link]() for H in C]
G = CyclicGroup(12)
C = G.composition_series()
[[Link]() for H in C]
[48]:
[12, 6, 3, 1]

10
[ ]:

11