The finite element method in one dimension

Most problems in structural mechanics, heat transfer and fluid flow are boundary value problems: a partial differential equation must hold everywhere inside a domain, and some conditions must hold on its boundary. For a prismatic beam or a thin-walled cylinder you can solve these on paper. For a bracket with a fillet, a boss and three bolt holes, you can’t.

The finite element method (FEM) gets around this in three moves:

1. Approximate. Cut the domain into small, simple pieces (elements) and describe the unknown field on each piece with a low-order polynomial, defined by its values at a few points (nodes).
2. Discretise. Instead of asking the differential equation to hold at every point, which a piecewise polynomial cannot do, ask it to hold on average (the weak form), weighted by a set of test functions. That turns the PDE into a linear system $\mathbf{K}\mathbf{u} = \mathbf{f}$.
3. Solve the linear system, then recover derived quantities such as strain and stress from $\mathbf{u}$.

The figure shows all three on a real part. Switch between geometry, mesh and stress. Using the method on real parts like this one is finite element analysis (FEA), and that is what this track teaches.

This chapter does all three for the simplest possible case: a straight bar under axial load. That is not a toy. Every idea you need for a 3D solid model (weak form, shape functions, element stiffness, assembly, supports, convergence) appears here unchanged. Going to 3D only adds bookkeeping.

Take a straight bar of length $L$, Young’s modulus $E$ and cross-section $A$, fixed at $x=0$. It carries a distributed axial load $q(x)$ (force per unit length, say its own weight hanging vertically) and a point load $F$ at the free end.

Cut out a slice of length $dx$. Equilibrium of the normal force $N$ across it gives $dN/dx + q = 0$. Hooke’s law with the small-strain definition $\varepsilon = du/dx$ gives $N = EA\, du/dx$. Combine them:

$$-\frac{d}{dx}\!\left(EA\,\frac{du}{dx}\right) = q(x), \qquad 0 < x < L \tag{1.1}$$

with two boundary conditions of different character:

• $u(0) = 0$ prescribes the unknown itself. This is an essential (Dirichlet) condition: a support.
• $EA\,u'(L) = F$ prescribes its derivative, which here is a force. This is a natural (Neumann) condition: a load.

For constant $EA$ and $q$, integrate twice:

$$\begin{aligned} N(x) &= F + q\,(L - x), \\[4pt] u(x) &= \frac{F x}{EA} + \frac{q}{EA}\left(Lx - \frac{x^2}{2}\right) \end{aligned} \tag{1.2}$$

This is the strong form: it asks $u$ to be twice differentiable and the equation to hold at every point. Keep this solution in mind; we’ll measure the finite element answer against it.

We want to approximate $u$ with something a computer can store: a finite list of numbers. The simplest continuous choice is a piecewise-linear function: split $[0, L]$ into elements of length $h$ and connect the values at the nodes with straight lines.

Now try to plug this $u_h$ into the strong form. Inside each element $u_h'' = 0$. At each node the slope jumps, so $u_h''$ is not a function at all: it is a Dirac delta whose strength equals the jump. The differential equation can’t be satisfied, not even approximately, in any pointwise sense.

There are two ways out. Use smoother functions (cubic splines, for example), which is costly and awkward in 2D and 3D. Or rewrite the problem so that it needs only first derivatives. FEM takes the second route.

Multiply the differential equation by an arbitrary test function $v(x)$ and integrate over the bar. Then integrate the left-hand side by parts to move one derivative from $u$ onto $v$:

$$\begin{gathered} \int_0^L EA\,u'\,v'\,dx \;=\; \int_0^L q\,v\,dx \;+\; F\,v(L) \\[2pt] \text{for all } v \text{ with } v(0) = 0 \end{gathered} \tag{1.3}$$

Three things happened:

• Only first derivatives remain. Piecewise-linear functions now qualify: their first derivative is piecewise constant, and the integral of a product of piecewise constants is perfectly well defined.
• The natural condition appeared by itself. The boundary term from integration by parts is $[EA\,u'\,v]_0^L$. At $x=L$, $EA\,u' = F$, which is where the $F\,v(L)$ on the right comes from. Loads enter the weak form automatically. That is why they are called natural.
• The essential condition did not. We must build $u(0)=0$ into the space of candidate functions and choose $v(0)=0$ so the unknown reaction at the support drops out.

Mechanically, this is the principle of virtual work. Read $v$ as a virtual displacement compatible with the supports. The left side is the internal virtual work of the normal force; the right side is the virtual work of the applied loads. Equilibrium means the two match for every admissible $v$. Equivalently, $u$ minimises the total potential energy

$$\Pi[u] = \tfrac12\int_0^L EA\,(u')^2\,dx - \int_0^L q\,u\,dx - F\,u(L). \tag{1.4}$$

For a smooth $u$ the strong and weak forms are equivalent. The weak form simply admits a larger family of candidate solutions, and that larger family includes our piecewise polynomials.

The Galerkin method restricts both the solution and the test functions to the same finite-dimensional space. For linear elements the natural basis of that space is the set of hat functions, the shape functions $N_i(x)$:

$$u_h(x) = \sum_{i=1}^{n} N_i(x)\,u_i \tag{1.5}$$

Each $N_i$ equals 1 at node $i$, 0 at all other nodes, and is linear in between. The coefficients $u_i$ are therefore not abstract weights: they are the nodal displacements, which is why FE output is “displacement at each node”.

Inside one element of length $h$, with local coordinate $\xi = (x - x_1)/h \in [0,1]$, only two hats are non-zero:

$$N_1(\xi) = 1 - \xi, \qquad N_2(\xi) = \xi \tag{1.6}$$

Three properties carry over to every element type you’ll meet later:

• Interpolation: $u_h(x_i) = u_i$.
• Partition of unity: $\sum_i N_i(x) = 1$, so a rigid translation (all $u_i$ equal) is represented exactly and produces zero strain.
• Compact support: $N_i$ is non-zero only on the elements touching node $i$. Two nodes interact only if they share an element. This is what makes the final matrix sparse.

Substitute the expansion (1.5) into the weak form (1.3) and choose the test function to be each basis function in turn, $v = N_j$. Every choice of $j$ gives one equation:

$$\begin{aligned} \sum_i K_{ji}\,u_i &= f_j, \\[2pt] K_{ji} &= \int_0^L EA\,N_i'\,N_j'\,dx, \\[2pt] f_j &= \int_0^L q\,N_j\,dx + F\,N_j(L). \end{aligned} \tag{1.7}$$

That is the linear system $\mathbf{K}\mathbf{u}=\mathbf{f}$: one row per test function, one column per nodal unknown. Because each $N_i$ lives on only two elements, the integral splits into a sum of per-element integrals. On one element $N_1' = -1/h$ and $N_2' = 1/h$, so

$$\mathbf{k}_e = \frac{EA}{h}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}, \qquad \mathbf{f}_e = \frac{q\,h}{2}\begin{bmatrix} 1 \\ 1 \end{bmatrix}. \tag{1.8}$$

Two observations worth keeping:

• $\mathbf{k}_e$ is the stiffness matrix of a spring with $k = EA/h$. Shorter elements are stiffer springs.
• $\mathbf{k}_e$ is singular: $\mathbf{k}_e [1, 1]^\mathsf{T} = \mathbf{0}$. A free element can translate without resistance. Every element matrix in every FE code has a null space made of its rigid-body motions (one here, three in 2D, six in 3D).

The distributed load is split equally between the two nodes. This is the consistent load vector, and it is what you get from the integral, not a convenient guess.

The global matrix is the sum of the element matrices, each scattered into the rows and columns of its own nodes. This is the direct stiffness method, or assembly, and it fits in a few lines:

K = zeros((n_nodes, n_nodes))
f = zeros(n_nodes)
for (i, j) in elements:            # node pair of each element
    dofs = [i, j]
    K[ix_(dofs, dofs)] += k_e      # 2×2 element matrix
    f[dofs]           += f_e

The assembled $\mathbf{K}$ is symmetric, banded (node $i$ only couples to $i\pm1$), and singular. The whole bar can still slide, so the supports have to go in before solving. The simplest way to impose $u_1 = 0$ is to delete row and column 1. What remains is symmetric positive definite and has a unique solution.

In a 3D solid model nothing changes except scale: three unknowns per node, a few hundred thousand to a few million rows, and a sparsity pattern set by the mesh connectivity. Sparse storage and sparse solvers are the only reason this is affordable.

Solve the 1D system and compare with the exact solution (1.2). Two results stand out.

The nodal displacements are exact, whatever the mesh. That is a special property of 1D problems with linear elements, and you should not expect it anywhere else. Between the nodes the straight lines miss the parabola by at most $q h^2 / 8EA$.

The normal force is worse. $N_h = EA\,u_h'$ is constant inside each element, so it approximates a linear function with a staircase. Its error is $q h/2$: halve the elements and you only halve the error. The staircase is exact at element midpoints. These superconvergent points are why FE codes evaluate stresses at interior integration points and then extrapolate or average to the nodes.

The general version of this observation holds in every dimension. For smooth solutions and elements of polynomial degree $p$:

$$\| u - u_h \| \le C\,h^{p+1}, \qquad \| \sigma - \sigma_h \| \le C\,h^{p} \tag{1.9}$$

Derived quantities converge one order slower than the primary unknown. In practice: a mesh that gives good displacements may still give poor stresses. Quadratic elements ($p = 2$) buy one extra order for both.

In plane stress the unknown becomes a vector field $(u_x, u_y)$ and the bar’s single strain becomes three: $\boldsymbol\varepsilon = (\varepsilon_x, \varepsilon_y, \gamma_{xy})$. The recipe does not change:

1. Interpolate each displacement component with shape functions on each element.
2. Strains are derivatives of displacements, so $\boldsymbol\varepsilon = \mathbf{B}\,\mathbf{u}_e$ with $\mathbf{B}$ built from $\partial N_i/\partial x$ and $\partial N_i/\partial y$.
3. Stresses follow from the material law: $\boldsymbol\sigma = \mathbf{D}\,\boldsymbol\varepsilon$.
4. The element stiffness is the same weak-form integral in matrix dress:

$$\mathbf{k}_e = \int_{A_e} \mathbf{B}^\mathsf{T}\,\mathbf{D}\,\mathbf{B}\; t\,dA \tag{1.10}$$

5. Assemble, apply supports, solve.

With the 3-node linear triangle, $\mathbf{B}$ is constant, so each element carries a single, constant strain. In bending, strain varies linearly through the depth. A constant-strain element can only follow that variation if there are many of them through the thickness. With few, the elements resist bending in ways the real material doesn’t, and the model comes out too stiff. Try it in the figure: one triangle through the depth recovers less than a third of the true deflection.

This is the first lesson that has no 1D counterpart: element choice matters as much as element size. It is why commercial codes default to quadratic elements, and it gets its own rule in the field notes.