Prev
Unit 10 · Boundary Value Problems and Partial Differential Equations
The Finite Element Method in One Dimension
70m100 XP
Next
🌡️Boundary Value Problems and Partial Differential Equations

The Finite Element Method in One Dimension

Energy and the Ritz method, the weak form and Galerkin, hat functions, element matrices and assembly, bars, beams and mass matrices

Göttingen, 1909
In 1787 Ernst Chladni, a musician and physicist from Wittenberg, scattered sand on a thin metal plate and stroked its edge with a violin bow. The plate sang, and the sand jumped away from the parts that vibrated and gathered along the lines that stayed still — a different pattern for every note. Chladni drew dozens of these figures and demonstrated them across Europe, in Paris in 1808. Explaining them took a century. Sophie Germain worked on the theory of vibrating plates for years after Chladni's visit to Paris. In 1850 Gustav Kirchhoff wrote down the correct equation — fourth order, with complicated conditions along a free edge — and could solve it only for a circular plate. For the square plates in Chladni's drawings, nobody could. Walther Ritz, a Swiss physicist, had been fighting tuberculosis since 1900. Back in Göttingen in 1908 after years of cures, he published a way around the equation. Do not try to satisfy the differential equation at every point. The vibrating plate takes the shape that makes its energy stationary; so write the shape as a combination of a few well-chosen functions, express the energy in terms of their coefficients, and minimise that — a problem in a handful of unknowns. With it, in 1909, he computed the Chladni figures of the free square plate and the pitches of their notes, and compared them with Chladni's measurements. The agreement was remarkable. Ritz died in Göttingen on 7 July 1909, aged 31. Göttingen, absorbed in existence proofs, hardly noticed his method. Engineers did: in 1915 Boris Galerkin in Petrograd recast it for rods and plates, and in 1941 Richard Courant proposed building the functions out of small local pieces. Ritz's idea with Courant's pieces is, almost exactly, the finite element method.
Take the simplest structure there is: an elastic bar along [0,L][0, L], axial stiffness EAEA, fixed at x=0x = 0, loaded by f(x)f(x) per unit length. There are two ways to say what its displacement u(x)u(x) is.
  • Equilibrium, point by point: −(EA u′)′=f-(EA\,u')' = f everywhere in the bar, with u(0)=0u(0) = 0 and no force at the free end, EA u′(L)=0EA\,u'(L) = 0.
  • Energy, all at once: among all displacements with u(0)=0u(0) = 0, the true one makes the total potential energy smallest,
Π(u)=∫0L12EA (u′)2 dx−∫0Lf u dx,\Pi(u) = \int_0^L \tfrac12 EA\,(u')^2\,dx - \int_0^L f\,u\,dx,
the strain energy stored minus the work done by the loads. You have met this principle twice already: solving Ax=bA\mathbf x = \mathbf b as minimising 12xTAx−bTx\tfrac12\mathbf x^{\mathsf T}A\mathbf x - \mathbf b^{\mathsf T}\mathbf x (lesson 4-3), and the draughtsman's spline settling into its least-energy shape (lesson 7-3).
Derivation
Minimising the energy gives the differential equation
Let uu be the minimiser and try a nearby displacement u+εvu + \varepsilon v, where vv is any function with v(0)=0v(0) = 0 (it must respect the support). Expand:
Π(u+εv)=Π(u)+ε[∫0LEA u′v′ dx−∫0Lf v dx]+ε22∫0LEA (v′)2 dx.\Pi(u + \varepsilon v) = \Pi(u) + \varepsilon\Big[\int_0^L EA\,u'v'\,dx - \int_0^L f\,v\,dx\Big] + \frac{\varepsilon^2}{2}\int_0^L EA\,(v')^2\,dx.
For uu to be the minimum, the term linear in ε\varepsilon must vanish — otherwise a small ε\varepsilon of the right sign would lower the energy. So for every admissible vv:
∫0LEA u′v′ dx=∫0Lf v dx.(∗)\int_0^L EA\,u'v'\,dx = \int_0^L f\,v\,dx. \qquad (\ast)
Integrate the left side by parts:
∫0LEA u′v′ dx=[EA u′ v]0L−∫0L(EA u′)′ v dx=EA u′(L) v(L)−∫0L(EA u′)′ v dx,\int_0^L EA\,u'v'\,dx = \Big[EA\,u'\,v\Big]_0^L - \int_0^L (EA\,u')'\,v\,dx = EA\,u'(L)\,v(L) - \int_0^L (EA\,u')'\,v\,dx,
using v(0)=0v(0) = 0. Then (∗)(\ast) becomes
∫0L[(EA u′)′+f] v dx−EA u′(L) v(L)=0for every v.\int_0^L \big[(EA\,u')' + f\big]\,v\,dx - EA\,u'(L)\,v(L) = 0 \quad\text{for every } v.
Choosing vv concentrated anywhere inside the bar forces the bracket to vanish: −(EA u′)′=f-(EA\,u')' = f, the equilibrium equation. Choosing v(L)≠0v(L) \ne 0 then forces EA u′(L)=0EA\,u'(L) = 0: the free end carries no force. The energy principle contains the equation and the free-end condition — nobody had to impose the second one. The ε2\varepsilon^2 term, half the strain energy of vv, is positive, which confirms a minimum.
💡The weak form and the Ritz–Galerkin method
Equation (∗)(\ast) is the weak form: 'for every virtual displacement vv, internal virtual work equals external virtual work' — the principle of virtual work. It contains only first derivatives, while the differential equation needs second ones, so it accepts functions with kinks. Ritz: choose a few functions φ1,…,φn\varphi_1, \dots, \varphi_n (each with φj(0)=0\varphi_j(0) = 0), look for u≈∑jcjφju \approx \sum_j c_j\varphi_j, and minimise Π\Pi over the coefficients cjc_j. Galerkin: require (∗)(\ast) only for v=φ1,…,φnv = \varphi_1, \dots, \varphi_n. For problems with an energy the two give the same equations; Galerkin's version also works when there is no energy to minimise (convection, for example). Either way the result is a linear system Kc=F,Kij=∫0LEA φi′φj′ dx,Fi=∫0Lf φi dx.K\mathbf c = \mathbf F, \qquad K_{ij} = \int_0^L EA\,\varphi_i'\varphi_j'\,dx, \qquad F_i = \int_0^L f\,\varphi_i\,dx.
Derivation
From the energy to Kc=FK\mathbf c = \mathbf F
Put u=∑jcjφju = \sum_j c_j\varphi_j into the energy. The derivative is u′=∑jcjφj′u' = \sum_j c_j\varphi_j', so (u′)2=∑i∑jcicjφi′φj′(u')^2 = \sum_i\sum_j c_ic_j\varphi_i'\varphi_j', and
Π(c)=12∑i∑jcicj∫0LEA φi′φj′ dx⏟Kij−∑ici∫0Lf φi dx⏟Fi=12 cTKc−cTF.\Pi(\mathbf c) = \tfrac12\sum_i\sum_j c_ic_j\underbrace{\int_0^L EA\,\varphi_i'\varphi_j'\,dx}_{K_{ij}} - \sum_i c_i\underbrace{\int_0^L f\,\varphi_i\,dx}_{F_i} = \tfrac12\,\mathbf c^{\mathsf T}K\mathbf c - \mathbf c^{\mathsf T}\mathbf F.
This is exactly the quadratic energy of lesson 4-3. KK is symmetric (swap ii and jj) and positive definite (it is twice the strain energy of a nonzero displacement), so the minimum is where the gradient vanishes: Kc=FK\mathbf c = \mathbf F. A differential equation has become a symmetric positive definite linear system — the kind Units 3 and 4 solve best.
💡Courant's idea: hat functions and elements
Ritz used global functions — polynomials or sines spread over the whole plate — so every KijK_{ij} was nonzero and each new problem needed new functions. Courant's choice makes the method automatic. Divide the bar into elements [xi−1,xi][x_{i-1}, x_i] and give each node a hat function: 1 at its own node, 0 at every other node, straight in between. Then
  • the coefficient cjc_j is simply the displacement at node jj;
  • a hat overlaps only its neighbours, so Kij=0K_{ij} = 0 unless ii and jj share an element: KK is sparse (here tridiagonal);
  • every integral splits into a sum over elements, each computed once from the same small recipe.
Inside one element of length ll, with ξ=(x−xi−1)/l\xi = (x - x_{i-1})/l running from 0 to 1, the displacement is u=ui−1N1(ξ)+uiN2(ξ)u = u_{i-1}N_1(\xi) + u_iN_2(\xi) with the shape functions N1=1−ξN_1 = 1 - \xi and N2=ξN_2 = \xi — the two halves of the hats that meet there.
Derivation
The bar element's stiffness and load
On one element, dN1/dx=−1/ldN_1/dx = -1/l and dN2/dx=+1/ldN_2/dx = +1/l (each changes by 1 over the length ll). The element's share of Kab=∫EA Na′Nb′ dxK_{ab} = \int EA\,N_a'N_b'\,dx is then a constant times the length:
ke=∫0lEA[−1/l1/l][−1/l1/l]dx=EAl[1−1−11].k_e = \int_0^l EA\begin{bmatrix} -1/l \\ 1/l \end{bmatrix}\begin{bmatrix} -1/l & 1/l \end{bmatrix}dx = \frac{EA}{l}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}.
That is a spring of stiffness EA/lEA/l between the two nodes — which is what a short bar is. A uniform load ww per metre gives
fe=∫0lw[1−ξξ]dx=wl2[11]:f_e = \int_0^l w\begin{bmatrix} 1 - \xi \\ \xi \end{bmatrix}dx = \frac{wl}{2}\begin{bmatrix} 1 \\ 1 \end{bmatrix}:
half of the element's load to each end. Assembly adds each kek_e into the rows and columns of its two nodes, and each fef_e into their entries — the COO accumulation of lesson 3-4. A node shared by two elements receives two contributions, so the global matrix has 2EA/l2EA/l on its interior diagonal and −EA/l-EA/l beside it: the tridiagonal matrix of lesson 3-3 again, now derived from energy rather than from difference quotients.
📐What converges — and what does not
A 30 m steel bar hanging under its own weight. The exact solution is u=wEA(Lx−x22)u = \frac{w}{EA}\big(Lx - \tfrac{x^2}{2}\big), with stress σ=w(L−x)/A\sigma = w(L - x)/A, largest at the support. One linear element already gives the exact tip displacement, 18.000 mm — in one dimension, linear elements with consistent loads are exact at the nodes for problems like this one. But stress is the derivative of displacement, constant within each element: the top element's stress is 120, 180, 210, 225 MPa with 1, 2, 4, 8 elements, against the exact 240 MPa at the support. A cantilever under uniform load with beam elements (below): the tip deflection is exact with one element (100.000 mm), but the root moment is 33.3, 38.3, 39.6 kN m for 1, 2, 4 elements against the exact 40.0. Displacements converge first; stresses and moments — the quantities that decide failure — converge more slowly.
💡Beam elements
A beam stores bending energy 12∫EI (w′′)2 dx\tfrac12\int EI\,(w'')^2\,dx, which contains second derivatives. A kinked deflection would have an infinite curvature at the kink, so the slope must be continuous between elements as well as the deflection. Each node therefore carries two unknowns, deflection ww and rotation θ=w′\theta = w', and each element uses the four cubic Hermite shape functions that match them: N1=1−3ξ2+2ξ3,N2=l(ξ−2ξ2+ξ3),N3=3ξ2−2ξ3,N4=l(−ξ2+ξ3).N_1 = 1 - 3\xi^2 + 2\xi^3,\quad N_2 = l(\xi - 2\xi^2 + \xi^3),\quad N_3 = 3\xi^2 - 2\xi^3,\quad N_4 = l(-\xi^2 + \xi^3). The element stiffness kab=∫0lEI Na′′Nb′′ dxk_{ab} = \int_0^l EI\,N_a''N_b''\,dx and the load for a uniform qq come out as ke=EIl3[126l−126l6l4l2−6l2l2−12−6l12−6l6l2l2−6l4l2],fe=ql12[6l6−l].k_e = \frac{EI}{l^3}\begin{bmatrix} 12 & 6l & -12 & 6l \\ 6l & 4l^2 & -6l & 2l^2 \\ -12 & -6l & 12 & -6l \\ 6l & 2l^2 & -6l & 4l^2 \end{bmatrix}, \qquad f_e = \frac{ql}{12}\begin{bmatrix} 6 \\ l \\ 6 \\ -l \end{bmatrix}.
Derivation
One entry of the beam stiffness, worked out
Check the Hermite functions first: at ξ=0\xi = 0, N1=1N_1 = 1 and the others vanish; at ξ=1\xi = 1, N3=1N_3 = 1 and the others vanish; and dN2/dx=1dN_2/dx = 1 at ξ=0\xi = 0, dN4/dx=1dN_4/dx = 1 at ξ=1\xi = 1, with all other slopes zero at the ends. So w=w1N1+θ1N2+w2N3+θ2N4w = w_1N_1 + \theta_1N_2 + w_2N_3 + \theta_2N_4 takes the right values and slopes at both nodes.
For the corner entry, differentiate twice with respect to x=lξx = l\xi (each d/dxd/dx brings a factor 1/l1/l):
N1′′=−6+12ξl2,k11=∫0lEI (N1′′)2 dx=EIl4∫01(12ξ−6)2 l dξ=EIl3(48−72+36)=12EIl3.N_1'' = \frac{-6 + 12\xi}{l^2}, \qquad k_{11} = \int_0^l EI\,(N_1'')^2\,dx = \frac{EI}{l^4}\int_0^1(12\xi - 6)^2\,l\,d\xi = \frac{EI}{l^3}(48 - 72 + 36) = \frac{12EI}{l^3}.
The other fifteen entries follow the same way. The load vector is ∫q Na dx\int q\,N_a\,dx; for example ∫0lqN2 dx=ql2∫01(ξ−2ξ2+ξ3) dξ=ql2(12−23+14)=ql212\int_0^l qN_2\,dx = ql^2\int_0^1(\xi - 2\xi^2 + \xi^3)\,d\xi = ql^2\big(\tfrac12 - \tfrac23 + \tfrac14\big) = \tfrac{ql^2}{12}. These are the matrices supplied in lesson 6-4's footbridge lab and used in this lesson's cantilever.
Derivation
Mass matrices, and why consistent mass overestimates frequencies
Consistent mass. Kinetic energy is 12∫ρA u˙2 dx\tfrac12\int\rho A\,\dot u^2\,dx. Interpolating the velocity with the same shape functions gives a mass matrix exactly as the strain energy gave KK: mab=∫ρA NaNb dxm_{ab} = \int\rho A\,N_aN_b\,dx. For the bar element, ∫01(1−ξ)2dξ=13\int_0^1(1 - \xi)^2d\xi = \tfrac13 and ∫01ξ(1−ξ) dξ=16\int_0^1\xi(1 - \xi)\,d\xi = \tfrac16, so
me=ρAl6[2112].m_e = \frac{\rho Al}{6}\begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}.
The beam element's consistent mass, by the same integrals with the Hermite functions, is μl420\frac{\mu l}{420} times a 4×44 \times 4 matrix with 156156 and 4l24l^2 on its diagonal — the one in lesson 6-4. The lumped alternative simply puts half of each element's mass at each node, ρAl2I\frac{\rho Al}{2}I: diagonal and cheap.
An upper bound. The lowest natural frequency squared is the smallest value of the Rayleigh quotient, strain energy over kinetic energy (lessons 6-1 and 6-2), over all possible shapes. The finite element model with consistent mass searches over a smaller set — only the shapes its elements can represent — so the minimum it finds can only be higher. Every computed frequency is above the true one, and refining the mesh brings it down. Lumped mass is not a Ritz approximation and usually errs below; the lab finds the first mode of a steel bar 0.64% high with consistent mass and 0.64% low with lumped mass on four elements.
When a change occurs in Nature, the quantity of action necessary for that change is as small as possible.
– Pierre Louis Maupertuis, Les Loix du Mouvement et du Repos (1746)
🕰From Galerkin to Boeing
Boris Galerkin, a structural engineer, published the version without energy — require the weak form for each basis function — in 1915, in a paper on rods and plates in the Vestnik Inzhenerov in Petrograd; in Russian texts it is often called the Bubnov–Galerkin method, after Ivan Bubnov's closely related work of 1913. Richard Courant, addressing the American Mathematical Society on 3 May 1941, proposed functions that are linear on each triangle of a mesh, and showed a torsion calculation with them — the first finite element computation that historians have found (lesson 10-6). Fifteen years later Turner, Clough, Martin and Topp at Boeing reached the same method from the engineering side (lesson 3-4), and in 1960 Clough named it.
⚠Finite element models are too stiff, never too flexible
Because the method minimises energy over a restricted set of shapes, a displacement-based finite element model can only be stiffer than the real structure: under a given load its displacements are, overall, too small and its natural frequencies too high. Coarse meshes therefore err on the unsafe side for deflection limits and resonance checks. (Stresses can be wrong in either direction locally — the hanging bar's top element is 25% low with two elements.)
💡Where you will use this
Structures: frame and truss programs are exactly these bar and beam elements, assembled in two and three dimensions. Rocket Propulsion and Aircraft Engines: rotor and shaft models for critical speeds, with consistent mass matrices. Spacecraft structures and Attitude Control: beam models of booms and solar arrays, reduced to a few modes for the controller. Geotechnics: piles and retaining walls as beams on springs.
📚Further reading
  • Martin Gander and Gerhard Wanner, 'From Euler, Ritz, and Galerkin to modern computing', SIAM Review 54(4), 627–666 (2012) — the full story of the cold open, with Ritz's own Chladni figures.
  • Gilbert Strang and George Fix, An Analysis of the Finite Element Method (Prentice-Hall, 1973; 2nd ed., Wellesley–Cambridge, 2008) — the Ritz–Galerkin view of finite elements, from one dimension up.
  • Thomas J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis (Prentice-Hall, 1987; Dover reprint 2000), chapter 1 — the one-dimensional model problem, weak form and element matrices.
↓But there is a problem
In one dimension every element is a line segment and every mesh is a row of them. The method's real power — and its real danger — appears in two and three dimensions, where elements are triangles and bricks fitted to curved walls, holes and fillets by a mesh generator, and the analyst must judge whether the mesh is fine enough where it matters. On 23 August 1991 the concrete base of the Sleipner A gas platform, a structure of twenty-four linked cells taller than a twenty-storey building, was being lowered in a controlled test in a Norwegian fjord. A wall cracked, water poured in, and the whole base sank to the floor of the fjord in minutes. The cause was traced to its computer model: a finite element analysis that had underestimated the stresses in the walls by 47%.