How To Use SymPy for Symbolic Mathematics in Python

Intermediate

Have you ever solved an algebra problem only to get a decimal approximation when you wanted the exact symbolic answer? Python’s SymPy library solves this problem by treating mathematics the way a mathematician does — symbolically. Instead of computing pi as 3.14159..., SymPy keeps it as the exact symbol pi. Instead of approximating a square root, it returns sqrt(2). This is symbolic computation, and it transforms Python into a full-featured computer algebra system.

SymPy is a pure Python library — no C extensions, no compiled code — so installation is straightforward with pip install sympy. It works alongside NumPy and SciPy but solves a different problem: those libraries compute numerical answers fast, while SymPy computes exact symbolic answers. You can use SymPy to solve equations, expand polynomials, compute derivatives and integrals, factor expressions, and even generate LaTeX output for publication-quality math.

In this article, we will cover the fundamentals of SymPy from the ground up: how to define symbolic variables, simplify and expand expressions, solve equations, compute limits and derivatives, evaluate integrals, and apply SymPy to a practical calculus problem. By the end, you will be able to use Python as a complete algebra and calculus tool.

SymPy Quick Example

Here is a self-contained example that shows SymPy solving a quadratic equation and computing a derivative — both with exact results:

# quick_sympy.py
from sympy import symbols, solve, diff, expand

x = symbols('x')

# Solve a quadratic equation
equation = x**2 - 5*x + 6
solutions = solve(equation, x)
print("Solutions:", solutions)

# Compute the derivative of x^3 - 2x + 1
expr = x**3 - 2*x + 1
derivative = diff(expr, x)
print("Derivative:", derivative)

# Expand a factored expression
expanded = expand((x + 2) * (x - 3))
print("Expanded:", expanded)

Output:

Solutions: [2, 3]
Derivative: 3*x**2 - 2
Expanded: x**2 - x - 6

These are exact symbolic results. solve() found the roots of the quadratic as integers, not decimals. diff() computed the derivative using the power rule. expand() multiplied out the factors. Every result is a SymPy expression that can be further manipulated, printed as LaTeX, or evaluated numerically.

What Is SymPy and Why Use It?

SymPy is a computer algebra system (CAS) written entirely in Python. A CAS is software that manipulates mathematical expressions in symbolic form, the way a human mathematician writes on paper, rather than computing numerical approximations. SymPy’s closest analogues are Mathematica and Maple, but SymPy is free, open source, and integrates naturally with the Python data science ecosystem.

Library	Purpose	Result Type	Best For
NumPy	Numerical computation	Float approximation	Fast array math
SciPy	Scientific algorithms	Float approximation	Optimization, stats
SymPy	Symbolic computation	Exact expression	Algebra, calculus proofs

Use SymPy when you need exact answers — factoring polynomials, solving equations analytically, or deriving formulas before implementing them numerically.

Defining Symbolic Variables

The foundation of SymPy is the Symbol class. Before using any variable in a symbolic expression, you declare it with symbols(). This tells SymPy “this letter stands for a mathematical unknown, not a Python value.”

# symbols_demo.py
from sympy import symbols, Symbol

# Single symbol
x = symbols('x')

# Multiple symbols at once
a, b, c = symbols('a b c')

# Symbol with assumptions (positive, real, integer, etc.)
n = symbols('n', integer=True, positive=True)
t = symbols('t', real=True)

print(type(x))        # 
print(x + x)              # 2*x
print(x * x)         # x**2
print(a + b + a)      # 2*a + b

Output:

<class 'sympy.core.symbol.Symbol'>
2*x
x**2
2*a + b

Assumptions like positive=True or integer=True help SymPy simplify expressions correctly. For example, sqrt(x**2) only simplifies to x when SymPy knows x is non-negative.

Simplifying and Manipulating Expressions

Once you have symbols, you can build expressions and simplify them. SymPy provides simplify(), expand(), factor(), cancel(), and collect() — each targeting different transformation goals.

# simplify_demo.py
from sympy import symbols, simplify, expand, factor, cancel, trigsimp, sin, cos

x, y = symbols('x y')

# simplify -- general-purpose simplification
expr1 = (x**2 - 1) / (x - 1)
print("cancel:", cancel(expr1))           # x + 1

# expand -- distribute multiplication
expr2 = (x + y)**3
print("expand:", expand(expr2))           # x**3 + 3*x**2*y + 3*x*y**2 + y**3

# factor -- reverse of expand
expr3 = x**3 + 3*x**2*y + 3*x*y**2 + y**3
print("factor:", factor(expr3))           # (x + y)**3

# trigsimp -- simplify trig expressions
trig_expr = sin(x)**2 + cos(x)**2
print("trigsimp:", trigsimp(trig_expr))   # 1

Output:

cancel: x + 1
expand: x**3 + 3*x**2*y + 3*x*y**2 + y**3
factor: (x + y)**3
trigsimp: 1

Notice that cancel() correctly identified that (x**2 - 1)/(x - 1) simplifies to x + 1 by canceling the common factor (x - 1). This is exact symbolic cancellation — floating-point arithmetic would have struggled near x = 1.

Solving Equations

The solve() function finds the values of unknowns that satisfy an equation. You can solve for one variable in terms of others, solve systems of equations, and even solve inequalities. Pass the expression (or a list of expressions) and the variable(s) to solve for.

# solve_demo.py
from sympy import symbols, solve, Eq, sqrt
	x, y = symbols('x y')

# Solve x^2 = 9
solutions = solve(x**2 - 9, x)
print("x^2 = 9:", solutions)           # [-3, 3]

# Use Eq() for equations with both sides
eq = Eq(2*x + 3, 11)
print("2x + 3 = 11:", solve(eq, x))      # [4]

# System of equations
eq1 = Eq(x + y, 7)
eq2 = Eq(2*x - y, 2)
print("System:", solve([eq1, eq2], [x, y]))  # {x: 3, y: 4}

# Solve quadratic formula symbolically
a, b, c = symbols('a b c')
quad = a*x**2 + b*x + c
print("Quadratic formula:", solve(quad, x))

Output:

x^2 = 9: [-3, 3]
2x + 3 = 11: [4]
System: {x: 3, y: 4}
Quadratic formula: [(-b - sqrt(-4*a*c + b**2))/(2*a), (-b + sqrt(-4*a*c + b**2))/(2*a)]

The last result is the closed-form quadratic formula — SymPy returns both roots as exact symbolic expressions. No numerical approximation, no round-off error. This is the kind of thing you’d otherwise look up in a reference; SymPy derives it from a*x**2 + b*x + c = 0 in one call.

Calculus: Derivatives, Integrals, Limits

SymPy’s calculus toolbox is what makes it indispensable for physics, engineering, and machine-learning gradient work. The four core functions you’ll use most are diff() (derivatives), integrate() (integrals), limit() (limits), and series() (Taylor expansion).

# calculus_demo.py
from sympy import symbols, diff, integrate, limit, series, sin, cos, oo, exp

x = symbols('x')

# Derivatives
f = x**3 + 2*x**2 - 5*x + 1
print("f(x)    =", f)
print("f'(x)   =", diff(f, x))         # first derivative
print("f''(x)  =", diff(f, x, 2))      # second derivative

# Partial derivatives
y = symbols('y')
g = x**2 * y + y**3
print("dg/dx  =", diff(g, x))          # treat y as constant
print("dg/dy  =", diff(g, y))

# Indefinite integral (antiderivative)
print("∫(2x + 3) dx =", integrate(2*x + 3, x))

# Definite integral
print("∫₀^π sin(x) dx =", integrate(sin(x), (x, 0, 3.14159)))

# Limits
print("lim x→0 sin(x)/x =", limit(sin(x)/x, x, 0))
print("lim x→∞ 1/x      =", limit(1/x, x, oo))

# Taylor series — first 5 terms of e^x around 0
print("e^x series:", series(exp(x), x, 0, 5))

Output:

f(x)    = x**3 + 2*x**2 - 5*x + 1
f'(x)   = 3*x**2 + 4*x - 5
f''(x)  = 6*x + 4
dg/dx  = 2*x*y
dg/dy  = x**2 + 3*y**2
∫(2x + 3) dx = x**2 + 3*x
∫₀^π sin(x) dx = 1.99999999...
lim x→0 sin(x)/x = 1
lim x→∞ 1/x      = 0
e^x series: 1 + x + x**2/2 + x**3/6 + x**4/24 + O(x**5)

oo is SymPy’s infinity. Notice limit(sin(x)/x, x, 0) returns exactly 1, not a numerical estimate — SymPy applies L’Hôpital’s rule symbolically. The Taylor series O(x**5) is the “big-O” remainder, exact to the term you asked for.

Linear Algebra with Matrices

SymPy ships with a Matrix class that works just like NumPy’s arrays but holds symbolic entries. You can compute determinants, inverses, eigenvalues, and reduced row echelon form symbolically:

# matrix_demo.py
from sympy import Matrix, symbols, eye

a, b = symbols('a b')

M = Matrix([
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 10],   # 10 instead of 9 so the matrix is invertible
])

print("Determinant:", M.det())
print("Inverse:")
print(M.inv())
print("Eigenvalues:", M.eigenvals())
print("Rank:", M.rank())

# Symbolic 2x2 matrix
S = Matrix([[a, b], [b, a]])
print("S squared:", S * S)
print("S inverse:", S.inv())

Output:

Determinant: -3
Inverse:
Matrix([[-2/3, -4/3, 1], [-2/3, 11/3, -2], [1, -2, 1]])
Eigenvalues: {16/3 + ...: 1, 16/3 - ...: 1, -1: 1}
Rank: 3
S squared: Matrix([[a**2 + b**2, 2*a*b], [2*a*b, a**2 + b**2]])
S inverse: Matrix([[a/(a**2 - b**2), -b/(a**2 - b**2)], [-b/(a**2 - b**2), a/(a**2 - b**2)]])

For numerical linear algebra at scale, use NumPy or SciPy — SymPy is slower because every cell carries its symbolic structure. SymPy shines when you want the closed form, not a number.

From Symbolic to Numerical: evalf() and lambdify()

Symbolic math is great, but eventually you want numbers. Two functions bridge the gap:

evalf() — evaluate a symbolic expression to a numeric value with arbitrary precision:

from sympy import pi, sqrt, E

print(pi.evalf())                  # 3.14159265358979
print(pi.evalf(50))                # 50 digits of pi
print(sqrt(2).evalf())             # 1.41421356237310
print((E**2).evalf())              # 7.38905609893065

lambdify() — convert a symbolic expression into a fast NumPy/SciPy-compatible function. This is the bridge to numerical work, plotting, and ML:

from sympy import symbols, lambdify, sin, cos
import numpy as np

x = symbols('x')
expr = sin(x) + cos(x) * x**2

# Build a numeric function from the symbolic expression
f = lambdify(x, expr, modules='numpy')

# Now call it like any vectorized function
xs = np.linspace(0, 6, 7)
print(f(xs))

lambdify compiles the symbolic expression into Python source and then turns it into a callable. The first call is the slowest because of compilation; subsequent calls are pure NumPy speed. This pattern is the standard way to derive a gradient symbolically with SymPy, then plug it into a numerical optimizer like scipy.optimize.

Common Pitfalls and Gotchas

• Forgetting to declare symbols. solve(x**2 - 1, x) fails with NameError if x isn’t declared first via symbols('x'). SymPy doesn’t treat Python variable names as symbolic — you have to opt in.
• Integer division surprises. 1/2 in Python 3 is 0.5 (a float), but inside SymPy expressions you often want the exact fraction 1/2. Wrap literals with Rational(1, 2) when exactness matters: integrate(x, (x, 0, Rational(1, 2))).
• Confusing = and Eq. solve(x**2 = 9, x) is a Python syntax error. Use either solve(x**2 - 9, x) (move everything to one side, SymPy assumes = 0) or solve(Eq(x**2, 9), x) (build an explicit equation object).
• Slowness on large expressions. SymPy is exact, not fast. Expressions with hundreds of nested symbols can take minutes to simplify. For hot loops, lambdify the expression to NumPy and run the loop there.
• Assumptions matter. sqrt(x**2) returns sqrt(x**2) by default because SymPy doesn’t know if x is positive. Declare it: x = symbols('x', positive=True) and now sqrt(x**2) simplifies to x.

FAQ

Q: When should I use SymPy vs NumPy?
A: Use SymPy when you need exact answers, closed-form expressions, derivatives, or to prove an identity. Use NumPy when you have numeric arrays and need speed. They complement each other — SymPy derives the math, NumPy runs the math.

Q: Why does integrate() return an unevaluated Integral?
A: SymPy returns the unevaluated form when it can’t find a closed-form antiderivative. Some integrals genuinely don’t have an elementary antiderivative (Gaussian, error function, etc.). Try a definite integral with explicit bounds, or use nintegrate for numerical integration via mpmath.

Q: How do I plot a SymPy expression?
A: Either use SymPy’s built-in sympy.plotting.plot() (good for quick plots), or lambdify the expression and pass it to matplotlib. The matplotlib path gives you full styling control.

Q: Can SymPy solve differential equations?
A: Yes — dsolve() handles ordinary differential equations symbolically. It works for linear ODEs, separable equations, and many classical forms. For PDEs or messy nonlinear ODEs, fall back to scipy.integrate.solve_ivp with a lambdified RHS.

Q: Is SymPy thread-safe?
A: Mostly yes for read-only operations, but the assumption system has global state. If you’re running symbolic math in multiple threads, give each thread its own set of symbols and don’t share simplification caches.

Wrapping Up

SymPy gives Python the same symbolic-math power that Mathematica and Maple have offered for decades, in pure Python with a clean API. Start with symbols(), simplify(), solve(), diff(), and integrate() — those five functions cover 80% of everyday symbolic work. When you need a number, reach for evalf() or lambdify() and hand off to NumPy.

The SymPy documentation has the complete reference, including the more specialized modules (statistics, geometry, combinatorics, physics, quantum). For tutorials on related Python topics, see below.

Related Python Tutorials

• How To Use SciPy for Scientific Computing in Python
• How To Use NumPy Arrays for Scientific Computing in Python
• How To Use Python Dask for Parallel and Distributed Computing