Stand up a finite-element solver for natural frequencies of a slender beam. Verify against the closed-form βL roots, then perturb until the closed form is gone.
An Euler–Bernoulli beam, length L, modulus E, area moment I, density ρ, cross-section A. Free transverse vibration is governed by
EI · ∂⁴w/∂x⁴ + ρA · ∂²w/∂t² = 0
Pick boundary conditions from {cantilever, simply-supported, clamped–clamped, free–free}. You want: the first few natural frequencies ωn and their mode shapes.
Finite-element discretization with Hermite cubic elements turns the PDE into a generalized matrix eigenproblem K φ = ω² M φ. Your brain picks the mesh, applies the BCs, judges whether the answer is believable. The numerics library does the eigensolve. The closed-form βL roots give you a check you can run by hand.
Start with a uniform beam, 1D, 10–20 elements, one BC family from the list. With your AI pair: derive the Hermite element K and M matrices, assemble them, apply the BCs, solve the generalized eigenproblem, extract the first 4–6 modes. Animate one mode if you want the dopamine.
Spec this out with your AI pair — governing equation, discretization scheme, BCs. Have it draft the solver. Your job is to specify, audit, and verify.
A WASM build of one reference solver runs in your browser. Each row is one mode of the cantilever, oscillating; dots are the FEM nodes, dashed line is the undeformed neutral axis, hatched edge is the fixed end. Vary N and watch the curves and the error column move.
| Mode | FEM ω | Closed form (βnL)² | Error |
|---|
Two unknowns per node: deflection w and slope θ = ∂w/∂x. The weak form of a fourth-order operator needs C¹ continuity between elements — the slope has to be continuous across nodes, not just the displacement.
Hermite cubic shape functions are the standard choice — four
per element, exactly matching the four DOFs (w, θ at
each of the two end nodes). Closed-form 4×4 element K and M
matrices are in any structural-dynamics text (Bathe, Cook…).
For the generalized eigenproblem use a real library:
scipy.linalg.eigh(K, M) in Python,
nalgebra::SymmetricEigen after Cholesky reduction in
Rust.
Use the consistent mass matrix (the integrated ∫ ρA NTN dx), not lumped (diagonal, putting half the element mass at each end). Lumped is faster to code and noticeably wrong for higher modes — you'll see the third and fourth frequencies drift low under refinement.
def element_k(le, EI): ... # 4x4 Hermite stiffness
def element_m(le, rhoA): ... # 4x4 consistent mass
def assemble(n_el): ... # scatter into global K, M
def apply_bcs(K, M, bc): ... # drop fixed rows/cols
def solve_eig(K, M, n_modes):
# generalized symmetric eigenproblem
return scipy.linalg.eigh(K, M, subset_by_index=[0, n_modes-1])
def plot_modes(modes): ... # optional
scipy.linalg.eigh (Python) or
nalgebra::SymmetricEigen after Cholesky (Rust). Half a
century of LAPACK heuristics for handling near-degenerate eigenvalues,
ill-conditioning, and scaling are already in there. Knowing when to
stop reinventing them is part of the judgement this site is here to
grow.
The reference solver is the same Rust crate driving the live panel above — uniform cantilever, Hermite cubics, consistent mass, Cholesky reduction onto a standard symmetric eigenproblem. solver/src/lib.rs (~80 lines).
Your solver probably looks different in shape — that's expected. What to compare:
Post-mortem from the author: shipped with consistent mass from the start because lumped is the kind of shortcut that hides until rung 3 refinement embarrasses you. The boring choice is the right one.