Tutorial 7

Interactive Notebooks

Click the badges above to run this tutorial interactively in your browser without installing anything!

• Binder: Free cloud-based Jupyter environment
• Colab: Google's free Jupyter notebook environment

Topics Covered

This tutorial covers the following topics using Python and SymPy:

• Eigenvalue / eigenvector problems for a 3-DOF mass-spring system
• Characteristic polynomial and the natural frequencies $\omega=\sqrt{\lambda}$
• Normalised and unscaled eigenvectors
• Diagonalisation: the eigenvector matrix $P$ and the spectral matrix $D$
• Powers of a matrix via $A^{k}=PD^{k}P^{-1}$
• Cayley-Hamilton theorem

Introduction

import sympy as sy
from sympy.abc import x, y, z, t
from sympy import Matrix, Rational, eye, zeros, diag, sqrt, pi, N, simplify, solve
sy.init_printing(use_latex=True)

lam = sy.Symbol('lambda')        # the eigenvalue symbol used in the characteristic polynomial

Question 1

The 3-DOF mass-spring system has coefficient matrix

$$A=\begin{bmatrix}\frac{k_{1}+k_{2}}{m_{1}} & -\frac{k_{2}}{m_{1}} & 0\\[4pt] -\frac{k_{2}}{m_{2}} & \frac{k_{2}+k_{3}}{m_{2}} & -\frac{k_{3}}{m_{2}}\\[4pt] 0 & -\frac{k_{3}}{m_{2}} & \frac{k_{3}+k_{4}}{m_{3}}\end{bmatrix} =\begin{bmatrix}\frac{100}{3} & -\frac{70}{3} & 0\\[4pt] -\frac{70}{3} & \frac{140}{3} & -\frac{70}{3}\\[4pt] 0 & -\frac{70}{3} & \frac{100}{3}\end{bmatrix}$$

with $k_{1}=k_{4}=15$, $k_{2}=k_{3}=35\ \text{N/m}$ and $m_{1}=m_{2}=m_{3}=1.5\ \text{kg}$.

A = Matrix([[Rational(100, 3), Rational(-70, 3), 0],
            [Rational(-70, 3), Rational(140, 3), Rational(-70, 3)],
            [0, Rational(-70, 3), Rational(100, 3)]])

display(A.charpoly(lam).as_expr())            # characteristic polynomial
display(solve(A.charpoly(lam).as_expr(), lam))    # exact eigenvalues

The characteristic equation is $-\lambda^{3}+\frac{340}{3}\lambda^{2}-\frac{9400}{3}\lambda+\frac{140000}{9}=0$, whose roots are $\lambda=40\mp\frac{10\sqrt{102}}{3}$ and $\lambda=\frac{100}{3}$.

ev = sorted(A.eigenvects(), key=lambda item: N(item[0]))   # ascending eigenvalue
display([N(item[0], 8) for item in ev])                    # [6.3349835, 33.333333, 73.665016]

lam1, _, vecs1 = ev[0]                                     # smallest eigenvalue
v1 = vecs1[0]
display(lam1, N(lam1, 8))                                  # 40 - 10*sqrt(102)/3 = 6.3350
display(v1.T, [N(c, 6) for c in v1])                       # [1, 1.1571, 1]

The smallest eigenvalue is $\lambda_{1}=\omega_{1}^{2}=6.3350$, with unscaled eigenvector $\{1,\ 1.1571,\ 1\}^{\mathsf{T}}$.

note

The printed solution quotes $\lambda_{1}=3839/606=6.3350$ and $\lambda_{3}=14954/203=73.6650$. Those are decimal approximations from a numerical root finder — the exact values are $\lambda_{1}=40-\frac{10\sqrt{102}}{3}$ and $\lambda_{3}=40+\frac{10\sqrt{102}}{3}$. SymPy's solve returns the exact radicals.

Question 2

The second largest eigenvalue and its normalised eigenvector.

lam2, _, vecs2 = ev[1]                     # middle eigenvalue
v2 = vecs2[0]
display(lam2, N(lam2, 8))                  # 100/3 = 33.3333
display(v2.T)                              # [-1, 0, 1]

display(v2.normalized().T)                 # [-sqrt(2)/2, 0, sqrt(2)/2]
display([N(c, 6) for c in v2.normalized()])    # [-0.7071, 0, 0.7071]

So $\lambda_{2}=\omega_{2}^{2}=\frac{100}{3}=33.3333$ with normalised eigenvector $\frac{1}{\sqrt{2}}\{-1,\ 0,\ 1\}^{\mathsf{T}}$.

Question 3

The largest eigenvalue and its unscaled eigenvector.

lam3, _, vecs3 = ev[2]                     # largest eigenvalue
v3 = vecs3[0]
display(lam3, N(lam3, 8))                  # 40 + 10*sqrt(102)/3 = 73.6650
display(v3.T, [N(c, 6) for c in v3])       # [1, -1.7285, 1]

So $\lambda_{3}=\omega_{3}^{2}=73.6650$ with unscaled eigenvector $\{1,\ -1.7285,\ 1\}^{\mathsf{T}}$.

Question 4

Combine the results to build the eigenvector matrix $P$ and diagonalise $A$.

vecs = [item[2][0] for item in ev]
vecs[1] = vecs[1].normalized()             # use the normalised middle eigenvector
P = Matrix.hstack(*vecs)                   # eigenvector (modal) matrix
D = diag(*[item[0] for item in ev])        # spectral matrix

display(N(P, 6))
display(N(D, 6))

# verify  P^-1 A P = D
display(simplify(P.inv()*A*P - D) == zeros(3, 3))     # True

$$P=\begin{bmatrix}1 & -\frac{1}{\sqrt{2}} & 1\\[4pt] 1.1571 & 0 & -1.7285\\[4pt] 1 & \frac{1}{\sqrt{2}} & 1\end{bmatrix}, \qquad D=\begin{bmatrix}6.335 & 0 & 0\\0 & 33.3333 & 0\\0 & 0 & 73.665\end{bmatrix}$$

The diagonal matrix $D$ is the eigenvalue matrix: its diagonal entries are exactly the eigenvalues of $A$, in the same order as the columns of $P$.

Question 5

Compute $A^{50}$ and comment on the eigenvalues and eigenvectors.

A50 = P * (D**50) * P.inv()
display(N(A50 / 10**93, 6))

$$A^{50}=10^{93}\begin{bmatrix} 0.4626 & -0.7995 & 0.4626\\ -0.7995 & 1.3820 & -0.7995\\ 0.4626 & -0.7995 & 0.4626\end{bmatrix}$$

The eigenvalues of $A^{50}$ are $\lambda_{i}^{50}$ (each eigenvalue is raised to the same power), while the eigenvectors of $A^{50}$ are unchanged from those of $A$.

Question 6

$B=\begin{bmatrix}1&2&3\\0&1&1\\0&0&2\end{bmatrix}$ with $\det(B)=2$ has an eigenvalue of $2$. Find the remaining eigenvalues and the characteristic equation.

B = Matrix([[1, 2, 3], [0, 1, 1], [0, 0, 2]])

display(B.trace(), B.det())                 # 4, 2
display(B.eigenvals())                      # {1: 2, 2: 1} -> 1, 1, 2
display(B.charpoly(lam).as_expr())          # lambda**3 - 4*lambda**2 + 5*lambda - 2

Using $\operatorname{trace}(B)=\sum\lambda_{i}=4$ and $\det(B)=\prod\lambda_{i}=2$ with $\lambda_{3}=2$ gives $\lambda_{1}+\lambda_{2}=2$ and $\lambda_{1}\lambda_{2}=1$, so $\lambda_{1}=\lambda_{2}=1$. The characteristic equation is $\lambda^{3}-4\lambda^{2}+5\lambda-2=0$ — obtained from the trace, the sum of the principal $2\times2$ minors and the determinant, without ever forming $\det(B-\lambda I)$.

Question 7

Verify the Cayley-Hamilton statements and compute $B^{5}$.

display(simplify(Rational(1, 2)*B**2 - 2*B + Rational(5, 2)*eye(3) - B.inv()))   # zero matrix

display(B**5)
display(26*B**2 - 47*B + 22*eye(3))         # same matrix -> formula verified

display(B**6)
display(57*B**2 - 108*B + 52*eye(3))       # same matrix -> formula verified

$$\mathbf{B}^{-1}=\tfrac{1}{2}\mathbf{B}^{2}-2\mathbf{B}+\tfrac{5}{2}\mathbf{I}, \qquad \mathbf{B}^{6}=57\mathbf{B}^{2}-108\mathbf{B}+52\mathbf{I}, \qquad \mathbf{B}^{5}=\begin{bmatrix}1&10&145\\0&1&31\\0&0&32\end{bmatrix}$$

Question 8

Why does $B^{5}=PD^{5}P^{-1}$ fail here?

P8 = Matrix([[1, 1, 5],
             [0, 0, 1],
             [0, 0, 1]])
display(P8.det())                           # 0
display(P8.rank())                          # 2 -> only 2 independent eigenvectors

The repeated eigenvalue $\lambda=1$ has algebraic multiplicity $2$ but only one independent eigenvector, so the eigenvector matrix has two identical columns and $\det(P)=0$. Since $P$ is singular, $P^{-1}$ does not exist and the diagonalisation formula cannot be used. The Cayley-Hamilton theorem has no such restriction — which is exactly why it was used in Q7.

Question 9

Advantage and disadvantage of each method.

	Diagonalisation $B^{k}=PD^{k}P^{-1}$	Cayley-Hamilton
Advantage	Very fast once $P$ and $D$ are known — a single formula gives any power.	Derived from the characteristic equation alone; no eigenvectors needed.
Disadvantage	Needs the complete eigenvector set, and fails when $P$ is not invertible (repeated roots, defective matrix).	The recurrence for each new power must be re-derived by hand.

Question 10

$C=\begin{bmatrix}0&1&1\\1&0&1\\1&1&0\end{bmatrix}$ — find all eigenvalues and normalised eigenvectors, then verify $CP=PD$.

C = Matrix([[0, 1, 1], [1, 0, 1], [1, 1, 0]])

display(C.charpoly(lam).as_expr())          # lambda**3 - 3*lambda - 2
display(C.eigenvals())                      # {-1: 2, 2: 1}

for val, mult, vecs in C.eigenvects():
    display(val, mult, [[N(c, 6) for c in v] for v in vecs],
            [list(v.normalized()) for v in vecs])

The eigenvalues are $\lambda_{1}=\lambda_{2}=-1$ (multiplicity 2) and $\lambda_{3}=2$, with normalised eigenvectors

$$P=\begin{bmatrix}-\frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{3}}\\[4pt] \frac{1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{3}}\\[4pt] 0 & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{3}}\end{bmatrix}, \qquad D=\begin{bmatrix}-1&0&0\\0&-1&0\\0&0&2\end{bmatrix}$$

P = Matrix([[-1/sqrt(2), -1/sqrt(2), 1/sqrt(3)],
            [ 1/sqrt(2),          0, 1/sqrt(3)],
            [         0,  1/sqrt(2), 1/sqrt(3)]])
D = diag(-1, -1, 2)

display(N(C*P, 6))                  # LHS
display(N(P*D, 6))                  # RHS
display(simplify(C*P - P*D) == zeros(3, 3))     # True -> verified

Since $CP=PD$, the eigenvalue matrix $D$ and eigenvector matrix $P$ satisfy the eigenvalue/eigenvector problem $(C-\lambda I)x=0$.

---

Extra Learning Resources

Key Concepts to Master

1. Characteristic equation: $\det(A-\lambda I)=0$
2. Trace and determinant checks: $\sum\lambda_{i}=\operatorname{trace}(A)$ and $\prod\lambda_{i}=\det(A)$
3. Normalisation: divide an eigenvector by its magnitude so that $\|v\|=1$
4. Diagonalisation: $D=P^{-1}AP$ exists only when $P$ (the eigenvector matrix) is invertible
5. Cayley-Hamilton: a matrix satisfies its own characteristic equation, so $A^{k}$ can always be reduced to a combination of $I,A,\ldots,A^{n-1}$

SymPy Resources

• SymPy charpoly — the characteristic polynomial
• SymPy eigenvects — eigenvalues with multiplicities and eigenvectors
• SymPy diagonalize — returns (P, D) directly

Common Pitfalls

• eigenvects() returns a list of tuples (eigenvalue, multiplicity, [vectors]) — unpack it before use
• A repeated eigenvalue does not guarantee repeated independent eigenvectors; when it does not, the matrix is defective and cannot be diagonalised
• Eigenvectors are only defined up to a scalar multiple — scaled and unscaled answers are both correct
• N(expr, n) gives n significant figures, not decimal places