SymPy — Symbolic Mathematics

Overview

SymPy is a Python library for symbolic mathematics that performs exact computation using mathematical symbols rather than numerical approximations. It covers algebra, calculus, equation solving, linear algebra, physics, and code generation — all within pure Python with no external dependencies.

When to Use

• Solving equations symbolically (algebraic, systems, differential equations)
• Performing calculus operations (derivatives, integrals, limits, series expansions)
• Simplifying and manipulating algebraic expressions
• Working with matrices symbolically (eigenvalues, determinants, decompositions)
• Converting symbolic expressions to fast numerical functions (lambdify → NumPy)
• Generating code from math expressions (C, Fortran, LaTeX)
• Needing exact results (e.g., sqrt(2) not 1.414...)
• For numerical computing (array operations, linear algebra on data), use numpy/scipy
• For statistical modeling (regression, hypothesis testing), use statsmodels

Prerequisites

pip install sympy
# Optional for numerical evaluation:
pip install numpy matplotlib

SymPy is pure Python — no compiled dependencies, installs everywhere.

Quick Start

from sympy import symbols, solve, diff, integrate, sqrt, pi

x = symbols('x')

# Solve equation
print(solve(x**2 - 5*x + 6, x))          # [2, 3]

# Derivative
print(diff(x**3 + 2*x, x))                # 3*x**2 + 2

# Integral
print(integrate(x**2, (x, 0, 1)))         # 1/3

# Exact arithmetic
print(sqrt(8))                              # 2*sqrt(2)
print(pi.evalf(30))                        # 3.14159265358979323846264338328

Core API

1. Symbols and Expressions

Create symbolic variables and manipulate expressions.

from sympy import symbols, Symbol, Rational, S, oo, pi, E, I
from sympy import simplify, expand, factor, collect, cancel, trigsimp

# Define symbols
x, y, z = symbols('x y z')

# With assumptions (improve simplification)
n = symbols('n', integer=True)
t = symbols('t', positive=True, real=True)
from sympy import sqrt
print(sqrt(t**2))   # t (not Abs(t), because t is positive)

# Exact fractions (avoid floats!)
expr = Rational(1, 3) * x + S(1)/7
print(expr)  # x/3 + 1/7

# Simplification
print(simplify(x**2 + 2*x + 1))          # (x + 1)**2
print(expand((x + 1)**3))                 # x**3 + 3*x**2 + 3*x + 1
print(factor(x**3 - x))                   # x*(x - 1)*(x + 1)
print(collect(x*y + x - 3 + 2*x**2 - z*x**2, x))  # x**2*(2 - z) + x*(y + 1) - 3

2. Calculus

Derivatives, integrals, limits, and series.

from sympy import symbols, diff, integrate, limit, series, oo, sin, cos, exp, log

x = symbols('x')

# Derivatives
print(diff(sin(x**2), x))                 # 2*x*cos(x**2)
print(diff(x**4, x, 3))                   # 24*x (third derivative)

# Partial derivatives
x, y = symbols('x y')
f = x**2 * y**3
print(diff(f, x, y))                      # 6*x*y**2

# Integrals
x = symbols('x')
print(integrate(x**2, x))                 # x**3/3 (indefinite)
print(integrate(exp(-x**2), (x, -oo, oo)))  # sqrt(pi) (Gaussian)
print(integrate(x * exp(-x), (x, 0, oo)))  # 1

# Limits
print(limit(sin(x)/x, x, 0))             # 1
print(limit((1 + 1/x)**x, x, oo))        # E

# Taylor series
print(series(exp(x), x, 0, 5))           # 1 + x + x**2/2 + x**3/6 + x**4/24 + O(x**5)

3. Equation Solving

Algebraic, transcendental, and differential equations.

from sympy import symbols, solve, solveset, Eq, S, linsolve, nonlinsolve, Function, dsolve

x, y = symbols('x y')

# Single equation
print(solve(x**2 - 4, x))                 # [-2, 2]
print(solveset(x**2 - 4, x, S.Reals))     # {-2, 2}

# System of linear equations
print(linsolve([x + y - 5, 2*x - y - 1], x, y))  # {(2, 3)}

# System of nonlinear equations
print(nonlinsolve([x**2 + y - 4, x + y**2 - 4], x, y))

# Differential equation: y'' + y = 0
f = Function('f')
ode = f(x).diff(x, 2) + f(x)
print(dsolve(ode, f(x)))                  # Eq(f(x), C1*sin(x) + C2*cos(x))

# With initial conditions
from sympy import Derivative
ics = {f(0): 1, f(x).diff(x).subs(x, 0): 0}
print(dsolve(ode, f(x), ics=ics))         # Eq(f(x), cos(x))

4. Matrices and Linear Algebra

Symbolic matrix operations.

from sympy import Matrix, eye, zeros, ones, diag, symbols

# Create matrices
M = Matrix([[1, 2], [3, 4]])
print(f"Det: {M.det()}")                   # -2
print(f"Inverse:\n{M**-1}")

# Symbolic matrices
a, b = symbols('a b')
M = Matrix([[a, b], [b, a]])
print(f"Eigenvalues: {M.eigenvals()}")     # {a - b: 1, a + b: 1}

# Eigenvectors and diagonalization
eigendata = M.eigenvects()
# [(eigenval, multiplicity, [eigenvectors]), ...]
P, D = M.diagonalize()
print(f"M = P*D*P^-1")

# Solve linear system Ax = b
A = Matrix([[1, 2], [3, 4]])
b = Matrix([5, 6])
x = A.solve(b)
print(f"Solution: {x.T}")

# Matrix calculus
t = symbols('t')
M_t = Matrix([[t, t**2], [1, t]])
print(f"dM/dt:\n{M_t.diff(t)}")

5. Code Generation

Convert symbolic expressions to fast numerical functions or compiled code.

import numpy as np
from sympy import symbols, lambdify, sin, exp, ccode, fcode, latex

x, y = symbols('x y')
expr = sin(x) * exp(-x**2 / 2)

# lambdify: symbolic → fast NumPy function
f = lambdify(x, expr, 'numpy')
x_vals = np.linspace(-5, 5, 1000)
y_vals = f(x_vals)
print(f"Shape: {y_vals.shape}, Max: {y_vals.max():.4f}")

# Multi-variable lambdify
expr2 = x**2 + y**2
f2 = lambdify((x, y), expr2, 'numpy')
print(f"f(3, 4) = {f2(3, 4)}")            # 25

# C code generation
print(ccode(expr))                          # sin(x)*exp(-1.0/2.0*pow(x, 2))

# Fortran code generation
print(fcode(expr))

# LaTeX output
print(latex(expr))                          # \sin{\left(x \right)} e^{- \frac{x^{2}}{2}}

6. Physics Module

Classical mechanics, vector analysis, and units.

from sympy import symbols, cos, sin, Function
from sympy.physics.mechanics import dynamicsymbols, LagrangesMethod, Particle, Point, ReferenceFrame
from sympy.physics.vector import dot, cross

# Vector analysis
N = ReferenceFrame('N')
v1 = 3*N.x + 4*N.y + 0*N.z
v2 = 1*N.x + 0*N.y + 2*N.z
print(f"Dot: {dot(v1, v2)}")               # 3
print(f"Cross: {cross(v1, v2)}")            # 8*N.x - 6*N.y - 4*N.z

# Simple pendulum via Lagrangian mechanics
q = dynamicsymbols('q')       # Generalized coordinate (angle)
m, g, l = symbols('m g l', positive=True)
T = Rational(1, 2) * m * (l * q.diff())**2          # Kinetic energy
V = m * g * l * (1 - cos(q))                         # Potential energy
L = T - V                                            # Lagrangian
print(f"Lagrangian: {L}")

Key Concepts

Exact vs Numerical Arithmetic

from sympy import Rational, S, sqrt, pi

# WRONG: introduces floating-point error
expr_bad = 0.5 * x           # Float 0.5, loses exactness

# CORRECT: exact symbolic arithmetic
expr_good = Rational(1, 2) * x    # Exact 1/2
expr_good = S(1)/2 * x            # Alternative exact syntax
expr_good = x / 2                  # Also exact

# Numerical evaluation when needed
print(sqrt(2).evalf())         # 1.41421356237310
print(pi.evalf(50))            # 50 digits of precision

Solver Selection Guide

Solver	Use When	Returns
solve(eq, x)	General purpose, legacy	List of solutions
solveset(eq, x, domain)	Algebraic equations (preferred)	Set (may be infinite)
linsolve(system, vars)	Linear systems	FiniteSet of tuples
nonlinsolve(system, vars)	Nonlinear systems	FiniteSet of tuples
dsolve(ode, f(x))	Ordinary differential equations	Equality (Eq)
nsolve(eq, x0)	Numerical root finding	Float approximation

Common Simplification Functions

Function	Does	Example
simplify()	General simplification (slow, tries everything)	sin(x)**2 + cos(x)**2 → 1
expand()	Distribute multiplication	(x+1)**2 → x**2+2*x+1
factor()	Factor into irreducibles	x**2-1 → (x-1)*(x+1)
collect()	Group by variable	Collect terms in x
cancel()	Cancel common factors in fractions	(x**2-1)/(x-1) → x+1
trigsimp()	Simplify trig expressions	Faster than simplify for trig
powsimp()	Simplify powers/exponentials	Combine x**a * x**b

Common Workflows

Workflow: Symbolic-to-Numeric Pipeline

from sympy import symbols, diff, integrate, lambdify, sin, cos
import numpy as np
import matplotlib.pyplot as plt

x = symbols('x')

# 1. Define expression symbolically
f_expr = sin(x) * cos(x)**2

# 2. Symbolic operations
f_prime = diff(f_expr, x)
F_expr = integrate(f_expr, x)
print(f"f(x) = {f_expr}")
print(f"f'(x) = {f_prime}")
print(f"F(x) = {F_expr}")

# 3. Convert to fast numerical functions
f_num = lambdify(x, f_expr, 'numpy')
f_prime_num = lambdify(x, f_prime, 'numpy')
F_num = lambdify(x, F_expr, 'numpy')

# 4. Evaluate and plot
x_vals = np.linspace(0, 2*np.pi, 500)
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
axes[0].plot(x_vals, f_num(x_vals)); axes[0].set_title('f(x)')
axes[1].plot(x_vals, f_prime_num(x_vals)); axes[1].set_title("f'(x)")
axes[2].plot(x_vals, F_num(x_vals)); axes[2].set_title('F(x)')
plt.tight_layout()
plt.savefig('symbolic_pipeline.png', dpi=150)
print("Saved symbolic_pipeline.png")

Workflow: Solve and Verify

from sympy import symbols, solve, simplify, Eq, sqrt

x = symbols('x')

# 1. Define equation
equation = x**3 - 6*x**2 + 11*x - 6

# 2. Solve symbolically
solutions = solve(equation, x)
print(f"Solutions: {solutions}")            # [1, 2, 3]

# 3. Verify each solution
for sol in solutions:
    result = simplify(equation.subs(x, sol))
    assert result == 0, f"Solution {sol} failed!"
    print(f"  x={sol}: f(x) = {result} ✓")

# 4. Factor the polynomial
from sympy import factor
print(f"Factored: {factor(equation)}")      # (x - 1)*(x - 2)*(x - 3)

Workflow: ODE System Analysis

1. Define the ODE using Function and dsolve()
2. Solve symbolically; apply initial conditions with ics={} parameter
3. Convert solution to numerical function with lambdify()
4. Plot the solution trajectory with matplotlib

Key Parameters

Parameter	Function	Default	Options	Effect
domain	solveset()	S.Complexes	S.Reals, S.Integers	Restrict solution domain
force	simplify()	False	True/False	Aggressive simplification
n	diff(expr, x, n)	1	1–∞	Order of derivative
Precision	evalf(n)	15	1–1000+	Digits of numerical precision
Backend	lambdify()	"math"	"numpy", "scipy", "mpmath"	Numerical backend for evaluation
rational	nsimplify()	True	True/False	Find exact rational approximation

Best Practices

1. Always use Rational() or S() for fractions: 0.5 * x introduces floats that break exact computation. Use Rational(1, 2) * x or S(1)/2 * x.

2. Add assumptions to symbols: symbols('x', positive=True) enables simplifications like sqrt(x**2) → x. Without assumptions, SymPy must handle the general complex case.

3. Use lambdify for numerical evaluation, not subs().evalf(): subs/evalf in a loop is 100-1000x slower than a single lambdify call.

   # Slow: [expr.subs(x, v).evalf() for v in values]
   # Fast: f = lambdify(x, expr, 'numpy'); f(np.array(values))
   

4. Anti-pattern — using simplify() as default: simplify() is slow because it tries many strategies. Use specific functions (factor, expand, trigsimp) when you know the desired form.

5. Prefer solveset over solve for algebraic equations: solveset returns proper mathematical sets and handles edge cases better. solve is legacy but still useful for general cases.

6. Anti-pattern — solving symbolically when numerical is sufficient: For equations with no closed-form solution, use nsolve(eq, x0) for numerical root finding instead of waiting for solve to fail.

7. Use init_printing() in Jupyter for readable output: from sympy import init_printing; init_printing() enables LaTeX rendering in notebooks.

Common Recipes

Recipe: Generate LaTeX Documentation

from sympy import symbols, Integral, Eq, latex, sqrt, pi

x = symbols('x')
integral = Integral(x**2 * sqrt(1 - x**2), (x, 0, 1))
result = integral.doit()

print(f"$$ {latex(integral)} = {latex(result)} $$")
# $$ \int\limits_{0}^{1} x^{2} \sqrt{1 - x^{2}}\, dx = \frac{\pi}{16} $$

Recipe: Parametric ODE Solution

from sympy import symbols, Function, dsolve, Eq, exp, lambdify
import numpy as np

x = symbols('x')
k, A = symbols('k A', positive=True)
f = Function('f')

# Solve y' = -ky with y(0) = A
ode = Eq(f(x).diff(x), -k * f(x))
solution = dsolve(ode, f(x), ics={f(0): A})
print(f"Solution: {solution}")              # f(x) = A*exp(-k*x)

# Evaluate for specific parameters
f_num = lambdify((x, k, A), solution.rhs, 'numpy')
x_vals = np.linspace(0, 5, 100)
y_vals = f_num(x_vals, k=0.5, A=10)
print(f"y(5) = {y_vals[-1]:.4f}")

Recipe: Symbolic Matrix Decomposition

from sympy import Matrix, symbols, pprint

a, b, c, d = symbols('a b c d')
M = Matrix([[a, b], [c, d]])

# Characteristic polynomial
lam = symbols('lambda')
char_poly = M.charpoly(lam)
print(f"Characteristic polynomial: {char_poly.as_expr()}")

# Eigenvalues (symbolic)
eigenvals = M.eigenvals()
print(f"Eigenvalues: {eigenvals}")

# Determinant and trace
print(f"det(M) = {M.det()}")               # a*d - b*c
print(f"tr(M) = {M.trace()}")              # a + d

Troubleshooting

Problem	Cause	Solution
NameError: name 'x' is not defined	Symbol not created	Define with x = symbols('x') before use
Unexpected float results	Using 0.5 instead of Rational(1,2)	Use Rational() or S() for exact fractions
simplify() very slow	Trying all strategies on complex expr	Use specific function: factor(), expand(), trigsimp()
solve() returns empty list	No closed-form solution exists	Use nsolve(eq, x0) for numerical approximation
sqrt(x**2) returns sqrt(x**2) not x	No assumption on x	Define x = symbols('x', positive=True)
lambdify wrong results	Expression has SymPy-specific functions	Specify backend: lambdify(x, expr, 'numpy') or 'scipy'
NotImplementedError in dsolve	ODE type not supported	Try numerical ODE solver (scipy odeint) instead

Related Skills

• matplotlib-scientific-plotting — plot symbolic results after lambdify conversion
• statsmodels-statistical-modeling — statistical inference; use when you need p-values, not exact algebra
• matlab-scientific-computing — MATLAB alternative for numerical (not symbolic) computing

References

• SymPy documentation — official API reference and tutorials
• SymPy tutorial — introduction for new users
• SymPy GitHub — source code, examples, and issues
• Meurer et al. (2017) "SymPy: symbolic computing in Python" — PeerJ Computer Science