A Brief Summary of the Finite Element Method for Differential Equations

The finite element (FE) method is a numerical technique for computing approximate solutions to complex mathematical problems described by differential equations. The method was developed in the 1950s to solve complicated problems in engineering, notably in elasticity and structural mechanics modeling involving elliptic partial differential equations and complicated geometries. But nowadays the range of applications is quite extensive. In particular, the FE method has been successfully applied to many problems such as fluid – structure interaction, thermomechanical, thermochemi-cal, thermo-chemo-mechanical problems, biomechanics, biomedical engineering, pie-zoelectric, ferroelectric, electromagnetics, and many others. This chapter contains a summary of the FE method. Since the remaining chapters of this textbook are based on the FE method, we present it in this chapter as a method for approximating solutions of ordinary differential equations (ODEs) and partial differential equations (PDEs).


An overview of the finite element method
Differential equations arise in many disciplines such as engineering, mathematics, sciences, economics, and many other fields. Unfortunately solutions to differential equations can rarely be expressed by closed formulas and numerical methods are needed to approximate their solutions. There are many numerical methods for approximating the solution to differential equations including the finite difference (FD), finite element (FE), finite volume (FV), spectral, and discontinuous Galerkin (DG) methods. These methods are used when the mathematical equations are too complicated to be solved analytically.
The FE method has become the standard numerical scheme for approximating the solution to many mathematical problems; see [1][2][3][4][5][6][7][8][9] and the references therein just to mention a few. In simple words, the FE method is a numerical method to solve differential equations by discretizing the domain into a finite mesh. Numerically speaking, a set of differential equations are converted into a set of algebraic equations to be solved for unknown at the nodes of the mesh. The FE method originated from the need to solve complex elasticity and structural analysis problems in civil and aeronautical engineering. The first development can be traced back to the work by Hrennikoff in 1941 [10] and Courant in 1943 [11]. Although these pioneers used different perspectives in their FE approaches, they each identified the one common and essential characteristic: mesh discretization of a continuous domain into a set of discrete sub-domains, usually called elements. Another fundamental mathematical contribution to the FE method is represented by Gilbert Strang and George Fix [12]. Since then, the FE method has been generalized for the numerical modeling of physical systems in many engineering disciplines including electromagnetism, heat transfer, and fluid dynamics.
The advantages of this method can be summarized as follows: 1. Numerical efficiency: The discretization of the calculation domain with finite elements yields matrices that are in most cases sparse and symmetric. Therefore, the system matrix, which is obtained after spatial and time discretization, is sparse and symmetric too. Both the storage of the system matrix and the solution of the algebraic system of equations can be performed in a very efficient way.
2. Treatment of nonlinearities: The modeling of nonlinear material behavior is well established for the FE method (e.g., nonlinear curves, hysteresis).
3. Complex geometry: By the use of the FE method, any complex domain can be discretized by triangular elements in 2D and by tetrahedra elements in 3D. 4. Applicable to many field problems: The FE method is suited for structural analysis, heat transfer, electrical/magnetical analysis, fluid and acoustic analysis, multi-physics, etc.
COMSOL Multiphysics (known as FEMLAB before 2005) is a commercial FE software package designed to address a wide range of physical phenomena. It is widely used in science and industry for research and development. It excels at modeling almost any multi-physics problem by solving the governing set of PDEs via the FE method. This software package is able to solve one, two and threedimensional problems. It comes with a modern graphical user interface to set up simulation models and can be scripted from Matlab or via its native Java API.
In this chapter, we introduce the FE method for several one-dimensional and two-dimensional model problems. Although the FE method has been extensively used in the field of structural mechanics, it has been successfully applied to solve several other types of engineering problems, such as heat conduction, fluid dynamics, seepage flow, and electric and magnetic fields. These applications prompted mathematicians to use this technique for the solution of complicated problems. For illustration, we will use simple one-dimensional and two-dimensional model problems to introduce the FE method.

The FE method for first-order linear IVPs
We first present the FE method as an approximation technique for solving the following first-order initial-value problem (IVP) using piecewise linear polynomials In order to apply the FE method to solve this problem, we carry out the following process.

Derive a weak form (variational formulation). This can be done by multiplying the ODE in (1) by a test function
integrating from a to b, using integration by parts, and applying v a ð Þ ¼ 0, to get 2. Generate a triangulation (also called a mesh) of the computational domain a, b ½ . For a one-dimensional problem, a mesh is a set of points in the interval a, b ½ , say, The point x i is called a node or nodal point. The length of the interval (called an element) (called a mesh size that measures how fine the partition is). If the mesh is uniformly distributed, then 3. Define a finite dimensional space over the triangulation: Let the solution u be in the space V. For the model problem (1), the solution space is V ¼ C 1 a, b ½ . We wish to construct a finite dimensional space (subspace) V h ⊂ V based on the mesh. When the FE space is a subspace of the solution space, the method is called conforming. It is known that in this case, the FE solution converges to the true solution provided the FE space approximates the given space in some sense [3]. Different finite dimensional spaces will generate different FE solutions.
Define the FE space as the set of all continuous piecewise linear polynomials ½ . An example of such a function is shown in Figure 1.
We remark that any function v ∈ V h is uniquely determined by its nodal values v x i ð Þ.
4. Construct a set of basis functions based on the triangulation. Since V h has finite dimension, we can find one set of basis functions. A basis for V h , where ϕ j ∈ V h are linearly independent. Then is the space spanned by the basis The simplest finite dimensional space is the piecewise continuous linear function space defined over the triangulation.
There are infinite number of sets of basis functions. We should choose a set of basis functions that are simple, have compact (minimum) support (that is, zero almost everywhere except for a small region), and meet the regularity requirement, that is, they have to be continuous, and differentiable except at nodal points. The simplest ones are the so-called hat functions satisfying ϕ i x i ð Þ ¼ 1 and The analytic form is (see Figure 2) This type of FE method (with similar trial and test space) is sometimes called a Galerkin method, named after the famous Russian mathematician and engineer Galerkin.
Implementation: The FE solution is a linear combination of the basis functions. since u h b ð Þ ¼ c N . Note that using the hat functions, we have u h x 0 ð Þ ¼ 0 and Thus, we get the following linear system Finally, we solve the linear system for c 1 , … , c N . We note that for i ¼ 1, 2, … , N À 1, we have Since it depends on f , we cannot generally expect to calculate it exactly. However, we can approximate it using a quadrature rule. Using the Trapezoidal rule

The FE method for first-order nonlinear IVPs
Here, we extend the FE method for the nonlinear IVP using piecewise linear polynomials The FE method consists of finding Finally, we solve the nonlinear system for c 1 , c 2 , … , c N using e:g:, Newton's method for systems of nonlinear equations. The system can be written as

The FE method for two-point BVPs
Here, we shall study the derivation and implementation of the FE method for two-point boundary-value problems (BVPs). For easy presentation, we consider the following model problem: Find u ∈ C 2 a, b ½ such that where u : Ω ¼ a, b ½ !  is the sought solution, q x ð Þ ≥ 0 is a continuous function on a, b ½ , and f ∈ L 2 a, b ½ . Under these assumptions, (3) has a unique solution u ∈ C 2 a, b ½ . For general q x ð Þ, it is impossible to find an explicit form of the solution. Therefore, our goal is to obtain a numerical solution via the FE method.

Different mathematical formulations for the 1D model
The model problem (3) can be reformulated into three different forms: (D)-form: the original differential equation (3). (V)-form: the variational form or weak form: The corresponding FE method is often called the Galerkin method. In other words, a Galerkin FE method is a FE method obtained from the variational form.
(M)-form: the minimization form: corresponding FE method is often called the Ritz method. Under some assumptions, the three different forms are equivalent, that is, they have the same solution as will be explained in the following theorem.

Galerkin method of the problem
To solve (3) using the FE method, we carry out the process described below. Usually, a FE method is always derived from the weak or variational formulation of the problem at hand.
Weak formulation of the problem: The Galerkin FE method starts by rewriting (3) in an equivalent variational formulation. To this end, let us define the vector by a test function v ∈ H 1 0 , integrating from a to b, and using integration by parts, we get Hence, the weak (or variational) form of (3) reads: Find We want to find u ∈ H 1 0 that satisfies (4). We note that a solution u to (4) is less regular than the solution u (3). Indeed, (4) has only u 0 whereas (3) contains u 00 . Furthermore, we can easily verify the following: 1. If u is strong solution (i:e:, solves (3)) then u is also weak solution (i:e:, solves (4)).
2. Conversely, if u is a weak solution with u ∈ C 2 a, b ½ , it is also strong solution.
3. Existence and uniqueness of weak solutions is obtained by the Lax-Milgram Theorem.
4. We can consider solutions with lower regularity using the weak formulation.

FE method gives an approximation of the weak solution.
From now on, we use the notation v

The FE formulation:
The FE method is based on the variational form (4). We note that the space H 1 0 contains many functions and it is therefore just as hard to find a function u ∈ H 1 0 which satisfies the variational Eq. (4) as it is to solve the original problem (3). Next, we study in details a special Galerkin method called the FE method. Let a ¼ The FE method consists of choosing a basis for the subspace V h that satisfies the following properties 1. The matrix A must be sparse (e:g: traditional or banded matrix). In this case, iterative methods for solving linear systems can be adapted to obtain an efficient solution.
2. u h must converge to the solution u of the original problem as h ! 0.
It is natural to obtain an approximation u h to u as follows: We call u h the FE approximation of u. We say that (5) is the Galerkin approximation of (4) and the method used to find u h ∈ V h is called Galerkin method.
FE approximation using Lagrange  1 elements: The simplest finite dimensional space is the piecewise continuous linear function space defined over the triangulation It is easy to show that V 1 h,0 has a finite dimension even although there are infinite number of elements in V 1 h,0 . The approximation of the FE method is therefore to look for an approximation u h within a small (finite dimensional) subspace Let V 1 h,0 be the space of all continuous piecewise linear functions, which vanish at the end points a and b. There are many types of basis functions ϕ i f g NÀ1 i¼1 . The simplest ones are the so-called hat functions satisfying ϕ i x j À Á ¼ δ ij , where δ ij is the Kronecker symbol. Note especially that there is no need to construct hat functions ϕ 0 and ϕ N since any function of V 1 h,0 must vanish at the end points x 0 ¼ a and The explicit expressions for the hat function ϕ i x ð Þ and its derivative ϕ 0 i x ð Þ are given by The FE approximation of (4) thus reads: We call u h the FE approximation of u. We say that (6) is the Galerkin approximation of (4) and the method used to find u h ∈ V 1 h,0 is called Galerkin method. It can be shown that (6) is equivalent to the N À 1 equations Derivation of the discrete system: Since u h ∈ V 1 h,0 , we can express it as a linear combination of hat functions i:e:, where c j are real numbers to be determined. We note that the coefficients c j , j ¼ 1, 2, … , N À 1 are the N À 1 nodal values of u h to be determined. Note that the index is only from 1 to N À 1, because of the zero boundary conditions. We We can use either the weak/variational form (V), or the minimization form (M), to derive a linear system of equations for the coefficients c j .
Substituting (8) into (7) yields The problem (7) is now equivalent to the following: Find the real numbers c 1 , c 2 , … , c NÀ1 that satisfy the linear system (9).
We note that the linear system (9) is equivalent to the system in matrix-vector form where Þmatrix, the so-called stiffness matrix when q ¼ 0, with entries and b ∈  NÀ1 , the so-called load vector, has entries To obtain the approximate solution we need to solve the linear system for the unknown vector c. We note that

Ritz method of the problem
The Ritz method is one of the earliest FE methods. However, not every problem has a minimization form. The minimization form for the model problem (3) is As before, we look for an approximate solution of the form (8). If we plug this into the functional form above, we get which is a multi-variable function of c 1 , c 2 , … , c NÀ1 and can be written as Taking the partial derivatives directly with respect to c i , we get Exchange the order of integration and the summation, we get which is exactly the same linear system (9) obtained using the Galerkin method.

Computer implementation
It is straightforward to calculate the entriesâ By symmetry, we also havê So, we can easily verify that Thus, the matrix A ¼â i,j þã i,j À Á is tridiagonal and has the form Finally, we obtain the following system:

Remark 2.2 Suppose that the partition is uniform
Then the stiffness matrix A and the load vector b have the form:  : Finally, we obtain the following system: c 0 ¼ c N ¼ 0 and which is the same system obtained using the finite difference method, where u 00 is approximated using the second-order midpoint formula We conclude that the above FE method using the composite trapezoidal rule is equivalent to the finite difference method of order 2.

Existence, uniqueness, and basic a priori error estimate
Theorem 2.2 The linear system (10) obtained using the FE method has a unique solution. Consequently, the FE method solution u h is unique.
Next, we state a general convergence result for the Galerkin method. We first define the following norm and semi-norm: Let u be the solution to (4) and u h be the solution to (6). Then there exists a constant C such that where C is given by (14) gives This inequality suggest that the error between u and u h is controlled by the interpolation error u À πu in the Á j j 1 -norm.
Let u be the solution to (4) and u h be the solution to (6). Then there exists a constant C such that Remark 2.4 1. If the partition is not uniform then we obtain the same error estimate with 2. The error is expressed in terms of the exact solution u. If it is expressed in terms of the computed solution u h it is an a posteriori error estimate (this yields a computable error bound).
4. u h is the best approximation within the space V 1 h,0 with respect to the v 0 k knorm. 5. The norm v 0 k k is referred to as the energy norm and has often a physical meaning.

Boundary conditions
In problem (3) we considered a homogeneous Dirichlet boundary conditions. Here, we extend the FE method to boundary conditions of different types. There are three important types of boundary conditions (BCs): 1. Dirichlet BCs: u a ð Þ ¼ α and u b ð Þ ¼ β for two real numbers α and β. This BC is also known as strong BC or essential BC.
2. Neumann BCs: u 0 a ð Þ ¼ α and u 0 b ð Þ ¼ β for two real numbers α and β. This BC is also known as natural BCs.

Robin BCs
Note that any combination is possible at the two boundary points. Nonhomogeneous Dirichlet boundary conditions: Let us consider the following two-point BVP: find u ∈ C 2 a, b ð Þ such that where α and β are given constants and f ∈ C a, b ð Þis a given function. In this case, the admissible function space the FE space V 1 h,0 defined earlier remain the same. Multiplying (15) by a test function v ∈ H 1 0 and integrating by parts gives Hence, the weak or variational form of (15) reads: Given Let V 1 h and V 1 h,0 , respectively, be the space of all continuous piecewise linear functions and the space of all continuous piecewise linear functions which vanish at the endpoints a and b. We also let a ¼ ½ . Moreover let ϕ i f g be the set of hat basis functions of V h associated with the N þ 1 nodes x j , j ¼ 0, 1, … , N, such that ϕ i x j À Á ¼ δ ij . The FE approximation of (16) thus reads: It can be shown that (17) is equivalent to the N À 1 equations Expanding u h as a linear combination of hat functions where the coefficients c j , j ¼ 1, 2, … , N À 1 are the N À 1 nodal values of u h to be determined.
Substituting (19) into (18) yields Þsystem of equations for c j . In matrix form we write where A is a N À 1 ð ÞÂ N À 1 ð Þmatrix, the so-called stiffness matrix, with entries Þvector containing the unknown coefficients c j , j ¼ 1, 2, … , N À 1, and b is a N À 1 ð Þvector, the so-called load vector, with entries Computer Implementation: The explicit expression for a hat function ϕ i x ð Þ is given by For simplicity we assume the partition is uniform so that h i ¼ h for i ¼ 1, 2, … , N. Hence the derivative ϕ 0 i x ð Þ is either À 1 h , 1 h , or 0 depending on the interval. It is straightforward to calculate the entries of the stiffness matrix. For |i À j| > 1, we have a i,j ¼ 0, since ϕ i and ϕ j lack overlapping support. However, if i ¼ j, then where we have used that Changing i to i À 1 we also have Thus the stiffness matrix is : The entries b i of the load vector must often be evaluated using quadrature, since they involve the function f which can be hard to integrate analytically. For example, using the trapezoidal rule one obtains the approximate load vector entries Assembly: We rewrite (20), (21), (22) as We note that u h a ð Þ ¼ α ¼ u a ð Þ and u h b ð Þ ¼ β ¼ u b ð Þ. Therefore, we see that the system matrix A remains the same, and only the first and last entries of the load vector b need to be modified because of the definition of the basis functions ϕ 0 , … , ϕ N f g . An alternative approach is to use all the basis functions ϕ 0 , … , ϕ N f gto form a larger system of equation, i:e:, and N þ 1 ð ÞÂ N þ 1 ð Þsystem. The procedure for inserting the boundary conditions into the system equation is: enter zeros in the first and N þ 1 ð Þ-th rows of the system matrix A except for unity in the main diagonal positions of these two rows, and enter α and β in the first and N þ 1 ð Þ-th rows of the vector b, respectively.
General boundary conditions: Let us consider the following two-point BVP: find u ∈ C 2 a, b ð Þ such that where α, β and γ are given numbers and f ∈ C a, b ð Þ is a given function. The boundary condition at x ¼ b is called a Robin boundary condition (combination and u and u 0 is prescribed at x ¼ b). In this case, the admissible function space is modified to

Multiplying (23) by a function v ∈ H 1 0 and integrating by parts gives
Hence, the weak or variational form of (23) reads: Given u a ð Þ ¼ α, find the approximate solution u ∈ H 1 0 , such that The FE space V 1 h is now the set of all continuous piecewise linear functions which vanish at the end point a. The FE approximation of (24) thus reads: Find the piecewise linear approximation u h to the solution u satisfies with u h a ð Þ ¼ α. As before, (25) can be formulated in matrix form.

Model problem with coefficient and general Robin BCs
Let us consider the following two-point BVP: where p ¼ p x ð Þ with p x ð Þ ≥ p 0 > 0, f ∈ L 2 I ð Þ, κ 0 , κ 1 ≥ 0, and α, β are given numbers. Let Multiplying (26) by a function v ∈ V and integrating by parts gives We gather all u-independent terms on the left and obtain The FE method consists of finding Implementation: We need to assemble a stiffness matrix A and a load vector b.
Let for simplification p ¼ 1. Then the matrix A and the vector b (when using the trapezoidal rule) are given by denotes the vector space of polynomials in one variable and of degree less than or equal to k. The FE method for Lagrange P 2 elements involves the discrete space: These spaces are composed of continuous, piecewise parabolic functions (polynomials of degree less than or equal to 2). The P 2 FE method consists in applying the internal variational approximation approach to these spaces.

Lemma 2.2 The space V
h is uniquely defined by its values at the mesh vertices x j , j ¼ 0, 1, … , N and at the midpoints x jþ 1 2 ¼ is the basis of the shape functions ϕ j defined as: Figure 3 shows the global shape functions for the space V 2 h and the three quadratic Lagrange P 2 shape functions on the reference interval À1, 1 ½ .

Finite Element Methods and Their Applications
Remark 2.5 Notice that we have: 0,h is uniquely defined by its values at the mesh vertices x j , j ¼ 1, 2, … , N À 1 and at the midpoints x jþ 1 2 is the basis of the shape functions ϕ j defined as: with ϕ ξ ð Þ and ψ ξ ð Þ are defined by (28).

Homogeneous boundary conditions
The variational formulation of the internal approximation of the Dirichlet BVP (3) consists now in finding u h ∈ V 2 0,h , such that: Here, it is convenient to introduce the notation x j 2 , j ¼ 1, … , 2N À 1 for the mesh points and ϕ j 2 , j ¼ 1, … , 2N À 1 for the basis of V 2 0,h . Using these notations, we have: ≈ u x j 2 are the unknowns coefficients. This formulation leads to solve in  2NÀ1 a linear system: is the unknown vector containing the coef- and load vector b ∈  2NÀ1 has entries Since the shape functions ϕ i have a small support, the matrix A is mostly composed of zeros. However, the main difference with the Lagrange P 1 FE method, the matrix A is no longer a tridiagonal matrix.
Computer Implementation: The coefficients of the matrix A can be computed more easily by considering the following change of variables, for ξ ∈ À1, 1 ½ : Hence, the shape functions can be reduced to only three basic shape functions (Figure 3): Their respective derivatives are This approach consists in considering all computations on an interval on the reference interval À1, 1 ½ . Thus, we have: In this case, the elementary contributions of the element I i to the stiffness matrix and to the mass matrix are given by the 3 Â 3 matrices K I i and M I i :

5:
Coefficients of the right-hand side b: Usually, the function f is only known by its values at the mesh points xi 2 , i ¼ 0, 1, … , 2N and thus, we use the decomposition of f in the basis of shape functions ϕi Each component bi 2 of the right-hand side vector is obtained as bi Using the previous decomposition of f , we obtain: Thus, the problem is reduced to computing the integrals It is easy to see that we obtain expressions very similar to that of the mass matrix. More precisely, the element I i ¼ x iÀ1 , x i ½ will contribute to only three components of indices i À 1, i À 1 2 and i as:

Nonhomogeneous boundary conditions
Consider the following two-point BVP: find u ∈ C 2 a, b ð Þ such that where α and β are given constants and f ∈ C a, b ð Þ is a given function.

Multiplying (29) by a function
and integrating by parts gives Hence, the weak or variational form of (29) reads: Given u a ð Þ ¼ α, u b ð Þ ¼ β, Let V 2 h and V 2 h,0 , respectively, be the space of all continuous piecewise quadratic functions and the space of all continuous piecewise quadratic functions which vanish at the end points a and b, on a uniform partition a ¼ ½ . The FE method scheme consists of finding u h ∈ V 2 h , such that: Introduce the notation x j 2 , j ¼ 0, 1, … , 2N À 1, 2N for the mesh points and ϕ j 2 , j ¼ 0, 1, … , 2N À 1, 2N for the basis of V 2 h and ϕ j 2 , j ¼ 1, … , 2N À 1 for the basis of V 2 0,h . Using these notations, we have: are the unknowns coefficients. We note that c 0 ¼ This formulation leads to solve in  2NÀ1 a linear system: and the load vector b ∈  2NÀ1 has entries Clearly, the only extra terms are given in the vector with entries

Meshes
Let Ω ⊂  2 bounded with ∂Ω assumed to be polygonal. A triangulation T h of Ω is a set of triangles T such that Ω ¼ ⋃ T ∈ T h T, and two triangles intersect by either a common triangle edge, or a corner, or nothing. Corners will be referred to as nodes. We let h T ¼ diam T ð Þ the length or the largest edge. Let T h have N nodes and M triangles. The data is stored in two matrices. The matrix P ∈  2ÂN describes the nodes ( x 1 , y 1 À Á , … , x N , y N À Á Þ and the matrix K ∈  3ÂM describes the triangles, i:e:, it describes which nodes (numerated from 1 to N) form a triangle T and how it is orientated: This means that triangle T i is formed by the nodes n α i , n β i , and n γ i (enumeration in counter-clockwise direction).
The Delaunay algorithm determine a triangulation with the given points as triangle nodes. Delaunay triangulations are optimal in the sense that the angles of all triangles are maximal.
Matlab has a built in toolbox called PDE Toolbox and includes a mesh generation algorithm.

Piecewise polynomial spaces
Let T be a triangle with nodes Given v i we compute c i by Theorem 3.1 Suppose that u ∈ C 2 T ð Þ. Then the following hold where C is a generic constant independent of h T and u, but it depends on the ratio between smallest and largest interior angle of the triangle T. Now, we consider the piecewise continuous interpolant πu ¼ ð Þ for all T ∈ T h . Then the following hold where C is a generic constant independent of h and u, but it depends on the ratio between smallest and largest interior angle of the triangles of T h . Here

L 2 -projection
Let Ω ⊂  2 . We consider the space The problem of finding P h u ∈ V h is equivalent to solve the following linear system ð The problem can be expressed as a linear system of equations Mc ¼ b, where c ¼ c 1 , c 2 , … , c N ½ t and the entries of the matrix M ∈  NÂN and the vector b ∈  N are given by In general, we use a quadrature rule to approximate integrals. The general form is where the ω j 0 s denote the weights and the N j À Á 0 s the quadrature points.

Lemma 3.1
The mass matrix M with entries m ij ¼ Ð Ω ϕ i ϕ j dxdy is symmetric and positive definite. Theorem 3.3 For any u ∈ L 2 Ω ð Þ the L 2 -projection P h u exists and is unique.

A priori error estimate
Theorem 3.4 Let u ∈ L 2 Ω ð Þ and let P h u be the L 2 -projection of u, then

The FE method for general elliptic problem
The FE method was designed to approximate solutions to complicated equations of elasticity and structural mechanics, usually modeled by elliptic type equations, with complicated geometries. It has been developed for other applications as well.
Consider the following two-dimensional elliptic problem: Find u such that where a > 0, b ≥ 0, κ ≥ 0, f ∈ L 2 Ω ð Þ and g ∈ C 0 ∂Ω ð Þ. We seek a weak solution u in In order to derive the weak formulation, we multiply (31) with v ∈ V, integrate over Ω and use Green's formula to obtain ð We obtain the weak form:

Finite Element Methods and Their Applications
We can formulate the method as in the 1D case by using the weak formulation (32). The FE method in 2D is defined as follows: Find u h ∈ V h such that where |E| is the length of E and δ ij is 1 for i ¼ j and 0 else. Assembly of load vector: We use a corner quadrature rule for approximating the integral. We obtain for

The Dirichlet problem
Consider the following Dirichlet Problem: Find u such that where f ∈ L 2 Ω ð Þ and g ∈ C 0 ∂Ω ð Þ. We seek a weak solution u in So the weak problem reads: Find u ∈ V g such that Assume that g is piecewise linear on ∂Ω with respect to the triangulation. Then the FE method in 2D is defined as follows: Assume that we have N nodes and J boundary nodes, then the matrix form of the FE method problem reads: Note that c 1 ∈  J is known (it contains the values of g in the boundary nodes). We can therefore solve the simplified problem reading: find c 0 ∈  NÀJ with A 0,0 c 0 ¼ b 0 À A 0,g c 1 .

The Neumann problem
Consider the following Neumann Problem: Find u such that where f ∈ L 2 Ω ð Þ and g ∈ C 0 ∂Ω ð Þ. Let us try to seek a solution to this problem in

Multiplying (35) by a test function
v ∈ V, integrating over Ω, and using Green's formula, we get ð

vgds:
Thus, the variational formulation reads: find u ∈ V such that ð In order to guarantee solvability, we note that if v ¼ 1 then we have

gds:
Therefore we need to assume the following compatibility condition ð Ω fdxdy þ to ensure that a solution can exist. Note that if u exists, it is only determined up to a constant, since u þ c is a solution if u is a solution and c ∈ . To fix this constant and obtain a unique solution a common trick is to impose the additional constraint Ð Ω udxdy ¼ 0. We therefore define the weak solution spacê which contains only functions with a zero mean value. This is a called a quotient space. This space guarantees a unique weak solution (with weak formulation as usual with test functions in V). So the weak problem reads: Find u ∈V such that ð Now, the FE method takes the form: find u h ∈V h ⊂V such that whereV h is the space of all continuous piecewise linear functions with a zero mean.

Finite elements for mixed Dirichlet-Neumann conditions
Here we describe briefly how Neumann conditions are handled in twodimensional finite elements. Suppose Ω is a domain in either  2 or  3 and assume that ∂Ω has been partitioned into two disjoint sets: ∂Ω ¼ Γ 1 ∪ Γ 2 . We consider the following BVP: where f ∈ L 2 Ω ð Þ. As for the 1-D case, Dirichlet conditions are termed essential boundary conditions because they must be explicitly imposed in the FE method, while Neumann conditions are called natural and need not be mentioned. We therefore define the space of test functions bŷ Multiplying (36) by a test function v ∈V and integrating over Ω, we get ð since v ¼ 0 on Γ 1 and ∇u Á n on Γ 1 . Thus the weak form of (36) is: We now restrict our discussion once more to two-dimensional polygonal domains. To apply the FE method, we must choose an approximating subspace of V. Since the boundary conditions are mixed, there are at least two points where the boundary conditions change from Dirichlet to Neumann. We will make the assumption that the mesh is chosen so that all such points are nodes (and that all such nodes belong to Γ 1 , that is, that Γ 1 includes its "endpoints"). We can then choose the approximating subspace ofV as follows: A basis for V h is formed by including all basis functions corresponding to interior boundary nodes that do not belong to Γ 1 . If the BVP includes only Neumann conditions, then the stiffness matrix will be singular, reflecting the fact that BVP either does not have a solution or has infinitely many solutions. Special care must be taken to compute a meaningful solution to the resulting linear system.

Inhomogeneous Dirichlet conditions on a rectangle
In a two-dimensional problem, inhomogeneous boundary conditions are handled just as in one dimension. Inhomogeneous Dirichlet conditions are addressed via the method of shifting the data (with a specially chosen piecewise linear function), while inhomogeneous Neumann conditions are taken into account directly when deriving the weak form. Both types of boundary conditions lead to a change in the load vector.
We have thus replaced each g i by a function h i which differs from g i by a linear function, and which has value zero at the two endpoints: The reader should notice how the second term interpolates between the boundary values on Γ 1 and Γ 3 , while the third term interpolates between the boundary values on Γ 2 and Γ 4 . In order for these two terms not to interfere with each other, it is necessary that the boundary data be zero at the corners. It was for this reason that we transformed the g i 0 s into the h i 0 s. The first term in the formula for w undoes this transformation. It is straightforward to verify that w satisfies the desired boundary conditions.

Inhomogeneous Neumann conditions on a rectangle
We can also apply the technique of shifting the data to transform a BVP with inhomogeneous Neumann conditions to a related BVP with homogeneous Neumann conditions. However, the details are somewhat more involved than in the Dirichlet case. Consider the following BVP with the Neumann conditions where Γ 1 , Γ 2 , Γ 3 , and Γ 4 are, respectively, the bottom, right, top, and left boundary edges of the rectangular domain Ω ¼ 0, a ð ÞÂ 0, b ð Þ. We first note that this is equivalent to We make the following observation: If there is a twice-continuously differentiable function u satisfying the given Neumann conditions, then, since u xy ¼ u yx , we have Àu xy x, 0 ð Þ ¼ g 0 1 x ð Þ, À u yx 0, y ð Þ ¼ g 4 0 y ð Þ, which together imply that g 0 1 0 ð Þ ¼ g 4 0 0 ð Þ. By similar reasoning, we have all of the following conditions: g 0 1 0 ð Þ ¼ g 4 0 0 ð Þ, g 0 1 0 ð Þ ¼ g 4 0 0 ð Þ, À g 0 1 a ð Þ ¼ g 0 2 0 ð Þ, g 0 2 b ð Þ ¼ g 0 3 a ð Þ: We will assume that (41) holds.
We obtain the weak form: Find u ∈ V such that ð Ω ∇u Á ∇vdxdy ¼ λ The FE method in 2D is defined as follows: Find λ h ∈  and u h ∈ V h such that ð where . Implementation: Substituting u h ¼ P N j¼1 c j ϕ j into (45) and picking v h ¼ ϕ i , we obtain This leads to an algebraic system of the form Ac ¼ λ h Mc, i:e: an algebraic eigenvalue problem.

Error analysis
Consider the following model Problem: Find u such that The weak form: Find u ∈ V 0 such that ð The FE approximation is defined as follows: Find u h ∈ V h,0 such that ð Now, let v j j j j j j 2 ¼ Ð Ω ∇v Á ∇vdxdy ¼ Ð Ω ∇v j j 2 dxdy be the energy norm on V 0 . There are two different kinds of error estimates, a priori estimates, where the error is bounded in terms of the exact solution, and a posteriori error estimates, where the error is bounded in terms of the computed solution.
Theorem 3.8 (A priori error bound) Let u ∈ V 0 denote the weak solution and u h ∈ V h,0 the corresponding FE method approximation. Then Theorem 3.9 Let u ∈ V 0 denote the weak solution and u h ∈ V h,0 the corresponding FE method approximation. If u ∈ C 2 Ω ð Þ, then there exists C independent of h T and u such that

The FE method for elliptic problems with a convection term
Consider the following convection-diffusion problem: Find u such that We seek a weak solution u in In order to derive the weak formulation, we multiply (46) with v ∈ V 0 , integrate over Ω and use Green's formula to obtain ð Ω fvdxdy ¼ À

cuvdxdy:
Note that there is no need to apply Green's formula to Ð Ω vb Á ∇udxdy. We obtain the weak form: Find u ∈ V 0 such that ð The FE method in 2D is defined as follows: where ð Þ :

Conclusion
In this chapter, we introduced the finite element (FE) method for approximation the solutions to ODEs and PDEs. More specifically, the FE method is presented for first-order initial-value problems for OEDs, second-order boundary-value problems for ODEs, second-order elliptic PDEs, second-order heat and wave equations. The remaining chapters of this textbook are based on the FE method. The derivation of the FE method for other problems is straightforward. In the remaining chapters, the FE method will developed to solve complicated problems in engineering, notably in elasticity and structural mechanics modeling involving elliptic partial differential equations and complicated geometries. For more details, we refer the reader to [1][2][3][4][6][7][8][9] and the references therein.

Author details
Mahboub Baccouch Department of Mathematics, University of Nebraska at Omaha, Omaha, NE, USA *Address all correspondence to: mbaccouch@unomaha.edu © 2021 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.