sympy cheat sheet

from sympy import *         # fine for scripts/notebooks
import sympy as sp          # safer in a package

x, y, z = symbols('x y z')  # Symbols must be DECLARED
n, k    = symbols('n k', integer=True)
f, g    = symbols('f g', cls=Function)
t       = Symbol('t', positive=True)

init_printing()             # pretty output in a notebook

• symbols('x y z')
  Space- or comma-separated. Returns a tuple.
• symbols('a0:5')
  Ranges: a0, a1, a2, a3, a4.
• from sympy.abc import x, y
  Shortcut for the obvious single letters.
• Function('f')
  An undefined function — f(x).diff(x) stays symbolic.

Gotchafrom sympy import * shadows built-ins — notably abs, pow, sum, and lambda-adjacent names. Fine in a notebook; use import sympy as sp in library code.

x = symbols('x')
expr = x + 1
x = 2            # rebinds the PYTHON name, not the Symbol
print(expr)      # x + 1     <- not 3

expr.subs(x, 2)  # 3         <- this is how you substitute

GotchaRebinding the Python variable does nothing to expressions already built from the Symbol. The tree holds the Symbol object, not the name. Verified: expr is still x + 1.

expr = x**2 + 2*x*y

expr.func           # Add            - the top node
expr.args           # (x**2, 2*x*y)  - its children
expr.args[0].args   # (x, 2)
srepr(expr)         # the exact constructor calls
expr.free_symbols   # {x, y}
expr.atoms(Number)  # {2}
count_ops(expr)     # 4
list(preorder_traversal(expr))

Tipsrepr() is the debugger for “why won't this simplify?”. It shows you what SymPy actually built, not what you meant.

Rational(1, 3)                      # 1/3   exact
S(1)/3                              # 1/3   S() = sympify
Integer(5) / 2                      # 5/2
S.Half                              # 1/2

x + 1/2                             # x + 0.5   <-- float leaked in!
x + Rational(1, 2)                  # x + 1/2
x + S(1)/2                          # x + 1/2

nsimplify(0.333333, rational=True)  # rescue a float -> 1/3
sympify("x**2 + 1")                 # string -> expression

Gotcha1/2 is evaluated by Python before SymPy ever sees it. By the time it reaches your expression it is already 0.5. Verified: x + 1/2 → x + 0.5.

• oo, -oo, zoo
  Infinity, negative infinity, complex infinity. Not inf.
• pi, E, I
  Exact π, e, and the imaginary unit. I**2 == -1.
• nan
  SymPy's NaN.
• S.Reals, S.Integers, S.Complexes
  Domains, for solveset.
• GoldenRatio, EulerGamma, Catalan
  Named constants, exact.

sqrt(8) auto-simplifies to 2*sqrt(2) — exactly, with no float anywhere. That is the whole point of the library.

sqrt(2).evalf()  # 1.41421356237310
pi.evalf(50)     # 50 digits, on demand
N(sqrt(2), 30)   # same thing
(x + pi).evalf(subs={x: 1})

float(sqrt(2))   # a real Python float, when you finally need one

Tipevalf uses mpmath, so precision is arbitrary. Ask for 1000 digits of π and you will get them.

• <code>a = b</code>
  Python assignment. Nothing symbolic happens.
• <code>a == b</code>
  Structural equality → a Python bool. (x+1)**2 == x**2+2*x+1 is False.
• <code>Eq(a, b)</code>
  A symbolic equation object. Feed it to solve.
• <code>simplify(a - b) == 0</code>
  The reliable test for mathematical equality.
• <code>a.equals(b)</code>
  Probes numerically at random points. Fast, and usually right.

GotchaThere is no fully general algorithm for symbolic equality (Richardson's theorem), which is why == only promises structure. simplify(a-b) is the pragmatic answer.

Eq(x, y);  Ne(x, y);  Lt(x, y);  Le(x, y);  Gt(x, y);  Ge(x, y)

Piecewise((x,  x > 0),
          (-x, True))  # True = the else branch

if x > 0:  ...         # TypeError!

Gotchaif x > 0: raises TypeError: cannot determine truth value of Relational. A symbolic inequality has no truth value until x is known. Verified.

x = symbols('x')                   # generic: could be complex
p = symbols('p', positive=True)
n = symbols('n', integer=True, nonzero=True)

sqrt(x**2)                         # sqrt(x**2)  <- cannot simplify: x might be negative
sqrt(p**2)                         # p           <- now it can

refine(sqrt(x**2), Q.positive(x))  # assume after the fact
ask(Q.positive(x**2 + 1))          # query the system
p.is_positive                      # True
x.is_positive                      # None  = "unknown", not False

• common
  real, positive, negative, nonnegative, integer, rational, nonzero, even, odd, prime, finite
• <code>.is_*</code> is three-valued
  True / False / None. None means unknown.

GotchaSymbol('x') != Symbol('x', positive=True) — they are different symbols. Mix them in one expression and nothing will cancel. Declare once, at the top.

simplify(sin(x)**2 + cos(x)**2)  # 1
simplify(expr, ratio=1.7)        # only accept if it got smaller

Gotchasimplify tries dozens of transformations and keeps whatever scores smallest by count_ops. It is slow, and “smallest” is often not the form you wanted. If you know the operation you need, call it by name.

• expand(expr)
  Multiply everything out.
• factor(expr)
  The inverse. Over the rationals by default.
• collect(expr, x)
  Group by powers of x.
• cancel(expr)
  Rational function → lowest terms, p/q.
• apart(expr, x)
  Partial fractions.
• together(expr)
  Back over one denominator.
• radsimp(expr)
  Rationalize the denominator.
• expand(expr, deep=False)
  Only the top level.

• trigsimp / expand_trig
  Contract with identities / expand angle sums.
• powsimp / expand_power_exp / powdenest
  Combine and split exponents.
• logcombine / expand_log
  Combine and split logs. Both usually need force=True.
• expr.rewrite(exp)
  sin(x) → complex exponentials.
• expr.rewrite(gamma)
  factorial(n) → gamma(n+1).
• nsimplify(0.5)
  Float → the exact thing it probably was.
• hyperexpand(expr)
  Evaluate hypergeometric functions.

Gotchaexpand_log(log(x*y)) returns it unchanged unless x, y are declared positive — because log(xy) = log x + log y is false for negative reals. Either declare the assumption or pass force=True.

expr.subs(x, 2)                                 # returns a NEW expression
expr.subs({x: 2, y: 3})                         # dict
expr.subs([(x, y), (y, 2)])                     # SEQUENTIAL -> 4
expr.subs([(x, y), (y, 2)], simultaneous=True)  # -> y + 2

expr.xreplace({x: 2})                           # exact structural, faster, no sympify
expr.replace(sin, cos)                          # pattern-level swap

GotchaA list of substitutions is applied one after another, so the second can hit the result of the first. Verified on x + y: sequential gives 4, simultaneous=True gives y + 2.

diff(sin(x)*x**2, x)  # 2*x*sin(x) + x**2*cos(x)
diff(expr, x, 2)      # second derivative
diff(expr, x, 2, y)   # mixed partials
expr.diff(x)          # method form

Derivative(f(x), x)   # UNEVALUATED
Derivative(f(x), x).doit()

f = Function('f')
f(x).diff(x)          # stays symbolic: Derivative(f(x), x)

• idiff(eq, y, x)
  Implicit differentiation.
• Matrix([...]).jacobian([x, y])
  Jacobian.
• hessian(expr, [x, y])
  Hessian.
• differentiate_finite(f(x).diff(x))
  Finite-difference approximation. From sympy.calculus.finite_diff.

integrate(x**2, x)                   # x**3/3   (no +C, ever)
integrate(x**2, (x, 0, 1))           # 1/3
integrate(exp(-x**2), (x, -oo, oo))  # sqrt(pi)
integrate(f, (x, 0, 1), (y, 0, 1))   # double

Integral(sin(x), (x, 0, pi))         # unevaluated
Integral(sin(x), (x, 0, pi)).doit()  # 2

GotchaTwo things integrate may hand you instead of an expression: an unevaluated Integral (it gave up), or a Piecewise. Verified: integrate(x**y, x) → Piecewise((x**(y+1)/(y+1), Ne(y,-1)), (log(x), True)). Check the type before you keep computing.

limit(sin(x)/x, x, 0)       # 1
limit(1/x, x, 0, '+')       # oo      <- direction matters
limit(1/x, x, 0, '-')       # -oo
limit((1 + 1/x)**x, x, oo)  # E

sin(x).series(x, 0, 6)      # x - x**3/6 + x**5/120 + O(x**6)
sin(x).series(x, 0, 6).removeO()
expr.series(x, oo, 3)       # asymptotic expansion
fourier_series(x, (x, -pi, pi)).truncate(3)

The O(x**6) term is a real object and it is contagious — it will swallow higher-order terms in later arithmetic. Strip it with .removeO() when you are done.

summation(k, (k, 1, n))         # n*(n + 1)/2   - closed form!
Sum(1/k**2, (k, 1, oo)).doit()  # pi**2/6
product(k, (k, 1, n))           # factorial(n)

Sum(k, (k, 1, n))               # unevaluated, printable

TipSymPy will find the closed form of a symbolic sum when one exists. That is something NumPy fundamentally cannot do.

solve(x**2 - 4, x)                   # [-2, 2]        - expression assumed = 0
solve(Eq(x**2, 4), x)                # same
solve([x + y - 2, x - y], [x, y])    # {x: 1, y: 1}
solve(x**2 - 4, x, dict=True)        # [{x: -2}, {x: 2}]  - stable shape

solveset(x**2 - 4, x)                # {-2, 2}        - a Set
solveset(sin(x), x, domain=S.Reals)  # ALL roots, as ImageSets
solveset(exp(x), x)                  # EmptySet

GotchaVerified: solve(sin(x), x) returns just [0, pi]. There are infinitely many roots. solveset returns all of them. If you are iterating over the result of solve, you may be iterating over a lie.

• linsolve([...], (x, y))
  Linear systems. Also takes an augmented Matrix.
• nonlinsolve([...], [x, y])
  Nonlinear systems → a set of tuples.
• nsolve(cos(x) - x, x, 1)
  Numeric root, needs a starting guess.
• roots(Poly(x**3 - 1, x))
  Roots with multiplicities, as a dict.
• real_roots / nroots
  Only the real ones / all of them numerically.
• solve_univariate_inequality(x**2 > 4, x)
  Inequalities → a set.
• reduce_inequalities([x > 2, x < 5], [x])
  Systems of inequalities.
• diophantine(x**2 + y**2 - 25)
  Integer solutions.
• linear_eq_to_matrix(eqs, [x, y])
  Equations → A, b.

f = Function('f')

dsolve(Eq(f(x).diff(x, 2) + f(x), 0), f(x))
#   Eq(f(x), C1*sin(x) + C2*cos(x))

dsolve(f(x).diff(x) - f(x), f(x), ics={f(0): 1})  # apply initial conditions
classify_ode(f(x).diff(x) - f(x), f(x))           # which methods apply
checkodesol(ode, sol)                             # verify it!

from sympy.solvers.pde import pdsolve

dsolve returns an Eq, with constants C1, C2, .... Pass ics= to pin them down — or solve for them yourself.

A = Matrix([[1, 2], [3, 4]])
eye(3);  zeros(2, 3);  ones(2, 2);  diag(1, 2, 3)

A + B;  A * B  # * IS matrix multiplication (unlike NumPy!)
A**2;  A.T;  A.inv();  A.det();  A.trace()
A.multiply_elementwise(B)

A[0, 1];  A[:, 0];  A.row(0);  A.shape
A.row_join(B);  A.col_join(B)

GotchaIn SymPy A * B is the matrix product. In NumPy it is elementwise. Same operator, opposite meaning — and both accept the other's arrays without complaint.

A.rref()                # (reduced form, pivot columns)
A.nullspace();  A.columnspace();  A.rank()

A.eigenvals()           # {eigenvalue: multiplicity}  - EXACT
A.eigenvects()          # [(val, mult, [vectors])]
P, D = A.diagonalize()  # A == P*D*P.inv()
A.charpoly();  A.jordan_form();  A.singular_values()

A.LUdecomposition();  A.QRdecomposition();  A.cholesky()
A.solve(b);  A.LUsolve(b)

TipSymbolic eigenvalues of a 2×2 with symbol entries come out as exact radicals. Try that in NumPy.

• MatrixSymbol('A', n, n)
  A matrix whose size is symbolic. A*B stays unevaluated.
• ImmutableMatrix
  Hashable — usable inside expressions and as dict keys.
• SparseMatrix(m, n, {(i,j): v})
  Sparse.
• M.jacobian([x, y]) / hessian(f, [x, y])
  Vector calculus.
• M.applyfunc(simplify)
  Map a function over every entry — often what you actually need.

GotchaMatrix is mutable and therefore unhashable. If you need it inside an expression or a set, use ImmutableMatrix.

import numpy as np
expr = sin(x)**2 + sqrt(x + 1)

f = lambdify(x, expr, 'numpy')      # -> a real Python function
f(np.linspace(0, 5, 100_000))       # vectorised, fast

g = lambdify((x, y), x*y, 'numpy')  # several args
h = lambdify(x, expr, cse=True)     # factor out repeated subexpressions
lambdify(x, expr, 'mpmath')         # arbitrary precision instead

GotchaNever loop over .subs() to evaluate an expression at many points. Measured on this machine: subs + float = 1263 µs per point; lambdify + NumPy = 0.22 µs per point. That is a 5,780× difference.

print(expr)          # plain text
pprint(expr)         # 2-D unicode art
latex(expr)          # for a paper
srepr(expr)          # the tree, exactly

ccode(expr);  fcode(expr);  rust_code(expr)
octave_code(expr);  mathematica_code(expr);  jscode(expr)
pycode(expr)

cse([expr1, expr2])  # common subexpression elimination

Tiplatex(Integral(x**2, (x, 0, 1))) gives you paper-ready LaTeX. This is, for many people, the single most-used feature in the library.

Interval(0, 1)                        # [0, 1]
Interval.open(0, 1)                   # (0, 1)     also .Lopen / .Ropen
FiniteSet(1, 2, 3)                    # {1, 2, 3}

Union(Interval(0,1), FiniteSet(3))
Intersection(Interval(0,2), Interval(1,3))
Complement(S.Reals, FiniteSet(0))

S.Reals;  S.Integers;  S.Naturals;  S.Complexes;  S.EmptySet
imageset(Lambda(n, 2*n), S.Integers)  # {..., -2, 0, 2, ...}
ConditionSet(x, x**2 > 4, S.Reals)    # "the x in R with ..."

• s.contains(2) / 2 in s
  Membership.
• s.is_subset(t) / s.measure
  Comparison and size.
• ImageSet
  How solveset expresses infinitely many roots.

• And(p, q) / Or / Not / Xor
  Also &, |, ~, ^.
• Implies(p, q) / Equivalent(p, q)
  Implication, biconditional.
• simplify_logic(expr)
  Minimize a boolean expression.
• to_cnf / to_dnf
  Normal forms.
• satisfiable(expr)
  A SAT solver, built in. Returns a model or False.

Gotcha^ is Xor, not exponentiation. x ^ y is a boolean, and it will not error — it will just be wrong. Powers are **.

p = Poly(x**3 - 2*x + 1, x)
p.degree();  p.all_coeffs();  p.LC();  p.eval(2)

factor_list(x**2 - 1)
gcd(x**2 - 1, x - 1);  lcm(x, x**2)
div(x**2 + 1, x - 1)                   # (quotient, remainder)
resultant(f, g);  discriminant(p)
groebner([x**2 - 1, y - x], x, y)
minimal_polynomial(sqrt(2) + sqrt(3), x)
interpolate([(1,1), (2,4), (3,9)], x)  # -> x**2

Poly is a dense internal representation — much faster than the generic expression tree for heavy polynomial work. Convert with .as_expr().

• isprime(97) / nextprime / prime(5)
  Primality and enumeration.
• factorint(360)
  {2: 3, 3: 2, 5: 1}
• primerange(1, 20) / divisors(28) / totient(10)
  The usual suspects.
• mod_inverse(3, 11) / gcd / lcm
  Modular arithmetic.
• binomial / factorial / fibonacci / catalan / bell
  Exact, arbitrarily large.
• from sympy.combinatorics import Permutation
  Permutation groups, Polyhedra, Prufer sequences.

from sympy.stats import Normal, Die, E, variance, density, P
X = Normal('X', 0, 1)
E(X);  variance(X);  density(X)(x);  P(X > 1)  # SYMBOLIC probability

from sympy.geometry import Point, Line, Circle, Triangle
Circle(Point(0,0), 1).intersection(Line(Point(0,0), Point(1,0)))

from sympy.physics.units import convert_to, meter, second, km, hour
convert_to(5*meter/second, km/hour)            # 18*km/hour

Tipsympy.stats gives you the exact density and moments, symbolically. It is a genuinely different tool from scipy.stats.

• lambdify, never subs-in-a-loop
  5,780× on the measured example. This is lever one, two and three.
• cse(exprs)
  Common subexpression elimination before you generate code.
• Poly(expr, x)
  For polynomial-heavy work — a dense representation beats the tree.
• xreplace over subs
  Exact structural swap: no sympify, no pattern matching. Much faster.
• Build unevaluated
  Integral, Sum, Derivative, then .doit() once.
• Declare assumptions early
  A positive=True symbol prunes whole branches of the simplifier.
• Avoid simplify() in a loop
  It is the most expensive call in the library.

TipProfile before you optimise: count_ops(expr) tells you how big the tree really is. Expression swell is usually the culprit, not SymPy itself.

• Symbols must be declared
  x + 1 → NameError until x = symbols('x').
• Symbols ≠ Python variables
  x = symbols('x'); e = x+1; x = 2 → e is still x + 1.
• <code>=</code> is not equality
  It is assignment. Eq(a, b) builds an equation.
• <code>==</code> is structural
  (x+1)**2 == x**2+2*x+1 → False. Use simplify(a-b)==0.
• <code>^</code> is Xor
  Powers are **. x ^ y silently gives you a boolean.
• <code>/</code> on two ints is a float
  x + 1/2 → x + 0.5. Use Rational(1,2) or S(1)/2.
• No implicit multiplication
  3x is a SyntaxError. Python, not SymPy.

Source: the official Tutorial → Gotchas and Explanations → Gotchas and Pitfalls pages.

• <code>solve(sin(x), x)</code>
  Returns [0, pi] — only two of infinitely many roots.
• Sequential <code>subs</code>
  (x+y).subs([(x,y),(y,2)]) → 4, not y+2. Pass simultaneous=True.
• <code>integrate</code> returns Piecewise
  integrate(x**y, x) → a Piecewise, not an expression.
• <code>if x &gt; 0:</code>
  TypeError: cannot determine truth value of Relational.
• <code>Symbol('x')</code> vs <code>Symbol('x', positive=True)</code>
  Different symbols. They will not cancel against each other.
• <code>A * B</code> on Matrices
  Matrix product — the opposite of NumPy's *.
• <code>Matrix</code> is mutable
  Therefore unhashable. Use ImmutableMatrix in expressions.
• <code>expand_log</code> does nothing
  Needs force=True or positive assumptions.
• <code>.is_positive</code> is <code>None</code>
  Three-valued logic. None means unknown, not False.
• Expressions are immutable
  expr.subs(x, 2) returns a new tree. Assign it.