Claude Code for SymPy: Symbolic Mathematics

SymPy is a pure-Python computer algebra system. pip install sympy. from sympy import *. Symbols: x, y, z = symbols("x y z"), n = symbols("n", integer=True, positive=True). Algebra: expand((x+y)**3), factor(x**3-1), simplify(sin(x)**2 + cos(x)**2), cancel((x**2-1)/(x-1)). Differentiation: diff(sin(x)*exp(x), x), diff(f(x,y), x, 2). Integration: integrate(x**2*exp(x), x), integrate(sin(x), (x, 0, pi)). Limits: limit(sin(x)/x, x, 0), limit(1/x, x, oo). Series: series(exp(x), x, 0, 6). Solve: solve(x**2 - 2, x), solve([x+y-2, x-y-0], [x, y]). Solveset: solveset(cos(x), x, domain=S.Reals). ODE: f = Function("f"), dsolve(f(x).diff(x,2) - f(x), f(x)). Matrix: M = Matrix([[1,2],[3,4]]), M.det(), M.eigenvals(), M.eigenvects(). Numerics: N(pi, 50), N(sqrt(2), 20). LaTeX: latex(Rational(1,3) + x**2) → \\frac{1}{3} + x^{2}. Lambdify: f = lambdify([x, y], expr, "numpy"). Assumptions: with assuming(x>0): print(simplify(sqrt(x**2))). Substitution: expr.subs(x, 2), expr.subs([(x,1),(y,2)]). Statistics: from sympy.stats import Normal, E, variance, X = Normal("X", 0, 1), E(X**2). Claude Code generates SymPy calculus solvers, ODE pipelines, matrix algebra scripts, and numeric-symbolic hybrid workflows.

CLAUDE.md for SymPy

## SymPy Stack
- Version: sympy >= 1.12
- Symbols: symbols("x y", real=True) | Function("f") for undefined funcs
- Algebra: expand | factor | simplify | cancel | apart | together
- Calculus: diff(expr, x) | integrate(expr, (x, a, b)) | limit(expr, x, val)
- Solve: solve(eq, x) | solveset(eq, x, S.Reals) | dsolve(ode, f(x))
- Matrix: Matrix([[...]]) | .det() | .inv() | .eigenvals() | .rref()
- Numeric: N(expr, n_digits) | lambdify([x], expr, "numpy") for fast eval
- LaTeX: latex(expr) | init_printing() in Jupyter for pretty output

SymPy Symbolic Mathematics Pipeline

# math/sympy_pipeline.py — symbolic mathematics with SymPy
from __future__ import annotations
from sympy import *
from sympy.stats import Normal, E, variance, density
import numpy as np


# ── 1. Calculus ───────────────────────────────────────────────────────────────

def diff_and_integrate(expr, var, order: int = 1) -> dict:
    """
    Compute derivative and antiderivative of an expression.
    Returns dict with derivative, integral, and simplified forms.
    """
    derivative  = diff(expr, var, order)
    integral    = integrate(expr, var)
    return {
        "expression": expr,
        "derivative":  simplify(derivative),
        "integral":    simplify(integral),
    }


def compute_limit(expr, var, point, direction: str = "+") -> Any:
    """
    Compute limit of expr as var → point.
    direction: "+" (right), "-" (left), "+-" (two-sided)
    """
    if direction == "+-":
        lim_plus  = limit(expr, var, point, "+")
        lim_minus = limit(expr, var, point, "-")
        if lim_plus == lim_minus:
            return lim_plus
        return {"left": lim_minus, "right": lim_plus, "exists": False}
    return limit(expr, var, point, direction)


def taylor_series(
    expr,
    var,
    point: int | float = 0,
    order: int = 6,
    remove_order: bool = True,
) -> Any:
    """
    Compute Taylor/McLaurin series expansion up to given order.
    remove_order=True strips the O(x^n) remainder term.
    """
    s = series(expr, var, point, order)
    return s.removeO() if remove_order else s


def partial_derivatives(expr, variables: list) -> dict:
    """
    Compute all first-order partial derivatives.
    Returns dict {var_name: partial_derivative}.
    """
    return {str(v): simplify(diff(expr, v)) for v in variables}


def gradient_and_hessian(expr, variables: list) -> tuple:
    """
    Compute gradient vector and Hessian matrix.
    Returns (gradient Matrix, Hessian Matrix).
    """
    grad = Matrix([diff(expr, v) for v in variables])
    hess = Matrix([[diff(expr, v1, v2) for v2 in variables] for v1 in variables])
    return grad, hess


def find_critical_points(expr, var) -> list:
    """
    Find critical points by solving f'(x) = 0.
    Returns list of critical x values.
    """
    fp = diff(expr, var)
    return solve(fp, var)


def classify_critical_points(expr, var) -> list[dict]:
    """
    Find critical points and classify as min/max/inflection via second derivative test.
    """
    fpp  = diff(expr, var, 2)
    cpts = find_critical_points(expr, var)
    results = []
    for cp in cpts:
        fpp_val = fpp.subs(var, cp)
        kind    = "minimum" if fpp_val > 0 else "maximum" if fpp_val < 0 else "inflection"
        results.append({
            "x":         cp,
            "f(x)":      expr.subs(var, cp),
            "f''(x)":    fpp_val,
            "type":      kind,
        })
    return results


# ── 2. Equation solving ───────────────────────────────────────────────────────

def solve_system(
    equations: list,
    variables: list,
) -> list[dict] | dict:
    """
    Solve a system of algebraic equations.
    Returns list of solution dicts or a single solution dict.
    """
    solutions = solve(equations, variables, dict=True)
    return solutions


def solve_ode(
    ode_expr,
    func,
    ics: dict = None,
) -> Any:
    """
    Solve an ODE symbolically.
    ode_expr: the ODE (set equal to 0 or use Eq)
    func:     the function to solve for, e.g. f(x)
    ics:      initial conditions dict, e.g. {f(0): 1, f'(0): 0}
    
    Example:
        x = symbols("x")
        f = Function("f")
        sol = solve_ode(f(x).diff(x,2) + f(x), f(x), {f(0): 0, f(x).diff(x).subs(x,0): 1})
    """
    sol = dsolve(ode_expr, func, ics=ics)
    return sol


def solve_inequality(expr, var) -> Any:
    """Solve a polynomial inequality (returns SymPy set)."""
    return solve_univariate_inequality(expr, var, relational=False)


def characteristic_polynomial(matrix) -> Any:
    """Compute characteristic polynomial of a matrix."""
    lam = symbols("lambda")
    n   = matrix.shape[0]
    I_n = eye(n)
    return expand(det(matrix - lam * I_n))


# ── 3. Linear algebra ─────────────────────────────────────────────────────────

def symbolic_matrix_analysis(M: Matrix) -> dict:
    """
    Full symbolic analysis of a square matrix:
    determinant, trace, rank, eigenvalues, reduced row echelon form.
    """
    eigenval_dict  = M.eigenvals()
    eigenvect_list = M.eigenvects()

    return {
        "det":        M.det(),
        "trace":      M.trace(),
        "rank":       M.rank(),
        "rref":       M.rref()[0],
        "eigenvals":  eigenval_dict,   # {eigenvalue: multiplicity}
        "eigenvects": [
            {"eigenvalue": ev, "multiplicity": mul, "vectors": vecs}
            for ev, mul, vecs in eigenvect_list
        ],
        "is_symmetric":    M == M.T,
        "is_invertible":   M.det() != 0,
    }


def gram_schmidt(vectors: list[Matrix]) -> list[Matrix]:
    """Gram-Schmidt orthogonalization of a list of column vectors."""
    orthogonal = GramSchmidt(vectors, orthonormal=True)
    return orthogonal


def solve_linear_system(A: Matrix, b: Matrix) -> Matrix | None:
    """Solve Ax = b symbolically. Returns solution vector or None if no unique solution."""
    aug = A.row_join(b)
    sol = Matrix(linsolve(aug, *symbols(f"x0:{A.shape[1]}")))
    return sol if sol else None


# ── 4. Series and approximations ──────────────────────────────────────────────

def fourier_coefficients(
    expr,
    var,
    n_terms: int = 5,
    period:  float = 2 * float(pi),
) -> dict:
    """
    Compute Fourier series coefficients a_n and b_n for the first n_terms.
    Assumes period L (default 2π → standard Fourier series).
    """
    L   = Rational(period) / 2
    a0  = (1 / L) * integrate(expr, (var, -L, L))
    a_n = {}
    b_n = {}
    n   = symbols("n", integer=True, positive=True)

    for k in range(1, n_terms + 1):
        a_n[k] = simplify((1/L) * integrate(expr * cos(k*pi*var/L), (var, -L, L)))
        b_n[k] = simplify((1/L) * integrate(expr * sin(k*pi*var/L), (var, -L, L)))

    return {"a0": simplify(a0), "a_n": a_n, "b_n": b_n}


# ── 5. Symbolic ↔ Numeric bridge ──────────────────────────────────────────────

def to_numpy_function(expr, variables: list, modules: str = "numpy"):
    """
    Convert a SymPy expression to a fast NumPy-compatible callable.
    The returned function accepts NumPy arrays.
    """
    return lambdify(variables, expr, modules=modules)


def evaluate_at(expr, substitutions: dict) -> float:
    """Substitute numeric values and return float result."""
    return float(expr.subs(substitutions).evalf())


def high_precision(expr, digits: int = 50) -> str:
    """Evaluate to arbitrary precision string (e.g. for irrational constants)."""
    return str(N(expr, digits))


# ── 6. Probability with SymPy.stats ──────────────────────────────────────────

def normal_distribution_moments(mu_val, sigma_val) -> dict:
    """Compute exact mean, variance, and moment-generating function of Normal(μ, σ)."""
    x   = symbols("x")
    t   = symbols("t")
    mu  = symbols("mu", real=True)
    sig = symbols("sigma", positive=True)
    X   = Normal("X", mu_val, sigma_val)

    return {
        "mean":     E(X),
        "variance": variance(X),
        "std":      sqrt(variance(X)),
        "mgf":      E(exp(t * X)),
    }


# ── Demo ──────────────────────────────────────────────────────────────────────

if __name__ == "__main__":
    print("SymPy Symbolic Mathematics Demo")
    print("=" * 50)

    x, y = symbols("x y", real=True)
    f    = Function("f")

    # Calculus demo
    expr = sin(x) * exp(-x**2)
    result = diff_and_integrate(expr, x)
    print(f"\nd/dx[sin(x)·e^(-x²)] = {result['derivative']}")

    # Taylor series
    s = taylor_series(sin(x) / x, x, 0, 8)
    print(f"\nsin(x)/x Taylor: {s}")

    # Critical points
    g = x**4 - 4*x**2 + 1
    cpts = classify_critical_points(g, x)
    print(f"\nCritical points of x⁴-4x²+1:")
    for cp in cpts:
        print(f"  x={cp['x']}, f={cp['f(x)']}, type={cp['type']}")

    # ODE
    ode = f(x).diff(x, 2) + 4*f(x)
    sol = solve_ode(Eq(ode, 0), f(x), ics={f(0): 0, f(x).diff(x).subs(x, 0): 2})
    print(f"\nODE f''(x)+4f(x)=0, f(0)=0, f'(0)=2: {sol}")

    # Matrix eigenvalues
    M = Matrix([[2, 1], [1, 3]])
    analysis = symbolic_matrix_analysis(M)
    print(f"\nMatrix [[2,1],[1,3]] eigenvals: {analysis['eigenvals']}")

    # lambdify
    expr_np = x**2 + 2*sin(x)
    fn = to_numpy_function(expr_np, [x])
    xs = np.linspace(0, 2*float(pi), 5)
    print(f"\nlambdify(x²+2sin(x)) at [0,π/2,π,3π/2,2π]: {fn(xs).round(4)}")

    # High-precision
    print(f"\nπ to 30 digits: {high_precision(pi, 30)}")
    print(f"e to 30 digits: {high_precision(E, 30)}")

For the NumPy/SciPy alternative for numerical computation — NumPy/SciPy compute quickly with floating-point approximations while SymPy performs exact rational arithmetic (Rational(1,3) + Rational(1,6) = Rational(1,2), not 0.4999...), produces closed-form antiderivatives that SciPy’s quad cannot, and lambdify bridges both worlds by converting verified symbolic expressions into vectorized NumPy functions, making SymPy the correct tool when exact mathematical derivation matters (textbook solutions, code generation, symbolic regression). For the Mathematica alternative when a full CAS notebook environment is needed — Mathematica is proprietary while SymPy is open-source, embeds directly in Python ML pipelines, and the dsolve ODE solver with initial conditions, GramSchmidt, and sympy.stats probability algebra all integrate with scikit-learn and PyTorch workflows, enabling automatic derivation of gradient expressions or loss function Jacobians that then get compiled to fast NumPy kernels. The Claude Skills 360 bundle includes SymPy skill sets covering symbolic differentiation and integration, Taylor and Fourier series, critical point classification, ODE solving with initial conditions, symbolic matrix analysis, Gram-Schmidt orthogonalization, arbitrary-precision evaluation, and lambdify NumPy bridges. Start with the free tier to try symbolic mathematics code generation.