NumPy and SciPy work with numerical matrices. SymPy's `Matrix` class works with symbolic matrices — matrices whose entries can be variables, fractions, or algebraic expressions. This lets you derive formulas for inverses, eigenvalues, or solutions to linear systems that hold for all valid parameter values, not just a specific set of numbers. You'd use it to derive the general solution to a parametric linear system, compute eigenvalues of a symbolic Jacobian, or verify matrix identities analytically. ### Creating and Inspecting Matrices `sympy.Matrix` accepts a nested list of values. Entries can be integers, fractions, or symbolic expressions.
from sympy import Matrix, symbols, sqrt, Rational
x, a = symbols('x a')
A = Matrix([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
B = Matrix([[a, 1],
[0, a**2]])
print("A =")
print(A)
print("\nShape:", A.shape)
print("Trace:", A.trace())
print("Rank:", A.rank())
print("\nSymbolic B =")
print(B)
print("det(B) =", B.det())- `Matrix([[...]])` stores entries exactly — fractions stay as fractions, symbols stay symbolic. - `A.trace()` sums the diagonal; `A.rank()` computes the rank using exact row reduction (no floating-point rounding). - `B.det()` returns `a³ − 0 = a³` — a symbolic expression showing how the determinant depends on `a`. ### Determinant and Inverse `det()` and `inv()` compute exact results. For symbolic matrices, the inverse is a matrix of rational functions of the symbolic parameters.
from sympy import Matrix, symbols
a, b = symbols('a b')
A = Matrix([[2, 1],
[5, 3]])
print("A:")
print(A)
print("\ndet(A) =", A.det())
print("\nA_inv =")
print(A.inv())
print("\nVerify A @ A_inv = I:")
print(A * A.inv())
# Symbolic matrix
B = Matrix([[a, b],
[1, a]])
print("\ndet(B) =", B.det())
print("\nB_inv =")
print(B.inv())
- `A.inv()` uses the adjugate formula for 2×2 matrices and exact Gaussian elimination for larger ones — no numerical error.
- `A * A.inv()` should return the identity matrix with exact integers, confirming correctness.
- `B.inv()` returns entries like `a/(a²−b)` — expressions that reveal when the inverse is undefined (when `a² = b`).
### Eigenvalues and Eigenvectors
`eigenvals()` returns a dictionary of `{eigenvalue: multiplicity}`. `eigenvects()` returns a list of `(eigenvalue, multiplicity, [eigenvectors])` tuples.
from sympy import Matrix, sqrt, symbols
A = Matrix([[3, 1],
[1, 3]])
vals = A.eigenvals()
vects = A.eigenvects()
print("Eigenvalues (value: multiplicity):", vals)
print()
for val, mult, vecs in vects:
print(f"λ = {val}, multiplicity = {mult}")
print(f" eigenvector: {vecs[0].T}")
# Verify: A @ v = λ * v
for val, mult, vecs in vects:
v = vecs[0]
print(f"A·v = {(A * v).T} | λ·v = {(val * v).T}")- `eigenvals()` returns eigenvalues as exact symbolic values — `2` and `4` for this symmetric matrix. - `eigenvects()` returns the full eigenvectors as column matrices (not normalized unless you call `.normalized()`). - Verifying `A @ v = λ * v` symbolically confirms correctness — the result should be the zero vector when subtracting. ### Solving a Symbolic Linear System `solve_linear_system` solves Ax = b symbolically, returning solutions in terms of any free variables or symbolic parameters.
from sympy import Matrix, symbols, linsolve
x, y, z = symbols('x y z')
# 3×3 system: unique solution
A = Matrix([[2, 1, -1],
[-3, -1, 2],
[-2, 1, 2]])
b = Matrix([8, -11, -3])
# Augmented matrix approach
aug = A.row_join(b)
sol = A.solve(b)
print("Solution:", sol.T)
print("Verify Ax - b:", (A * sol - b).T)
# Symbolic RHS
p, q = symbols('p q')
A2 = Matrix([[1, 2], [3, 4]])
b2 = Matrix([p, q])
sol2 = A2.solve(b2)
print("\nSymbolic solution:")
print("x =", sol2[0])
print("y =", sol2[1])- `A.solve(b)` uses exact Gaussian elimination — the result is an exact rational or symbolic expression, not a floating-point approximation. - For the symbolic RHS, the solution expresses x and y as linear functions of p and q — this is the general solution valid for any p, q. - `A.solve` raises `ValueError` if A is singular. Use `A.nullspace()` to analyze the homogeneous system when solutions exist but are not unique. ### Matrix Functions and Decompositions SymPy can compute diagonalization, the matrix exponential, and other matrix functions exactly.
from sympy import Matrix, exp, symbols, eye
A = Matrix([[1, 1],
[0, 2]])
# Diagonalization: A = P @ D @ P_inv
P, D = A.diagonalize()
print("P =")
print(P)
print("\nD (diagonal eigenvalue matrix) =")
print(D)
print("\nVerify P @ D @ P_inv == A:", P * D * P.inv() == A)
# Matrix exponential via diagonalization: expm(A) = P @ expm(D) @ P_inv
from sympy import exp as sym_exp, zeros
expD = zeros(2)
for i in range(2):
expD[i, i] = sym_exp(D[i, i])
expmA = P * expD * P.inv()
print("\nexpm(A) =")
print(expmA)- `A.diagonalize()` returns `(P, D)` where D is diagonal and `A = P @ D @ P⁻¹`. It raises an error if A is not diagonalizable. - The matrix exponential via diagonalization works by exponentiating each diagonal entry of D individually — exact for diagonalizable matrices. - `P * D * P.inv() == A` returns `True` — SymPy verifies matrix equality symbolically. Symbolic matrix algebra is particularly useful when combined with [solving equations with SymPy](/tutorials/solving-equations-with-sympy) — set the characteristic polynomial `det(A - λI) = 0` and solve for the eigenvalues analytically. For large numerical linear systems, see [linear algebra with scipy.linalg](/tutorials/linear-algebra-with-scipy-linalg).