A different question to ask the equation
The earlier guides in this rung all asked the same kind of question: make the PDE hold, approximately, at the points of a grid. Replace each derivative by a stencil, demand u_t = k u_xx at every crossing, and march. That works beautifully on a rectangle, but real problems live on awkward shapes — a turbine blade, a coastline, a region with a hole in it — where a tidy grid simply will not fit. The finite element method is born from asking the equation a different question entirely, one that bends gracefully to any shape.
Take the model elliptic problem of this whole ladder, Poisson's equation -Laplacian u = f on a region, with u = 0 on the boundary. Instead of insisting the equation hold pointwise, multiply both sides by an arbitrary test function v that also vanishes on the boundary, and integrate over the region. A single integration by parts then moves one derivative off u and onto v, and the boundary term dies because v = 0 there. What survives is the integral of grad u dotted with grad v equals the integral of f times v — and crucially, this must hold for every admissible v. That is the weak formulation: the equation no longer asked at points, but tested against a whole space of functions.
strong form: -Laplacian u = f, u = 0 on boundary
multiply by v, integrate, integrate by parts once:
weak form: integral( grad u . grad v ) = integral( f v )
must hold for EVERY test function v
a(u, v) = L(v) for all v
\___/ \__/
bilinear linearWhy the weak form is the right thing to discretise
Two gifts come from this rewrite. First, the weak form asks less of the solution: only one derivative of u appears, not two, so u no longer needs to be twice differentiable. It only needs a single derivative that is square-integrable — exactly the membership condition of the space H^1_0, the natural home of functions with finite energy that vanish on the boundary. This is the same relaxation you met when weak formulations let kinked, corner-bearing functions count as honest solutions. Second, the left side is a bilinear form a(u,v) — linear in each slot separately — and the right side is a linear functional L(v). The whole problem becomes: find u in H^1_0 so that a(u,v) = L(v) for all v.
Why does such a u exist, and why is it unique? Here is one of the quiet triumphs of the theory: the Lax-Milgram theorem guarantees it, provided the bilinear form is bounded and coercive — meaning a(v,v) is at least a positive multiple of the energy of v, so the form cannot collapse in any direction. For Poisson's equation a(v,v) is literally the integral of |grad v|^2, the Dirichlet energy, and coercivity follows from a Poincare inequality. There is even a beautiful physical reading: among all functions with the right boundary values, the true solution is the one that minimises that energy — this is Dirichlet's principle, and the weak form is just the statement that its first variation vanishes.
Tent functions: making the space finite
The weak problem still lives in an infinite-dimensional space H^1_0 — we cannot search all of it on a computer. The decisive step, the Galerkin method, is to replace that vast space by a small, finite-dimensional subspace V_h spanned by a handful of simple, chosen functions, and to demand a(u_h, v) = L(v) only for v inside V_h. Geometrically you are projecting the true solution onto a manageable slab of the function space; the smaller the slab, the cruder the answer, but the more refined V_h is, the closer u_h presses against the truth.
What makes it the finite element method is the clever choice of basis for V_h. Chop the region into small simple pieces — intervals in 1D, triangles in 2D — a mesh that can hug any shape, fine where the solution varies fast and coarse where it is lazy. On this mesh put one basis function phi_j per node: the hat function (or tent function) that equals 1 at node j and slides linearly down to 0 at every neighbouring node, so each phi_j is nonzero only on the few elements touching node j. Any candidate u_h is then a weighted sum, u_h = sum of c_j phi_j, and the unknowns are just the nodal values c_j — the height of the surface at each mesh point.
Assembling the stiffness matrix
Now turn the Galerkin condition into linear algebra. It suffices to test against each basis function, so plug u_h = sum of c_j phi_j into a(u_h, phi_i) = L(phi_i) for every node i. Because a is bilinear, the sum comes out and you are left with sum over j of a(phi_j, phi_i) c_j = L(phi_i). Collect the coefficients a(phi_j, phi_i) into a matrix K with entries K_ij = integral of grad phi_i dotted with grad phi_j — the famous stiffness matrix, named from its origin in structural mechanics. Collect the right-hand sides into a vector F with F_i = integral of f phi_i. The continuous PDE has become the linear system K c = F.
Here is where the local support of the tents pays off enormously. K_ij is nonzero only when phi_i and phi_j overlap — that is, only when nodes i and j share an element. So in 1D each row of K has just three nonzero entries, and in 2D a handful: the matrix is sparse, mostly zeros, and the same fast solvers that handled implicit schemes in the previous guide make short work of it. Better still, K is symmetric and positive-definite for Poisson, a direct inheritance of the coercive, symmetric bilinear form, which guarantees the system has exactly one solution and lets the most efficient solvers loose.
- Mesh the region: split it into elements (intervals or triangles), refining where you expect the solution to vary fast.
- Choose basis functions: one hat function phi_j per interior node, each linear on every element and 1 at its own node, 0 at the others.
- Build the stiffness matrix K with K_ij = integral of grad phi_i . grad phi_j, computed element by element — most entries are zero.
- Build the load vector F with F_i = integral of f phi_i, usually by a small quadrature rule on each element.
- Solve the sparse linear system K c = F; the coefficients c are the nodal values, and u_h = sum of c_j phi_j is the approximate solution.
How good is it? The error estimate
An honest method must say how wrong it is. The single most important fact about the Galerkin solution is Cea's lemma: u_h is, up to a fixed constant, the best possible approximation to the true u that the subspace V_h can offer, measured in the energy norm. In words, the method does not merely produce some answer in V_h — it produces essentially the closest point in V_h to the truth. That reduces the whole error question to a pure approximation question: how well can piecewise-linear tents on a mesh of size h approximate a smooth function?
For linear elements, interpolation theory answers cleanly: the error in the energy norm (which controls the first derivative) shrinks like C h, first-order in the mesh size, while the error in the plain L^2 norm (the size of u_h - u itself) shrinks like C h^2, second-order — the celebrated extra power coming from a duality argument. As you halve the mesh spacing, the function values roughly quarter their error. The hidden constant C again hides higher derivatives of the true solution, so the estimate is honest only when u is smooth enough; near a re-entrant corner, where u is merely H^1, the rate genuinely degrades, and this is exactly the symptom that adaptive mesh refinement chases — pour small elements only where the error indicator flags trouble.
Step back and notice how different the bookkeeping felt from finite differences. There was no Taylor expansion, no stencil, no per-point truncation error; convergence flowed instead from coercivity (existence and stability, via Lax-Milgram), the best-approximation property (Cea), and interpolation (how well tents fit smooth functions). It is the same trinity of consistency, stability, and convergence you saw in the Galerkin picture and in the equivalence theorem two guides ago — just dressed in the language of energy and function spaces rather than grids and steps. The next guide leaves elliptic comfort behind and confronts the hardest case of all: conservation laws that develop shocks, where finite volumes and entropy conditions take over.