Tutorials

Solving Equations with SymPy

Numerical solvers like `scipy.optimize.brentq` find approximate roots; SymPy's `solve` finds exact ones. When the answer to "find x where f(x) = 0" is `sqrt(5)/2` or `(1 + sqrt(5))/2`, SymPy returns that exact form — no rounding, no loss of information. You'd use it to derive closed-form solutions, verify analytical results, or solve parametric equations where the unknowns are symbolic. SymPy handles polynomials, rational functions, transcendental equations (with trig, exp, log), and systems of multiple equations.

### Solving a Single Equation

`solve(expr, var)` finds all values of `var` that make `expr` equal to zero. Pass an `Eq` object if you need to express `f(x) = g(x)` directly.

from sympy import symbols, solve, sqrt, Rational

x = symbols('x')

# Polynomial: x² - 5x + 6 = 0
sols = solve(x**2 - 5*x + 6, x)
print("x² - 5x + 6 = 0:", sols)

# Irrational roots: x² - 2 = 0
sols2 = solve(x**2 - 2, x)
print("x² - 2 = 0:", sols2)

# Quadratic with no rational roots: x² + x - 1 = 0
sols3 = solve(x**2 + x - 1, x)
print("x² + x - 1 = 0:", sols3)
x² - 5x + 6 = 0: [2, 3]
x² - 2 = 0: [-sqrt(2), sqrt(2)]
x² + x - 1 = 0: [-1/2 + sqrt(5)/2, -sqrt(5)/2 - 1/2]
- `solve(expr, x)` treats `expr` as equal to zero and solves for `x` — pass `Eq(lhs, rhs)` to solve `lhs = rhs` explicitly.
- The roots of `x² - 2` are returned as `[-sqrt(2), sqrt(2)]` — exact irrational numbers, not decimals.
- `(−1 + sqrt(5))/2` is the golden ratio φ − 1; SymPy returns the exact symbolic form rather than 0.618.

### solveset: A More Modern Interface

`solveset` returns a SymPy `Set` object, which handles infinite solution sets (like `sin(x) = 0`) more gracefully than `solve`'s list.

from sympy import symbols, solveset, sin, cos, Eq, S, pi

x = symbols('x')

# Algebraic
print(solveset(x**2 - 4, x, domain=S.Reals))

# Transcendental: sin(x) = 1/2
sols = solveset(Eq(sin(x), 0), x, domain=S.Reals)
print("sin(x) = 0 over reals:", sols)

# Restrict to [0, 2π]
from sympy import Interval
restricted = solveset(Eq(sin(x), Rational(1, 2)), x, domain=Interval(0, 2*pi))
print("sin(x) = 1/2 on [0, 2π]:", restricted)
{-2, 2}
sin(x) = 0 over reals: Union(ImageSet(Lambda(_n, 2*_n*pi), Integers), ImageSet(Lambda(_n, 2*_n*pi + pi), Integers))
sin(x) = 1/2 on [0, 2π]: {pi/6, 5*pi/6}
- `domain=S.Reals` restricts solutions to real numbers; use `S.Complexes` to include complex roots.
- `solveset(sin(x) = 0, domain=S.Reals)` returns an `ImageSet` describing all integer multiples of π — the infinite solution set.
- Passing an `Interval` domain restricts to finitely many solutions within that range.

### Systems of Equations

`solve` handles systems by passing a list of equations and a list of unknowns. It returns a list of dictionaries (or a single dictionary) mapping each variable to its solution.

from sympy import symbols, solve, Eq, sqrt

x, y = symbols('x y')

# Linear system: x + y = 5, 2x - y = 1
sol_linear = solve([x + y - 5, 2*x - y - 1], [x, y])
print("Linear system:", sol_linear)

# Nonlinear: circle and line intersection
# x² + y² = 4, y = x - 1
sol_nonlinear = solve([x**2 + y**2 - 4, y - x + 1], [x, y])
print("Circle ∩ Line:", sol_nonlinear)
Linear system: {x: 2, y: 3}
Circle ∩ Line: [(1/2 - sqrt(7)/2, -sqrt(7)/2 - 1/2), (1/2 + sqrt(7)/2, -1/2 + sqrt(7)/2)]
- Pass equations as expressions equal to zero (`x + y - 5` means `x + y = 5`), or as `Eq(lhs, rhs)` objects.
- The nonlinear system returns two solutions — the two intersection points of the circle and line — as a list of `(x, y)` tuples.
- SymPy may not find solutions for transcendental systems; `nsolve` provides a numerical fallback.

### Equations with Parameters

When an equation contains symbolic parameters, `solve` returns the solution in terms of those parameters — a closed-form expression you can analyze or specialize.

from sympy import symbols, solve, sqrt, simplify

x, a, b, c = symbols('x a b c')

# Quadratic formula: ax² + bx + c = 0
sols = solve(a*x**2 + b*x + c, x)
print("General quadratic solutions:")
for s in sols:
    print(" ", s)

# Verify the discriminant
discriminant = (sols[0] - sols[1])**2
print("\nDiscriminant (expanded):", simplify(discriminant.expand()))
General quadratic solutions:
  (-b - sqrt(-4*a*c + b**2))/(2*a)
  (-b + sqrt(-4*a*c + b**2))/(2*a)

Discriminant (expanded): (-4*a*c + b**2)/a**2
- `solve(a*x**2 + b*x + c, x)` returns the quadratic formula exactly: `(-b ± sqrt(b² - 4ac)) / (2a)`.
- `(sols[0] - sols[1])**2` computes the square of the difference between the two roots — expanding and simplifying gives `(b² - 4ac) / a²`.
- This pattern of "derive the formula once, substitute parameters later" avoids repeated numerical solving.

Combine symbolic solving with [symbolic differentiation](/tutorials/symbolic-differentiation-with-sympy): differentiate a function, set the derivative to zero, and call `solve` to find exact critical points analytically.