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

Common Simplification Functions

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

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.

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

1. 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))

1. 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.

1. 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.

1. 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.

1. 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

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