P1 and P2 are ready to be discretized, which leads to a common sub-problem (3). The basic idea is to replace the infinite-dimensional linear problem: :Find u \in H_0^1 such that :\forall v \in H_0^1, \; -\phi(u,v)=\int fv with a finite-dimensional version: where V is a finite-dimensional
subspace of H_0^1. There are many possible choices for V (one possibility leads to the
spectral method). However, we take V as a space of piecewise polynomial functions for the finite element method.
For problem P1 We take the interval (0,1), choose n values of x with 0=x_0 and we define V by: V = \{v:[0,1] \to \mathbb R\;: v\text{ is continuous, } v|_{[x_k,x_{k+1}]} \text{ is linear for } k = 0,\dots,n \text{, and } v(0)=v(1)=0 \} where we define x_0=0 and x_{n+1}=1. Observe that functions in V are not differentiable according to the elementary definition of calculus. Indeed, if v \in V then the derivative is typically not defined at any x = x_k, k = 1,\ldots,n. However, the derivative exists at every other value of x, and one can use this derivative for
integration by parts.
For problem P2 We need V to be a set of functions of \Omega. In the figure on the right, we have illustrated a
triangulation of a 15-sided
polygonal region \Omega in the plane (below), and a
piecewise linear function (above, in color) of this polygon which is linear on each triangle of the triangulation; the space V would consist of functions that are linear on each triangle of the chosen triangulation. One hopes that as the underlying triangular mesh becomes finer and finer, the solution of the discrete problem (3) will, in some sense, converge to the solution of the original boundary value problem P2. To measure this mesh fineness, the triangulation is indexed by a real-valued parameter h > 0 which one takes to be very small. This parameter will be related to the largest or average triangle size in the triangulation. As we refine the triangulation, the space of piecewise linear functions V must also change with h. For this reason, one often reads V_h instead of V in the literature. Since we do not perform such an analysis, we will not use this notation.
Choosing a basis To complete the discretization, we must select a
basis of V. In the one-dimensional case, for each control point x_k we will choose the piecewise linear function v_k in V whose value is 1 at x_k and zero at every x_j,\;j \neq k, i.e., v_{k}(x) = \begin{cases} {x-x_{k-1} \over x_k\,-x_{k-1}} & \text{ if } x \in [x_{k-1},x_k], \\ {x_{k+1}\,-x \over x_{k+1}\,-x_k} & \text{ if } x \in [x_k,x_{k+1}], \\ 0 & \text{ otherwise}, \end{cases} for k = 1,\dots,n; this basis is a shifted and scaled
tent function. For the two-dimensional case, we choose again one basis function v_k per vertex x_k of the triangulation of the planar region \Omega. The function v_k is the unique function of V whose value is 1 at x_k and zero at every x_j,\;j \neq k. Depending on the author, the word "element" in the "finite element method" refers to the domain's triangles, the piecewise linear basis function, or both. So, for instance, an author interested in curved domains might replace the triangles with curved primitives and so might describe the elements as being curvilinear. On the other hand, some authors replace "piecewise linear" with "piecewise quadratic" or even "piecewise polynomial". The author might then say "higher order element" instead of "higher degree polynomial". The finite element method is not restricted to triangles (tetrahedra in 3-d or higher-order simplexes in multidimensional spaces). Still, it can be defined on quadrilateral subdomains (hexahedra, prisms, or pyramids in 3-d, and so on). Higher-order shapes (curvilinear elements) can be defined with polynomial and even non-polynomial shapes (e.g., ellipse or circle). Examples of methods that use higher degree piecewise polynomial basis functions are the
hp-FEM and
spectral FEM. More advanced implementations (adaptive finite element methods) utilize a method to assess the quality of the results (based on error estimation theory) and modify the mesh during the solution aiming to achieve an approximate solution within some bounds from the exact solution of the continuum problem. Mesh adaptivity may utilize various techniques; the most popular are: • moving nodes (r-adaptivity) • refining (and unrefined) elements (h-adaptivity) • changing order of base functions (p-adaptivity) • combinations of the above (
hp-adaptivity).
Small support of the basis L of the discretized linear system The primary advantage of this choice of basis is that the inner products \langle v_j,v_k \rangle = \int_0^1 v_j v_k\,dx and \phi(v_j,v_k)=\int_0^1 v_j' v_k'\,dx will be zero for almost all j,k. (The matrix containing \langle v_j, v_k \rangle in the (j,k) location is known as the
Gramian matrix.) In the one dimensional case, the
support of v_k is the interval [x_{k-1},x_{k+1}]. Hence, the integrands of \langle v_j, v_k \rangle and \phi(v_j,v_k) are identically zero whenever |j-k|>1. Similarly, in the planar case, if x_j and x_k do not share an edge of the triangulation, then the integrals \int_{\Omega} v_j v_k\,ds and \int_{\Omega} \nabla v_j \cdot \nabla v_k\,ds are both zero.
Matrix form of the problem If we write u(x) = \sum_{k=1}^n u_k v_k(x) and f(x) = \sum_{k=1}^n f_k v_k(x) then problem (3), taking v(x) = v_j(x) for j = 1, \dots, n, becomes {{NumBlk|:|-\sum_{k=1}^n u_k \phi (v_k,v_j) = \sum_{k=1}^n f_k \int v_k v_j dx for j = 1,\dots,n. |}} If we denote by \mathbf{u} and \mathbf{f} the column vectors (u_1,\dots,u_n)^t and (f_1,\dots,f_n)^t, and if we let L = (L_{ij}) and M = (M_{ij}) be matrices whose entries are L_{ij} = \phi (v_i,v_j) and M_{ij} = \int v_i v_j dx then we may rephrase (4) as {{NumBlk||-L \mathbf{u} = M \mathbf{f}.|}} It is not necessary to assume f(x) = \sum_{k=1}^n f_k v_k(x). For a general function f(x), problem (3) with v(x) = v_j(x) for j = 1, \dots, n becomes actually simpler, since no matrix M is used, {{NumBlk||-L \mathbf{u} = \mathbf{b},|}} where \mathbf{b} = (b_1, \dots, b_n)^t and b_j = \int f v_j dx for j = 1, \dots, n. As we have discussed before, most of the entries of L and M are zero because the basis functions v_k have small support. So we now have to solve a linear system in the unknown \mathbf{u} where most of the entries of the matrix L, which we need to invert, are zero. Such matrices are known as
sparse matrices, and there are efficient solvers for such problems (much more efficient than actually inverting the matrix.) In addition, L is symmetric and positive definite, so a technique such as the
conjugate gradient method is favored. For problems that are not too large, sparse
LU decompositions and
Cholesky decompositions still work well. For instance,
MATLAB's backslash operator (which uses sparse LU, sparse Cholesky, and other factorization methods) can be sufficient for meshes with a hundred thousand vertices. The matrix L is usually referred to as the
stiffness matrix, while the matrix M is dubbed the
mass matrix.
General form of the finite element method In general, the finite element method is characterized by the following process. • One chooses a grid for \Omega . In the preceding treatment, the grid consisted of triangles, but one can also use squares or curvilinear polygons. • Then, one chooses basis functions. We used piecewise linear basis functions in our discussion, but it is common to use piecewise polynomial basis functions. Separate consideration is the smoothness of the basis functions. For second-order
elliptic boundary value problems, piecewise polynomial basis function that is merely continuous suffice (i.e., the derivatives are discontinuous.) For higher-order partial differential equations, one must use smoother basis functions. For instance, for a fourth-order problem such as u_{xxxx} + u_{yyyy} = f, one may use piecewise quadratic basis functions that are
C^1. Another consideration is the relation of the finite-dimensional space V to its infinite-dimensional counterpart in the examples above H_0^1. A
conforming element method is one in which space V is a subspace of the element space for the continuous problem. The example above is such a method. If this condition is not satisfied, we obtain a
nonconforming element method, an example of which is the space of piecewise linear functions over the mesh, which are continuous at each edge midpoint. Since these functions are generally discontinuous along the edges, this finite-dimensional space is not a subspace of the original H_0^1. Typically, one has an algorithm for subdividing a given mesh. If the primary method for increasing precision is to subdivide the mesh, one has an
h-method (
h is customarily the diameter of the largest element in the mesh.) In this manner, if one shows that the error with a grid h is bounded above by Ch^p, for some C and p > 0, then one has an order
p method. Under specific hypotheses (for instance, if the domain is convex), a piecewise polynomial of order d method will have an error of order p = d+1. If instead of making
h smaller, one increases the degree of the polynomials used in the basis function, one has a
p-method. If one combines these two refinement types, one obtains an
hp-method (
hp-FEM). In the hp-FEM, the polynomial degrees can vary from element to element. High-order methods with large uniform
p are called spectral finite element methods (
SFEM). These are not to be confused with
spectral methods. For vector partial differential equations, the basis functions may take values in \mathbb{R}^n. ==Various types of finite element methods==