The Finite Element Method as an Implementation of the Variational Approach

Finite Element Method

Motivation: Numerical Implementation of Weak Solutions

The Lax–Milgram theorem guarantees the existence and uniqueness of a weak solution. The Galerkin method approximates the infinite-dimensional problem by a finite-dimensional one: a finite-dimensional subspace $V_h \subset V$ is chosen. The finite element method (FEM) is a concrete implementation: the domain is partitioned into elements (triangles/tetrahedra) and a piecewise-polynomial basis is taken. This is the standard in engineering CAE computations.

Galerkin Method

Galerkin problem: Find $u_h \in V_h$ such that $a(u_h, v_h) = F(v_h)$ for all $v_h \in V_h$.

Basis ${\varphi_1, \ldots, \varphi_n}$. Decomposition $u_h = \sum_j u_j \varphi_j$. Substitution yields the linear system $Ku = f$:

• $K_{ij} = a(\varphi_j, \varphi_i)$ — stiffness matrix.
• $f_i = F(\varphi_i)$ — load vector.

Céa's Theorem: $|u - u_h|V \leq (M/\alpha) \cdot \inf{v_h \in V_h} |u - v_h|_V$. The Galerkin approximation is optimal up to the constant $M/\alpha$.

FEM P1 error estimate: $|u - u_h|{H^1} \leq C h |u|{H^2}$, where $h$ is the element size. The finer the mesh, the more accurate the solution.

Finite Elements

P1 elements (piecewise linear): $\varphi_i(x_j) = \delta_{ij}$ (“hat” function). On each triangle — a linear function. $\nabla \varphi_i$ is a piecewise-constant vector.

Stiffness matrix: $K_{ij} = \int_\Omega \nabla\varphi_i \cdot \nabla\varphi_j, dx$. Sparse: $K_{ij} \neq 0$ only for neighboring nodes. The number of nonzero elements is $O(N)$ with $N = \dim(V_h) \to$ efficient iterative solvers.

Numerical Example

Problem: FEM with $N=4$ equal elements for $-u'' = 1$ on $[0,1]$, $u(0) = u(1) = 0$.

Step 1. Nodes: $x_0=0$, $x_1=0.25$, $x_2=0.5$, $x_3=0.75$, $x_4=1$. Step size $h=0.25$. Interior nodes: $1,2,3$. Basis: $\varphi_1, \varphi_2, \varphi_3$ — “hats” peaking at $x_i$.

Step 2. Stiffness matrix (1D P1): $K_{ii} = 2/h = 8$, $K_{i,i+1} = -1/h = -4$. $K = \begin{bmatrix} 8 & -4 & 0 \ -4 & 8 & -4 \ 0 & -4 & 8 \end{bmatrix}$.

Step 3. Load vector: $f_i = \int_0^1 1 \cdot \varphi_i, dx = h = 0.25$. $f = [0.25, 0.25, 0.25]^T$.

Step 4. Solve $Ku = f$. By symmetry $u_1 = u_3$. From the 3rd row: $-4u_2 + 8u_3 = 0.25$. From the 1st: $8u_1 - 4u_2 = 0.25$. From the 2nd: $-4u_1 + 8u_2 - 4u_3 = 0.25 \to -8u_1 + 8u_2 = 0.25$ (using $u_1 = u_3$) $\to u_2 = u_1 + 1/32$.

From the 1st: $8u_1 - 4(u_1 + 1/32) = 0.25 \to 4u_1 = 0.25 + 0.125 = 0.375 \to$ $u_1 = 3/32 = 0.09375$. $u_2 = 4/32 = 0.125$. $u_3 = 3/32 = 0.09375$.

Step 5. Exact solution: $u(x) = x(1-x)/2$. $u(0.25) = 0.25 \cdot 0.75 / 2 = \mathbf{0.09375}$ ✓. $u(0.5) = 0.5 \cdot 0.5 / 2 = \mathbf{0.125}$ ✓. FEM P1 yields the exact nodal value for a quadratic exact solution (superconvergence)!

Step 6. As $h \to 0$: $|u - u_h|_{H^1} \to 0$ as $O(h)$. For $u \in H^2$ error in the $C[0,1]$ norm $\leq C h^2$. Just $N=4$ elements suffice for machine accuracy in this problem.

Real-World Application

Aerospace and automotive engineering: ANSYS, COMSOL, Abaqus are CAE systems based on FEM. In designing the aerodynamic components of the Boeing 787 airplane, thousands of FEM simulations of the strength and aerodynamics of carbon composite structures were performed.

Additional Aspects

The finite element method (FEM) is a numerical implementation of the weak formulation. The domain $\Omega$ is divided into simple elements (triangles in 2D, tetrahedra in 3D); the solution is sought in a finite-dimensional subspace $V_h \subset H^1_0$ consisting of piecewise-polynomial functions. Substitution into the variational form yields a linear system $K \cdot u = f$, where $K$ is a sparse stiffness matrix. Céa's lemma guarantees that the error $|u-u_h|$ is estimated by the distance from the exact solution to $V_h$, which provides quasi-optimal estimates. FEM is the standard in engineering computations: COMSOL, ANSYS, Abaqus, FEniCS use exactly this mathematics for structural analysis, heat transfer, electromagnetism, hydrodynamics, and biomechanics.

Connection with Other Areas of Mathematics

The finite element method is naturally embedded in the general theory of elliptic boundary value problems. On one hand, it relies on functional analysis: Riesz representation theorems for linear functionals, the theory of Sobolev spaces by Ladyzhenskaya, Uraltseva, Agmon–Douglis–Nirenberg. On the other hand, each finite element approximation leads to a large sparse linear algebraic system, for which Krylov subspace methods are used (Hestenes–Stiefel conjugate gradients, Saad–Schultz GMRES), and the stability of preconditioners is tied to the theory of spectral decompositions.

Topological aspects emerge through finite elements on de Rham forms. The works of Georges de Rham, then Nédélec and Raviart, showed how to construct elements compatible with the gradient–curl–divergence complex. Modern theory of homology-compatible finite elements (finite element exterior calculus by Arnold, Falk, Winther) is based on de Rham's theorem and commutative diagrams between continuous and discrete complexes.

Connections with probabilistic methods are expressed via stochastic coefficients and random loads. In stochastic FEM, ideas of polynomial chaos (Wiener, Ghanem) and generalized Hermite, Legendre polynomials are used; weak convergence here intertwines with the theory of stochastic processes and the law of large numbers for empirical measures. For nonlinear problems (plasticity, contact), variational inequalities and Brezis' theory of monotone operators come to the fore.

Numerical methods for differential equations form a unified “ecosystem”: comparison of FEM and finite difference methods relies on the Eckmann–Babuška theorem on stability and norm equivalence; mixed finite elements (Raviart–Thomas, Brezzi–Fortin–Marini) are closely tied to the theory of saddle-point problems and the Ladyzhenskaya–Babuška–Brezzi condition.

Historical Note and Evolution of the Idea

The roots of FEM go back to the works of Rayleigh and Ritz in the late 19th and early 20th centuries, where the variational principle of minimum potential energy was used for vibration and elasticity problems. In 1943, R. Courant, in an article in the Quarterly of Applied Mathematics, proposed piecewise linear functions on triangles for approximating solutions of the Laplace equation, effectively describing a prototype of finite elements. The term "finite element method" became established in the 1950s–1960s in aerospace industry works: A. Clough, J. Argyris, R. Klan. The classic monograph is O. Zienkiewicz's 1967 “The Finite Element Method in Structural and Continuum Mechanics”, which systematized engineering applications. In parallel, the rigorous mathematical foundation was developed by O. Axel, I. Babuška, F. Brézis, later the monographs by C. Ciarlet, P. G. Ciarlet, R. Rannacher. In the 1970s–1980s, tools for a posteriori estimates and adaptive meshes were formulated (Nödelman, Verfürth), and mixed and hybrid schemes emerged. The works of Babuška and Rogaš provided a general analysis of convergence and stability. In the late 20th century, emphasis shifted toward three-dimensional unstructured meshes, automatic generation of tetrahedral meshes (Delaunay, Bowyer–Watson algorithms), and parallel implementation. In the 21st century, isogeometric elements (Hughes et al., 2005) have been vigorously developed, combining CAD representation and FEM, as well as high-order methods and discrete differential geometry.