Eigenvalue Analysis: Vibration Modes¶

This notebook applies eigenvalue analysis to a mechanical vibration problem. We derive the mass and stiffness matrices symbolically, then compute natural frequencies and mode shapes numerically with np_eig (LAPACK). The workflow demonstrates the symbolic-numeric bridge: exact physics setup, numeric eigensolve, and visual interpretation.

Key functions: eigenvalues, np_eig, np_inv, np_matmul, np_transpose.

In [127]:
load("numerics")$
load("ax-plots")$

Mass-Spring System [S+N]¶

Consider three masses connected in a line by four springs, anchored to walls at both ends:

wall ---[k1]--- m1 ---[k2]--- m2 ---[k3]--- m3 ---[k4]--- wall

The equations of motion for free vibration ($F_{\text{ext}} = 0$) are:

$$M \ddot{x} = -K x$$

where $M$ is the diagonal mass matrix and $K$ is the tridiagonal stiffness matrix. We derive $K$ symbolically from the spring forces.

In [128]:
/* Symbolic mass and stiffness matrices */
M_sym : matrix(
  [m1, 0, 0],
  [0, m2, 0],
  [0, 0, m3]
);

/* Stiffness: F_i = sum of spring forces on mass i */
/* Mass 1: -k1*x1 + k2*(x2 - x1) = -(k1+k2)*x1 + k2*x2           */
/* Mass 2: -k2*(x2-x1) + k3*(x3-x2) = k2*x1 - (k2+k3)*x2 + k3*x3 */
/* Mass 3: -k3*(x3-x2) - k4*x3 = k3*x2 - (k3+k4)*x3              */
K_sym : matrix(
  [k1+k2,  -k2,      0],
  [-k2,    k2+k3,   -k3],
  [0,      -k3,     k3+k4]
);

print("Mass matrix M:")$
M_sym;
print("Stiffness matrix K:")$
K_sym;

matrix([m1,0,0],[0,m2,0],[0,0,m3])
matrix([k2+k1,-k2,0],[-k2,k3+k2,-k3],[0,-k3,k4+k3])
Mass matrix M:
matrix([m1,0,0],[0,m2,0],[0,0,m3])
Stiffness matrix K:
Out[128]:
In [129]:
/* Assign numeric values */
params : [m1=1, m2=2, m3=1, k1=10, k2=10, k3=10, k4=10]$

M_num : float(subst(params, M_sym));
K_num : float(subst(params, K_sym));

print("M (numeric):")$
M_num;
print("K (numeric):")$
K_num;

matrix([1.0,0.0,0.0],[0.0,2.0,0.0],[0.0,0.0,1.0])
matrix([20.0,-10.0,0.0],[-10.0,20.0,-10.0],[0.0,-10.0,20.0])
M (numeric):
matrix([1.0,0.0,0.0],[0.0,2.0,0.0],[0.0,0.0,1.0])
K (numeric):
Out[129]:

Eigenvalue Problem [S+N]¶

For harmonic motion $x(t) = \phi\, e^{i\omega t}$, the equation $M\ddot{x} = -Kx$ becomes the generalised eigenvalue problem:

$$K \phi = \omega^2 M \phi$$

We convert to the standard form $M^{-1}K\,\phi = \omega^2\,\phi$ and solve with np_eig. The eigenvalues are $\lambda_i = \omega_i^2$ and the eigenvectors $\phi_i$ are the mode shapes.

In [130]:
/* Convert to ndarrays */
M_nd : ndarray(M_num)$
K_nd : ndarray(K_num)$

/* Standard eigenvalue problem: M^{-1} K */
MinvK : np_matmul(np_inv(M_nd), K_nd)$

print("M^{-1} K:")$
np_to_matrix(MinvK);

M^{-1} K:
Out[130]:
In [131]:
/* Eigendecomposition */
[eigenvals, eigenvecs] : np_eig(MinvK)$

lambdas : np_to_list(eigenvals)$
print("Eigenvalues (omega^2):")$
for i : 1 thru length(lambdas) do
  print(sconcat("  lambda_", i, " = ", lambdas[i]))$

/* Natural frequencies */
print("Natural frequencies (omega = sqrt(lambda)):")$
for i : 1 thru length(lambdas) do
  print(sconcat("  omega_", i, " = ", sqrt(abs(lambdas[i])),
                " rad/s  (f = ", sqrt(abs(lambdas[i]))/(2*float(%pi)), " Hz)"))$

Eigenvalues (omega^2):
  lambda_1 = 3.819660112501051
  lambda_2 = 20.000000000000004
  lambda_3 = 26.180339887498956
Natural frequencies (omega = sqrt(lambda)):
  omega_1 = 1.9543950758485478 rad/s  (f = 0.3110516370757561 Hz)
  omega_2 = 4.47213595499958 rad/s  (f = 0.7117625434171772 Hz)
  omega_3 = 5.116672736016928 rad/s  (f = 0.8143437581206265 Hz)

Mode Shapes¶

Each eigenvector describes a mode shape -- the relative displacement pattern of the masses at that natural frequency. We extract and plot each mode.

  • Mode 1 (lowest frequency): all masses move in the same direction
  • Mode 2: the centre mass stays still while outer masses move opposite
  • Mode 3 (highest frequency): alternating motion, maximum relative displacement
In [132]:
/* Extract eigenvectors (columns of the eigenvector matrix) */
V_mat : np_to_matrix(eigenvecs)$

print("Eigenvector matrix (columns = mode shapes):")$
V_mat;

Eigenvector matrix (columns = mode shapes):
Out[132]:
In [133]:
/* Plot mode shapes as bar charts */
positions : [1, 2, 3]$
mode_labels : ["Mode 1 (lowest freq)", "Mode 2 (middle freq)", "Mode 3 (highest freq)"]$

/* Extract each mode shape from the eigenvector matrix columns */
modes : makelist(
  makelist(V_mat[row][col], row, 1, 3),
  col, 1, 3
)$

/* Sort modes by eigenvalue magnitude (ascending frequency) */
indexed : makelist([abs(lambdas[i]), i], i, 1, 3)$
indexed : sort(indexed, lambda([a, b], a[1] < b[1]))$
sorted_indices : makelist(indexed[i][2], i, 1, 3)$

for k : 1 thru 3 do block(
  [idx : sorted_indices[k], mode_data],
  mode_data : makelist([positions[j], modes[idx][j]], j, 1, 3),
  print(sconcat(mode_labels[k],
                "  (omega = ", sqrt(abs(lambdas[idx])), " rad/s)")),
  print("  Displacements:", modes[idx])
)$

Mode 1 (lowest freq)  (omega = 1.9543950758485478 rad/s)
  Displacements: [-0.4653411271949863,-0.752937760164676,-0.46534112719498616]
Mode 2 (middle freq)  (omega = 4.47213595499958 rad/s)
  Displacements:
                [0.7071067811865469,3.6721064886523107e-16,
                 -0.7071067811865481]
Mode 3 (highest freq)  (omega = 5.116672736016928 rad/s)
  Displacements: [0.6479361632942993,-0.4004465714560786,0.6479361632942977]
In [134]:
/* Visualise all three mode shapes */
colors : ["#377eb8", "#e41a1c", "#4daf4a"]$

/* Build displacement data for plotting: include wall anchors at positions 0 and 4 */
mode_plots : []$
for k : 1 thru 3 do block(
  [idx : sorted_indices[k],
   disp, pts],
  disp : modes[idx],
  /* Normalize so max displacement = 1 */
  disp : disp / lmax(map(abs, disp)),
  pts : [[0, 0], [1, disp[1]], [2, disp[2]], [3, disp[3]], [4, 0]]
)$

/* Plot mode shapes with walls anchored at zero */
ax_draw2d(
  color=colors[1], name=sconcat("Mode 1 (w=", float(sqrt(abs(lambdas[sorted_indices[1]]))), ")"),
  lines(
    [0, 1, 2, 3, 4],
    block([d : modes[sorted_indices[1]]],
          d : d / lmax(map(abs, d)),
          append([0], d, [0]))
  ),
  color=colors[2], name=sconcat("Mode 2 (w=", float(sqrt(abs(lambdas[sorted_indices[2]]))), ")"),
  lines(
    [0, 1, 2, 3, 4],
    block([d : modes[sorted_indices[2]]],
          d : d / lmax(map(abs, d)),
          append([0], d, [0]))
  ),
  color=colors[3], name=sconcat("Mode 3 (w=", float(sqrt(abs(lambdas[sorted_indices[3]]))), ")"),
  lines(
    [0, 1, 2, 3, 4],
    block([d : modes[sorted_indices[3]]],
          d : d / lmax(map(abs, d)),
          append([0], d, [0]))
  ),
  title="Vibration Mode Shapes (normalised)",
  xlabel="Mass position", ylabel="Displacement",
  grid=true
)$
No description has been provided for this image

Symbolic Eigenvalues: 2-Mass System [S+N]¶

For a simpler 2-mass, 3-spring system we can compute eigenvalues symbolically to obtain closed-form natural frequencies.

wall ---[k]--- m ---[k]--- m ---[k]--- wall

With equal masses $m$ and equal springs $k$, the system is symmetric and the eigenvalues have a clean closed form.

In [135]:
/* 2-mass system: symbolic */
M2 : matrix(
  [m, 0],
  [0, m]
);

K2 : matrix(
  [2*k, -k],
  [-k,  2*k]
);

print("M:")$
M2;
print("K:")$
K2;

matrix([m,0],[0,m])
matrix([2*k,-k],[-k,2*k])
M:
matrix([m,0],[0,m])
K:
Out[135]:
In [136]:
/* Compute M^{-1} K symbolically */
MinvK2 : invert(M2) . K2;
MinvK2 : ratsimp(MinvK2);

print("M^{-1} K =")$
MinvK2;

matrix([(2*k)/m,-(k/m)],[-(k/m),(2*k)/m])
matrix([(2*k)/m,-(k/m)],[-(k/m),(2*k)/m])
M^{-1} K =
Out[136]:
In [137]:
/* Symbolic eigenvalues */
evals2 : eigenvalues(MinvK2);

print("Eigenvalues (omega^2):")$
print(evals2[1])$

print("Natural frequencies:")$
for ev in evals2[1] do
  print(sconcat("  omega = sqrt(", ev, ") = ", sqrt(ev)))$

[[k/m,(3*k)/m],[1,1]]
Eigenvalues (omega^2):
[k/m,(3*k)/m]
Natural frequencies:
  omega = sqrt(k/m) = sqrt(k/m)
  omega = sqrt((3*k)/m) = sqrt(3)*sqrt(k/m)
In [138]:
/* Verify numerically: k=10, m=1 */
M2_num : float(subst([m=1, k=10], M2))$
K2_num : float(subst([m=1, k=10], K2))$

M2_nd : ndarray(M2_num)$
K2_nd : ndarray(K2_num)$

[evals2_nd, evecs2_nd] : np_eig(np_matmul(np_inv(M2_nd), K2_nd))$

print("Numeric eigenvalues:", np_to_list(evals2_nd))$
print("Expected: k/m = 10, 3k/m = 30")$

/* Symbolic evaluation for comparison */
print("Symbolic at k=10, m=1:",
      float(subst([m=1, k=10], evals2[1])))$

Numeric eigenvalues: [30.0,10.000000000000002]
Expected: k/m = 10, 3k/m = 30
Symbolic at k=10, m=1: [10.0,30.0]

The symbolic eigenvalues are $\omega_1^2 = k/m$ and $\omega_2^2 = 3k/m$. For $k=10$, $m=1$: $\omega_1 = \sqrt{10} \approx 3.16$ and $\omega_2 = \sqrt{30} \approx 5.48$ rad/s. The numeric result confirms this exactly.

Verifying Orthogonality¶

For a symmetric positive-definite system, the eigenvectors should be $M$-orthogonal:

$$\phi_i^T M \phi_j = 0 \quad \text{for } i \neq j$$

This property is fundamental to modal analysis: it means the modes are independent and can be excited separately. We verify this numerically for the 3-mass system.

In [139]:
/* M-orthogonality check: V^T M V should be diagonal */
VtMV : np_matmul(np_transpose(eigenvecs), np_matmul(M_nd, eigenvecs))$

print("V^T M V (should be diagonal):")$
np_to_matrix(VtMV);

V^T M V (should be diagonal):
Out[139]:
In [140]:
/* Check individual inner products */
V_cols : makelist(np_reshape(np_col(eigenvecs, j), [3, 1]), j, 0, 2)$

print("M-orthogonality check (phi_i^T M phi_j):")$
for i : 1 thru 3 do
  for j : i thru 3 do block(
    [val],
    val : np_to_list(
      np_matmul(np_transpose(V_cols[i]),
                np_matmul(M_nd, V_cols[j]))
    ),
    if i = j then
      print(sconcat("  <phi_", i, ", M phi_", j, "> = ", val[1], "  (modal mass)"))
    else
      print(sconcat("  <phi_", i, ", M phi_", j, "> = ", val[1],
                    if abs(val[1]) < 1e-10 then "  (orthogonal)" else "  (NOT orthogonal!)"))
  )$

M-orthogonality check (phi_i^T M phi_j):
  <phi_1, M phi_1> = 1.5669152706817988  (modal mass)
  <phi_1, M phi_2> = -1.2046279426456211e-16  (orthogonal)
  <phi_1, M phi_3> = 3.3306690738754696e-16  (orthogonal)
  <phi_2, M phi_2> = 1.0  (modal mass)
  <phi_2, M phi_3> = -7.343864078089007e-18  (orthogonal)
  <phi_3, M phi_3> = 1.160357456590928  (modal mass)

The off-diagonal entries are zero (to machine precision), confirming $M$-orthogonality. The diagonal entries are the modal masses -- the effective mass participating in each mode.

Summary: This notebook showed the full eigenvalue analysis pipeline:

  1. Derive mass and stiffness matrices symbolically from physics
  2. Assign numeric parameters and convert to ndarrays
  3. Solve the eigenvalue problem with np_eig (LAPACK)
  4. Interpret eigenvalues as natural frequencies, eigenvectors as mode shapes
  5. Verify with symbolic eigenvalues for simple cases
  6. Validate orthogonality numerically

The symbolic-numeric bridge lets us keep the physics exact while leveraging LAPACK's robust eigensolvers for the heavy numerical work.