From an infinite menu to a finite one
The previous guide rewrote our model problem — say -u''(x) = f(x) on the interval [0, 1] with u(0) = u(1) = 0 — in its weak form: find a function u so that the integral of u'(x)*v'(x) equals the integral of f(x)*v(x), for every admissible test function v. That single line is exact, but it asks for u inside an infinite-dimensional space of functions — far too big to store. The whole finite element method is one disciplined compromise: stop searching the infinite space, and search only a finite-dimensional subspace we can actually write down with numbers.
This is the Galerkin method from the last guide, now made concrete. We pick a finite set of building-block functions phi_1, phi_2, ..., phi_N — the basis — and insist the approximate solution be a combination u_h(x) = c_1*phi_1(x) + ... + c_N*phi_N(x). The unknowns are no longer a whole function but just N numbers, the coefficients c_1, ..., c_N. Galerkin then says: require the weak form to hold not for every test function, but only when v is each basis function phi_i in turn. That is N equations for N unknowns — a finite linear system we can hand to the solvers from earlier rungs.
The mesh: chopping the domain into elements
Before we can build local functions we need somewhere to hang them. Cut the interval [0, 1] at nodes 0 = x_0 < x_1 < ... < x_N < x_{N+1} = 1. Each little sub-interval [x_{j}, x_{j+1}] is an element, and the whole collection of nodes and elements is the finite element mesh. In 1D a mesh is just a string of beads; in 2D it is usually a patchwork of triangles tiling a region, and in 3D a packing of tetrahedra. The mesh need not be uniform — that freedom is the method's superpower. You can crowd tiny elements where the solution is sharp (a crack tip, a boundary layer) and let them sprawl where it is calm, fitting awkward shapes a regular grid could never follow.
Now the basis. Over this mesh, attach to each interior node x_j one shape function phi_j: a continuous tent (or hat) function that equals 1 at its own node x_j, falls linearly to 0 at the two neighbouring nodes, and is exactly 0 everywhere else. Picture a circus tent pitched at x_j and pegged down at x_{j-1} and x_{j+1}. Because each tent is nonzero over only its two adjacent elements, the basis functions barely overlap — and that locality is the whole game. A function u_h built from these tents is automatically a piecewise-linear interpolant: straight-line segments between nodes, with c_j the height (and, conveniently, the value u_h(x_j)) at node j.
Where the matrix comes from
Substitute u_h = c_1*phi_1 + ... + c_N*phi_N into the weak form and let v run through each phi_i. The integral of u_h'*phi_i' becomes the sum over j of c_j times the integral of phi_j'*phi_i'. Collect those integrals into a matrix A with entries A_ij = integral of phi_i'(x)*phi_j'(x) dx, and the right-hand side into a vector b with b_i = integral of f(x)*phi_i(x) dx. The infinite-dimensional weak problem has collapsed into the familiar finite system A x = b, where x = (c_1, ..., c_N) holds the nodal values. The matrix A is the stiffness matrix — a name borrowed from its origin in structural engineering, where A_ij really measured how stiffly node i resists a push at node j.
Here the local tents pay off spectacularly. The entry A_ij is an integral of phi_i'*phi_j', and that product is nonzero only where both tents are nonzero — only when nodes i and j sit in a shared element, i.e. when they are neighbours. So A_ij = 0 unless |i - j| <= 1. In 1D the stiffness matrix is tridiagonal; in 2D each node touches only its handful of mesh neighbours, giving a beautifully sparse matrix with a few nonzeros per row out of N. We never store, nor want, a dense N-by-N array. For the uniform mesh with spacing h, a short calculation gives the diagonal 2/h and each off-diagonal -1/h — and dividing through by h you may notice it is exactly the (-1, 2, -1) finite-difference stencil from the previous rung. The two methods meet on this simple problem; in general they differ, but the kinship is real.
weak form, with u_h = sum_j c_j phi_j and v = phi_i :
A_ij = integral_0^1 phi_i'(x) phi_j'(x) dx
b_i = integral_0^1 f(x) phi_i(x) dx
A x = b , x = (c_1, ..., c_N)
uniform mesh, spacing h, hat functions:
[ 2 -1 ] tridiagonal,
A = 1 [ -1 2 -1 ] symmetric,
--- [ -1 2 -1 ] positive definite
h [ . . .] (sparse: |i-j|<=1)Assembly, and the matrix's good manners
In practice nobody computes A_ij by one big global integral. Instead we assemble the matrix element by element. Loop over each element, compute a tiny local 2-by-2 (in 1D) or 3-by-3 (for triangles) element matrix of overlap integrals among only the shape functions living on that element, then add those small numbers into the right slots of the global A. Each element contributes a little, and the contributions pile up where elements share a node. This element-local loop is the engine of every real FEM code: it handles unstructured meshes and curved elements uniformly, and it is where the isoparametric mapping (warping a reference triangle onto a curved one) does its work.
The stiffness matrix from a symmetric weak form is symmetric (A_ij = A_ji, since the integrand phi_i'*phi_j' does not care about order) and positive definite (it descends from an energy that is genuinely minimized). That is a gift: a symmetric positive-definite system is exactly the comfortable case from the linear-algebra rungs — solvable by a Cholesky factorization for a direct solve, or by the conjugate gradient method for a large sparse one, with no pivoting drama. Two cautions stay honest, though. Refining the mesh makes A larger and worse conditioned — its condition number grows like O(1/h^2), so a finer mesh quietly eats some of your double-precision digits and slows iterative solvers. And the bare A is singular until the boundary conditions are imposed, the subject of the next paragraph.
Boundary conditions, the mass matrix, and where we land
Not all boundary conditions enter the same way — a genuinely subtle point worth pinning down. A condition fixing the value of u, like u(0) = 0, is essential: it cannot be expressed by the weak-form integral, so we enforce it directly by removing that node's unknown (it is not free, its coefficient is just 0). A condition on the flux, like u'(1) = g, is natural: it falls out of the integration-by-parts that built the weak form and simply appears in the right-hand side b — you do nothing special, it takes care of itself. This essential-versus-natural split is a recurring surprise for newcomers, and getting it backwards is a classic FEM bug.
One sibling matrix deserves a name. If the problem also carries an undifferentiated u term — a time-dependent heat equation u_t = u_xx, or an eigenvalue problem — the same assembly produces a second matrix M with entries M_ij = integral of phi_i(x)*phi_j(x) dx (the shape functions themselves, not their derivatives). This is the mass matrix, equally sparse and symmetric positive definite, and it is what couples FEM in space to the time-stepping methods of the ODE rungs: discretize space with elements and you are left with a system M*c'(t) = -A*c(t) of ordinary differential equations in the nodal coefficients — the method of lines again, now with a tasteful mass matrix on the time derivative.
Step back and see the bridge you have crossed. Last guide's abstract weak form has become a concrete, sparse, symmetric positive-definite system A x = b whose unknowns are the node values of a piecewise-linear solution. Every number in it is still an approximation: we restricted the search to tent functions (a modelling choice the error theory, and Cea's lemma, will tell us costs O(h^2) in the right norm), the integrals carry quadrature error, and the conditioning O(1/h^2) erodes a few digits as we refine. The next guide actually solves this system and reads its answer; the one after weighs its accuracy and refines the mesh where it hurts. The hard conceptual leap — turning calculus into a matrix — is now behind you.