Symbolic Linear Algebra
While NumPy is the king of numerical linear algebra, SymPy allows us to perform linear algebra with variables. This is essential for deriving formulas, checking proofs, and solving systems where some parameters are not yet known.
Creating Symbolic Matrices
In SymPy, matrices are created using the Matrix class. Unlike NumPy arrays, these matrices are mutable by default and can contain any SymPy expression.
Interactive Lab
import sympy as sp
x, y = sp.symbols('x y')
# Create a 2x2 matrix
M = sp.Matrix([[x, y], [y, x]])
print("Symbolic Matrix M:\n")
sp.pprint(M)
# Basic operations
print("\nMatrix Squared:\n")
sp.pprint(M**2)
Expected output
Symbolic Matrix M:
⎡x y⎤
⎢ ⎥
⎣y x⎦
Matrix Squared:
⎡ 2 2 ⎤
⎢x + y 2⋅x⋅y ⎥
⎢ ⎥
⎣ 2⋅x⋅y 2 2 ⎦
x + y
1
2import sympy as sp
3x, y = sp.symbols('x y')
4
5
6M = sp.Matrix([[x, y], [y, x]])
7print("Symbolic Matrix M:\n")
8sp.pprint(M)
9
10
11print("\nMatrix Squared:\n")
12sp.pprint(M**2)
13
Determinants and Inversion
Calculating a determinant symbolically is a common task in stability analysis and geometry.
Interactive Lab
import sympy as sp
a, b, c, d = sp.symbols('a b c d')
# The general 2x2 matrix
A = sp.Matrix([[a, b], [c, d]])
det_A = A.det()
print(f"Determinant of A: {det_A}")
inv_A = A.inv()
print("\nInverse of A:")
sp.pprint(inv_A)
Expected output
Determinant of A: a*d - b*c
Inverse of A:
⎡ d -b ⎤
⎢─────── ───────⎥
⎢ad - bc ad - bc⎥
⎢ ⎥
⎢ -c a ⎥
⎢─────── ───────⎥
⎣ad - bc ad - bc⎦
1
2import sympy as sp
3a, b, c, d = sp.symbols('a b c d')
4
5
6A = sp.Matrix([[a, b], [c, d]])
7
8det_A = A.det()
9print(f"Determinant of A: {det_A}")
10
11inv_A = A.inv()
12print("\nInverse of A:")
13sp.pprint(inv_A)
14
Eigenvalues and Eigenvectors
One of the most powerful features of SymPy is calculating exact eigenvalues. Numerical methods might lose precision for nearly singular matrices or complex clusters; SymPy finds the exact roots of the characteristic polynomial.
Interactive Lab
import sympy as sp
lam = sp.symbols('lambda')
M = sp.Matrix([[3, -2], [4, -1]])
# Characteristic polynomial
poly = M.charpoly(lam)
print(f"Characteristic Polynomial: {poly.as_expr()}")
# Eigenvalues
eigenvals = M.eigenvals()
print(f"Eigenvalues: {eigenvals}") # Returns {value: multiplicity}
# Eigenvectors
eigenvects = M.eigenvects()
print("\nEigenvectors (Value, Multiplicity, Basis):")
for v in eigenvects:
print(v)
Expected output
Characteristic Polynomial: lambda**2 - 2*lambda + 5
Eigenvalues: {1 - 2*I: 1, 1 + 2*I: 1}
Eigenvectors (Value, Multiplicity, Basis):
(1 - 2*I, 1, [Matrix([
[1/2 + I/2],
[ 1]])])
(1 + 2*I, 1, [Matrix([
[1/2 - I/2],
[ 1]])])
1
2import sympy as sp
3lam = sp.symbols('lambda')
4M = sp.Matrix([[3, -2], [4, -1]])
5
6
7poly = M.charpoly(lam)
8print(f"Characteristic Polynomial: {poly.as_expr()}")
9
10
11eigenvals = M.eigenvals()
12print(f"Eigenvalues: {eigenvals}")
13
14
15eigenvects = M.eigenvects()
16print("\nEigenvectors (Value, Multiplicity, Basis):")
17for v in eigenvects:
18 print(v)
19
Matrix Decompositions
SymPy supports several decompositions, including LU, QR, and Diagonalization.
Interactive Lab
import sympy as sp
M = sp.Matrix([[1, 2], [2, 1]])
# Diagonalization: M = P * D * P^-1
P, D = M.diagonalize()
print("Diagonal Matrix D:")
sp.pprint(D)
print("\nTransformation Matrix P:")
sp.pprint(P)
# Verify
assert P * D * P.inv() == M
print("\nVerification successful: P * D * P^-1 == M")
Expected output
Diagonal Matrix D:
⎡-1 0⎤
⎢ ⎥
⎣ 0 3⎦
Transformation Matrix P:
⎡-1 1⎤
⎢ ⎥
⎣ 1 1⎦
Verification successful: P * D * P^-1 == M
1
2import sympy as sp
3M = sp.Matrix([[1, 2], [2, 1]])
4
5
6P, D = M.diagonalize()
7
8print("Diagonal Matrix D:")
9sp.pprint(D)
10print("\nTransformation Matrix P:")
11sp.pprint(P)
12
13
14assert P * D * P.inv() == M
15print("\nVerification successful: P * D * P^-1 == M")
16
Solving Linear Systems ()
You can solve a system of linear equations by passing a matrix and a column vector to the LUsolve method or by using solve_linear_system.
Interactive Lab
import sympy as sp
x1, x2 = sp.symbols('x1 x2')
# System:
# x1 + 2*x2 = 5
# 3*x1 + 4*x2 = 11
A = sp.Matrix([[1, 2], [3, 4]])
b = sp.Matrix([5, 11])
sol = A.LUsolve(b)
print(f"Solution vector: {sol}")
Expected output
Solution vector: Matrix([[1], [2]])
1
2import sympy as sp
3x1, x2 = sp.symbols('x1 x2')
4
5
6
7
8
9A = sp.Matrix([[1, 2], [3, 4]])
10b = sp.Matrix([5, 11])
11
12sol = A.LUsolve(b)
13print(f"Solution vector: {sol}")
14
Applications in Engineering
Symbolic matrices are often used to define Stiffness Matrices in structural analysis or Jacobian Matrices in robot kinematics. Because the variables are preserved, we can calculate the Jacobian once and then substitute specific joint angles thousands of times during a simulation.
Interactive Lab
import sympy as sp
theta = sp.symbols('theta')
# Rotation matrix around Z-axis
R = sp.Matrix([
[sp.cos(theta), -sp.sin(theta), 0],
[sp.sin(theta), sp.cos(theta), 0],
[0, 0, 1]
])
# Derivative of rotation matrix with respect to theta
dR_dtheta = R.diff(theta)
print("Rotation Matrix Derivative:")
sp.pprint(dR_dtheta)
Expected output
Rotation Matrix Derivative:
⎡-sin(theta) -cos(theta) 0⎤
⎢ ⎥
⎢ cos(theta) -sin(theta) 0⎥
⎢ ⎥
⎣ 0 0 0⎦
1
2import sympy as sp
3theta = sp.symbols('theta')
4
5
6R = sp.Matrix([
7 [sp.cos(theta), -sp.sin(theta), 0],
8 [sp.sin(theta), sp.cos(theta), 0],
9 [0, 0, 1]
10])
11
12
13dR_dtheta = R.diff(theta)
14
15print("Rotation Matrix Derivative:")
16sp.pprint(dR_dtheta)
17
In the next module, we will conclude our deep dive into the scientific stack by looking at how to bridge these symbolic results back into the numerical world of SciPy.