SymPy Symbolic Math — Deep Dive

How SymPy Represents Expressions

SymPy builds expressions as tree structures (DAGs). The expression x**2 + 2*x + 1 becomes:

    Add
   / | \
Pow  Mul  1
/ \  / \
x  2 2  x

Every node is an immutable Python object inheriting from Basic. This design enables pattern matching, substitution, and transformation via tree traversal.

Key Internal Classes

• Symbol — Leaf node representing an unknown variable.
• Number — Leaf node for constants (Integer, Rational, Float).
• Add, Mul, Pow — Internal nodes for operations.
• Function — Nodes for sin, cos, exp, and user-defined functions.

You can inspect any expression’s tree:

from sympy import symbols, srepr
x = symbols('x')
print(srepr(x**2 + 1))
# Add(Pow(Symbol('x'), Integer(2)), Integer(1))

Advanced Solving

Systems of Nonlinear Equations

from sympy import symbols, solve, sin, cos

x, y = symbols('x y')
solutions = solve([x**2 + y**2 - 1, y - sin(x)], [x, y])

SymPy attempts analytical solutions. For transcendental systems like this one, it may return an implicit form or use nsolve for numerical approximation.

Differential Equations

from sympy import Function, dsolve, Eq, symbols

t = symbols('t')
y = Function('y')

# y'' + y = 0 (simple harmonic oscillator)
ode = Eq(y(t).diff(t, 2) + y(t), 0)
solution = dsolve(ode, y(t))
# y(t) = C1*sin(t) + C2*cos(t)

SymPy handles first and second-order ODEs, linear systems, and many special types (Bernoulli, Riccati, exact equations). For PDEs, support is more limited but includes separation of variables.

Matrix Operations

from sympy import Matrix, symbols

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

M.det()           # a*d - b*c
M.inv()           # Symbolic inverse
M.eigenvals()     # Eigenvalues as dict {eigenvalue: multiplicity}
M.diagonalize()   # (P, D) where M = P*D*P^(-1)

Symbolic matrices are invaluable for deriving closed-form solutions in control theory, quantum mechanics, and statistics.

Assumptions System

SymPy’s assumptions engine lets you declare properties of symbols, which affects simplification:

x = symbols('x', positive=True)
from sympy import sqrt

sqrt(x**2)  # x (not Abs(x), because x is positive)

Available assumptions include real, positive, integer, prime, finite, commutative, and many more. Without assumptions, SymPy must be conservative:

x = symbols('x')  # No assumptions
sqrt(x**2)  # sqrt(x**2) — can't simplify without knowing x's sign

Code Generation

SymPy can export expressions to optimized code in multiple languages.

Python/NumPy via lambdify

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

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

f = lambdify(x, expr, modules='numpy')
f(np.linspace(0, 2*np.pi, 1000))  # Fast NumPy evaluation

C Code Generation

from sympy.utilities.codegen import codegen

x, y = symbols('x y')
expr = x**2 + y**2 + 2*x*y

[(c_name, c_code), (h_name, h_code)] = codegen(
    ('distance_squared', expr), 'C', 'output'
)
print(c_code)

This generates a C function with proper headers. SymPy also supports Fortran, Rust, Julia, JavaScript, and Octave.

Optimized Expression Trees

from sympy import cse

x, y, z = symbols('x y z')
exprs = [x**2 + y**2, x**2 + z**2, x**2 + y**2 + z**2]

replacements, reduced = cse(exprs)
# replacements: [(x0, x**2), (x1, y**2), (x2, z**2)]
# reduced: [x0 + x1, x0 + x2, x0 + x1 + x2]

Common Subexpression Elimination (cse) identifies repeated computations — critical for generating efficient numerical code from complex symbolic expressions.

Performance Strategies

SymPy’s pure-Python implementation means it’s orders of magnitude slower than compiled CAS like Mathematica for large expressions. Practical strategies:

1. Minimize Expression Complexity

# Slow: operate on expanded form
expanded = expand((x + y + z)**20)  # 1771 terms

# Fast: keep factored form when possible
factored = (x + y + z)**20  # Single node

2. Use evaluate=False for Construction

from sympy import Add

# SymPy normally auto-simplifies during construction
# For building large expressions, defer evaluation:
expr = Add(*terms, evaluate=False)
# Then simplify once at the end

3. Cache Intermediate Results

from sympy import Symbol, diff, cache

x = Symbol('x')
expr = heavy_computation(x)

# SymPy caches internally, but you can also pickle results
import pickle
with open('expr_cache.pkl', 'wb') as f:
    pickle.dump(expr, f)

4. Use SymEngine for Speed-Critical Paths

symengine is a C++ symbolic math library with a Python wrapper that’s 100-1000x faster than SymPy for expression manipulation:

import symengine
x = symengine.Symbol('x')
expr = (x + 1)**100
expanded = symengine.expand(expr)  # Milliseconds vs seconds in SymPy

SymEngine is API-compatible with SymPy for core operations, making migration straightforward.

Real-World Applications

Physics: Deriving Equations of Motion

from sympy import symbols, Function, Lagrangian
from sympy.physics.mechanics import LagrangesMethod, dynamicsymbols

# Define generalized coordinates
q = dynamicsymbols('q')
m, g, l = symbols('m g l', positive=True)

# Lagrangian for a simple pendulum
T = m * l**2 * q.diff()**2 / 2  # Kinetic energy
V = -m * g * l * cos(q)          # Potential energy
L = T - V

Engineering: Transfer Function Analysis

from sympy import symbols, apart, inverse_laplace_transform

s, t = symbols('s t')

# Second-order system transfer function
H = 1 / (s**2 + 2*s + 5)

# Partial fraction decomposition
apart(H, s)

# Inverse Laplace transform (time domain response)
from sympy import inverse_laplace_transform
h_t = inverse_laplace_transform(H, s, t)

Generating Optimized Numerical Kernels

A complete workflow: derive a gradient symbolically, optimize it, and export:

from sympy import symbols, Matrix, exp, cse, lambdify

x1, x2, w1, w2, b = symbols('x1 x2 w1 w2 b')
inputs = Matrix([x1, x2])
weights = Matrix([w1, w2])

# Sigmoid model
z = weights.dot(inputs) + b
sigma = 1 / (1 + exp(-z))

# Gradient w.r.t. weights
grad = Matrix([sigma.diff(w) for w in [w1, w2, b]])

# Optimize
replacements, optimized = cse(grad)

# Export to fast NumPy function
grad_fn = lambdify([x1, x2, w1, w2, b], optimized, modules='numpy')

Gotchas

Equality testing — expr1 == expr2 checks structural equality, not mathematical equality. Use simplify(expr1 - expr2) == 0 or expr1.equals(expr2) for mathematical comparison.

Float vs Rational — sympy.Float(0.1) introduces floating-point representation. Use sympy.Rational(1, 10) for exact arithmetic.

Import conflicts — SymPy redefines sin, cos, sqrt, etc. Import carefully to avoid shadowing NumPy or math module functions.

Memory — Large symbolic expressions can consume significant memory. An expanded polynomial of degree 20 in 10 variables has millions of terms. Keep expressions factored until you need the expanded form.

One Thing to Remember

SymPy turns Python into a computer algebra system — derive exact symbolic results, then bridge to numerical code with lambdify or code generation, giving you mathematical rigor and computational speed in the same workflow.