A Hands-On Python Study
When an engineer asks “how much will this beam sag?”, the textbook answer usually comes from integrating a differential equation. But there is a more physical way to arrive at the same shape: nature settles into the deflection that makes the total potential energy as small as possible. This is the Principle of Minimum Potential Energy, and it is the theoretical backbone of the Finite Element Method used throughout structural engineering.
In this article we take a concrete example — a simply supported beam under a uniform load — and find its deflection curve three different ways: an exact analytical formula, a Rayleigh-Ritz energy minimization with only two unknowns, and a full finite-element energy minimization with dozens of unknowns. All three will agree almost perfectly, which is exactly the point: minimizing energy really does reproduce the equilibrium shape of the structure.
1. The Physics: Total Potential Energy of a Bent Beam
For an Euler-Bernoulli beam with flexural rigidity $EI$, deflection $w(x)$, and a distributed load $q(x)$, the strain energy stored by bending is
$$
U[w] = \int_0^L \frac{EI}{2}\left(\frac{d^2 w}{dx^2}\right)^2 dx
$$
and the work done by the external load as the beam deflects is
$$
W[w] = \int_0^L q(x), w(x), dx
$$
The total potential energy functional is
$$
\Pi[w] = U[w] - W[w]
$$
The Principle of Minimum Potential Energy states that among all deflection shapes $w(x)$ satisfying the geometric boundary conditions, the true equilibrium shape is the one that minimizes $\Pi[w]$. Setting the first variation $\delta \Pi = 0$ recovers the familiar governing equation
$$
EI,\frac{d^4 w}{dx^4} = q(x)
$$
So instead of solving this fourth-order ODE directly, we can instead search for the function that minimizes $\Pi[w]$ — which is exactly what we will do numerically below.
2. The Example Problem
Simply supported steel beam:
$$
L = 4\ \text{m}, \qquad E = 210\ \text{GPa}, \qquad I = 8.5\times 10^{-6}\ \text{m}^4, \qquad q_0 = 15\ \text{kN/m (uniform)}
$$
The exact Euler-Bernoulli solution for this classic case is
$$
w(x) = \frac{q_0,x\left(L^3 - 2Lx^2 + x^3\right)}{24,EI}, \qquad
w_{\max} = w!\left(\frac{L}{2}\right) = \frac{5,q_0 L^4}{384,EI}
$$
We will use this closed-form result as the “ground truth” to check our two energy-minimization approaches against.
3. Two Ways to Minimize the Energy
Rayleigh-Ritz (few unknowns). We approximate the deflection with a truncated sine series that already satisfies $w(0)=w(L)=0$ and $w’’(0)=w’’(L)=0$:
$$
w(x) \approx a_1 \sin\frac{\pi x}{L} + a_2 \sin\frac{3\pi x}{L}
$$
Substituting into $\Pi[w]$ and using the orthogonality of sine functions turns the functional minimization into an ordinary 2-variable minimization of $\Pi(a_1,a_2)$, which scipy.optimize.minimize solves in a handful of iterations.
Finite-element energy minimization (many unknowns). For a more general, higher-fidelity picture we discretize the beam into 40 two-node Hermite beam elements. Because $\Pi$ is a quadratic function of the nodal deflections/rotations, its minimum is found by solving one linear system $\mathbf{K}\mathbf{a} = \mathbf{F}$ instead of running an iterative optimizer — this is dramatically faster and is, in fact, exactly what commercial FEM software does under the hood.
4. Full Source Code (Google Colab, single cell)
1 | import numpy as np |
5. Code Walkthrough
Section 1 – Parameters. Standard steel beam constants are defined once and reused everywhere via keyword arguments, so changing L, EI, or q0 at the top automatically propagates through every method below.
Section 2 – Analytical solution. w_analytical implements the closed-form Euler-Bernoulli formula directly; this becomes our reference for measuring accuracy of the two energy-based methods.
Section 3 – Rayleigh-Ritz. The key trick is that once you plug the sine trial function into $\Pi[w]$, the integrals collapse (thanks to sine orthogonality) into a simple closed-form quadratic in $a_1, a_2$ — the potential_energy function evaluates it in vectorized NumPy without ever performing a numerical integration. scipy.optimize.minimize (BFGS) then locates the minimum in only a handful of iterations, because the surface is a smooth convex bowl. The same closed-form expression is re-evaluated over a 2D meshgrid to build the energy-landscape surface for plotting — this is fully vectorized, so even a $120\times120$ grid evaluates instantly.
Section 4 – Finite-element energy minimization. This is the “fast version” for a large number of unknowns. Rather than handing scipy an 80+ dimensional nonlinear minimization (which would be slow and could converge poorly), we exploit the fact that $\Pi$ is exactly quadratic in the nodal unknowns for this linear-elastic problem. That means its minimizer is the unique solution of the linear system $\mathbf{K}\mathbf{a}=\mathbf{F}$, where $\mathbf{K}$ is the assembled global stiffness matrix and $\mathbf{F}$ the assembled load vector — solved instantly with np.linalg.solve. The element stiffness matrix and consistent load vector used here are the standard 2-node Euler-Bernoulli beam (Hermite cubic) formulas. Boundary conditions are applied simply by deleting the rows/columns corresponding to the two deflection DOFs that are pinned (fixed), leaving the rotations free, exactly matching a simply-supported beam.
Section 5 – Deflection surface. w_analytical is called with a 2D array of load values via NumPy broadcasting, producing the whole family of deflection curves for varying $q_0$ in one vectorized call — no Python loop needed.
Section 6 – Plots. Three figures are generated: the Rayleigh-Ritz energy bowl in 3D, a 2D comparison of all three deflection curves, and a 3D deflection-vs-load surface.
=== Beam parameters === L = 4.0 m EI = 1.7850e+06 N*m^2 q0 = 1.5000e+04 N/m Analytical maximum deflection at midspan: 28.0112 mm === Rayleigh-Ritz result === a1 = 2.811924e-02 m, a2 = 1.157096e-04 m Converged: True, iterations: 5 Max relative error of Ritz solution: 0.0379 % === FEM (energy minimization) result === Number of elements: 40 Max relative error of FEM solution: 0.074233 % FEM max deflection: 28.0112 mm All computations finished successfully.
6. Reading the Graphs
Figure 1 — the energy landscape. This 3D surface plots $\Pi(a_1,a_2)$ over a range of trial coefficients. It should look like a smooth, slightly elongated bowl (paraboloid), since $\Pi$ is quadratic. The red marker sits exactly at the bottom of the bowl — this is the point the optimizer converged to, and it is the physical deflection shape the beam actually adopts. Everywhere else on the surface represents a kinematically possible but non-equilibrium shape that would require some external agent to hold it there; only at the bottom does the structure balance on its own. This single picture is the clearest visual proof of the minimum-potential-energy principle: equilibrium literally sits at the lowest point of the energy surface.
Figure 2 — deflection curve comparison. Three predictions are overlaid: the exact analytical curve (solid black), the 2-term Rayleigh-Ritz curve (dashed blue), and the finite-element nodal solution (red dots). All three should be essentially indistinguishable, sagging maximally at midspan (around 28 mm for the given parameters) and returning to zero at both supports. The console output prints the relative error of each approximate method against the exact solution — expect the Ritz error to be under 0.1% and the FEM error to be smaller still, which shows that even a very small number of well-chosen basis functions, or a modest FEM mesh, is enough to essentially reproduce the exact elasticity solution once you frame the problem as energy minimization.
Figure 3 — deflection surface across load intensities. This second 3D plot sweeps the uniform load $q_0$ from light to heavy and stacks the resulting deflection curves into one continuous surface. Position along the beam runs along one horizontal axis, load intensity along the other, and deflection is the height. Because $w(x)$ is linear in $q_0$ for a linear-elastic beam, the surface is a smooth, uniformly tilted “wedge” — every cross-section at a fixed $x$ is a straight ramp in $q_0$, and every cross-section at a fixed $q_0$ reproduces the familiar sagging curve from Figure 2. This surface is a convenient way to read off, at a glance, how sensitive the beam’s midspan deflection is to load magnitude without re-running the calculation for each case.
7. Takeaways
Framing beam bending as an energy-minimization problem does more than reproduce the textbook formula — it generalizes cleanly. The same $\Pi[w] = U[w] - W[w]$ recipe extends directly to beams with variable cross-sections, non-uniform loads, multiple spans, or elastic foundations, simply by changing what goes into the integrals. The Rayleigh-Ritz method shows the idea in its purest, lowest-dimensional form, while the finite-element approach shows how the same principle scales up to realistic structural models — and, because the energy is quadratic, always boils down to nothing more exotic than solving $\mathbf{K}\mathbf{a}=\mathbf{F}$.
