Open access peer-reviewed chapter

A Brief Summary of the Finite Element Method for Differential Equations

Written By

Mahboub Baccouch

Submitted: June 9th, 2020 Reviewed: December 8th, 2020 Published: February 17th, 2021

DOI: 10.5772/intechopen.95423

From the Edited Volume

Finite Element Methods and Their Applications

Edited by Mahboub Baccouch

Chapter metrics overview

View Full Metrics

Abstract

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, thermochemical, thermo-chemo-mechanical problems, biomechanics, biomedical engineering, piezoelectric, 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).

Keywords

• the finite element method
• initial-value problems
• boundary-value problems
• Laplace equation
• heat equation
• wave equation

1. Introduction

1.1 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 three-dimensional 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.

2. The FE method for ODEs

2.1 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

u=fx,xab,ua=u0.E1

In order to apply the FE method to solve this problem, we carry out the following process.

1. Derive a weak form (variational formulation). This can be done by multiplying the ODE in (1) by a test function vxV0=vL2ab:v2+v2<,va=0, where v2=abv2xdx, integrating from a to b, using integration by parts, and applying va=0, to get

abfvdx=abuvdx=abuvdx+ubvbuava=abuvdx+ubvb.

2. Generate a triangulation (also called a mesh) of the computational domain ab. For a one-dimensional problem, a mesh is a set of points in the interval ab, say, a=x0x1xN=b. The point xi is called a node or nodal point. The length of the interval (called an element) Ii=xi1xi is hi=xixi1. Let h=max1iNhi (called a mesh size that measures how fine the partition is). If the mesh is uniformly distributed, then xi=a+ih, i=0,1,,N, where h=baN.

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=C1ab. We wish to construct a finite dimensional space (subspace) VhV 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 Vh={v:vIiP1Ii,i=1,2,,N,va=0}, where P1Ii is the space of polynomials of degree 1 on Ii. Functions in Vh are linear on each Ii, and continuous on the whole interval ab. An example of such a function is shown in Figure 1.

We remark that any function vVh is uniquely determined by its nodal values vxi.

4. Construct a set of basis functions based on the triangulation. Since Vh has finite dimension, we can find one set of basis functions. A basis for Vh is ϕjj=0N, where ϕjVh are linearly independent. Then Vh=vhxVvhx=j=0Ncjϕjx is the space spanned by the basis functions ϕii=0N. The simplest finite dimensional space is the piecewise continuous linear function space defined over the triangulation.

Vh=vhxVvhxis piecewise continuous linear overabwithvha=0.

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 ϕixi=1 and ϕixj=0 for ij. The analytic form is (see Figure 2)

ϕ0x=x1xh,xI1,0,else,,ϕNx=xxN1h,xIN,0,else,,ϕix=xxi1h,xIi,xi+1xh,xIi+1,0,else.

1. Approximate the exact solution u by a continuous piecewise linear function uhx. The FE method consists of finding uhVh such that

abuhvdx+uhbvb=abfvdx,vVh.

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. Writing uhx=j=0Ncjϕjx, where c0,c1,,cN are unknowns, and choosing v=ϕi,i=1,2,,N to get

j=0Ncjabϕjϕidx+cNϕib=abfϕidx,i=1,2,,N,

since uhb=cN. Note that using the hat functions, we have uhx0=0 and uhxi=j=0Ncjϕjxi=ciϕixi=ci for i=1,2,,N. Thus, we get the following linear system

j=1Ncjabϕjϕidx+cNϕib=abϕ0ϕidx,i=1,2,,N.

Finally, we solve the linear system for c1,,cN. We note that for i=1,2,,N1, we have

abϕiϕidx=xi1xi+1ϕiϕidx=1hixi1xixxi1hidx1hixixi+1xi+1xhidx=0.

However, for i=N, we have

abϕNϕN'dx=xN1xNϕNϕN'dx=xN1xNxxN1hNxxN1hNdx=1hNxN1xNxxN1hNdx=12.

Similarly, for i=1,2,,N, we have

abϕi1ϕidx=xi1xiϕi1ϕidx=xi1xixixhixxi1hidx=1hixi1xixixhidx=12,
abϕi+1ϕidx=xixi+1ϕi+1ϕidx=xixi+1xxihi+1xi+1xhi+1dx=1hi+1xixi+1xxihi+1dx=12.

We next calculate abfϕidx. 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 abfxdxba2fa+fb and using ϕixi1=ϕixi+1=0 and ϕixi=1, we get

abfϕidx=xi1xifϕidx+xixi+1fϕidxhi+hi+12fxi,i=1,2,,N1,
abfxϕNdx=xN1xNfxϕNdxhN2fxN1ϕNxN1+fxNϕNxN=hN2fxN.

Thus, we obtain the following linear system of equations

01200120120012012001212c1c2cN1cN=h1+h22fx1h2+h32fx2hN1+hN2fxN1hN2fxN.

The determinant of the above matrix is 12N. Thus, the system has a unique solution c1,c2,,cN.

Remark 2.1Suppose thatua=u0, then we letuhx=j=0Ncjϕjx. Sinceu0=uhx0=j=1Ncjϕjx0=c0ϕ0x0=c0, we only need to findc1,c2,,cN. Choosingv=ϕi,i=1,2,,N, we get the following linear system

j=1Ncjabϕjϕidx+cNϕib=abfϕidx+u0abϕ0ϕidx,i=1,2,,N.

Finally, we solve the linear system for c1,,cN. We note that abϕ0ϕidx=0 for i=2,,N and

abϕ0ϕ1'dx=x0x1x1xh1xx0h1dx=1h1x0x1x1xh1dx=12.

Following the same steps used for the case ua=0, we obtain the following linear system of equations

01200120120012012001212c1c2cN1cN=h1+h22fx1+u02h2+h32fx2hN1+hN2fxN1hN2fxN.

2.2 The FE method for first-order nonlinear IVPs

Here, we extend the FE method for the nonlinear IVP using piecewise linear polynomials

u=fxu,xab,ua=u0.E2

The FE method consists of finding uhVh={v:vIiP1Ii,i=1,2,,N,va=0}, such that

uhbvbabuhvdx=abfxuhvdx,vVh.

Writing uhx=j=0Ncjϕjx and choosing v=ϕi,i=1,2,,N, we get

cNϕiabϕNϕidxj=0N1cjabϕjϕidxabfxj=0Ncjϕjϕidx=0,i=1,2,,N,

where uhx0=c0=u0. Finally, we solve the nonlinear system for c1,c2,,cN using e.g., Newton’s method for systems of nonlinear equations. The system can be written as Fic1c2cN=0,i=1,2,,N, where

Fi=cNϕiabϕNϕidxj=0N1cjabϕjϕidxabfxj=0Ncjϕjϕidx,i=1,2,,N.

Let αi=j=0Ncjabϕjϕidx and βi=abfxj=0Ncjϕjϕidx. Then, for i=1,2,,N1,

αi=ci1xi1xiϕi1ϕidx+cixi1xiϕiϕidx+xixi+1ϕiϕidx+ci+1xixi+1ϕi+1ϕidx=ci1xi1xixixhi2dx+cixi1xixxi1hi2dxxixi+1xi+1xhi+12dxci+1xixi+1xxihi+12dx=12ci1+ci121212ci+1=12ci112ci+1,αN=cN1xN1xNϕN1ϕN'dx+cNxN1xNϕNϕN'dx=cN1xN1xNxNxhN2dx+cNxN1xNxxN1hN2dx=12cN1+12cN.

Similarly,

βi=xi1xi+1fxj=0Ncjϕjϕidx=xi1xifxj=0Ncjϕjϕidx+xixi+1fxj=0Ncjϕjϕidx.

Using Simpson’s Rule abfxdxba6fa+4fa+b2+fb, and using ϕixi1=ϕixi+1=0, ϕixi=1, j=0Ncjϕjxi1+hi2=ci1+ci2, ϕixi1+hi2=12, j=0Ncjϕjxi=ci, we have, for i=1,2,,N1,

βihi3fxi1+hi2ci1+ci2+hi+hi+16fxici+hi+13fxi+hi+12ci+ci+12.

However, for i=N, we have

βNhN62fxN1+h2cN1+cN2+fxNcN.

Next, we compute the Jacobian matrix with entries

Ji,j=Ficj=abϕjϕidxabfuxj=0Ncjϕjϕjϕidx=ai,jbi,j,i=1,2,,N.

We already computed the entries ai,j as

ai,i1=abϕi1ϕidx=12,ai,i=abϕiϕidx=0,i=1,2,,N1,aN,N=abϕNϕN'dx=12,ai,i+1=abϕi+1ϕidx=12.

Using Simpson’s Rule, we get

bi,i1=xi1xiϕi1ϕifuxj=0Ncjϕjdxhi6fuxi1+hi2ci1+ci2,bi,i+1=xixi+1ϕi+1ϕifuxj=0Ncjϕjdxhi+16fuxi+hi+12ci+ci+12,bi,i=xi1xiϕi2fuxj=0Ncjϕjdx+xixi+1ϕi2fuxj=0Ncjϕjdxhi6fuxi1+hi2ci1+ci2+hi+hi+16fuxici+hi+16fuxi+hi+12ci+ci+12,bN,N=xN1xNϕN2fuxj=0NcjϕjdxhN6fuxN1+h2cN1+cN2+fuxNcN.

2.3 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 uC2ab such that

u+qxu=fx,xΩ=ab,ua=ub=0,E3

where u:Ω¯=abR is the sought solution, qx0 is a continuous function on ab, and fL2ab. Under these assumptions, (3) has a unique solution uC2ab. For general qx, it is impossible to find an explicit form of the solution. Therefore, our goal is to obtain a numerical solution via the FE method.

2.3.1 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: abuvdx+abquvdx=abfvdx, for any test function v in the Sobolev space H01ab=vL2ab:v2+v2<va=vb=0, where v2=abv2xdx. 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: minvxH01abab12v2+12qv2fvdx. The 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.

Theorem 2.1 (Mathematical equivalences)Suppose thatuexists and continuous onab. Then we have the following mathematical equivalences.

(D) is equivalent to (V), (V) is equivalent to (M), and (M) is equivalent to (D).

2.3.2 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 space H01=vL2ab:v2+v2<va=vb=0. Multiplying (3) by a test function vH01, integrating from a to b, and using integration by parts, we get

abfvdx=abuvdx+abquvdx=abuvdx+abquvdx,

since va=vb=0. Hence, the weak (or variational) form of (3) reads: Find uH01, such that

abuvdx+abquvdx=abfvdx,vH01.E4

We want to find uH01 that satisfies (4). We note that a solution u to (4) is less regular than the solution u (3). Indeed, (4) has only u whereas (3) contains u. 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 uC2ab, 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.

5. FE method gives an approximation of the weak solution.

From now on, we use the notation v=vΩ, where Ω=ab.

The FE formulation: The FE method is based on the variational form (4). We note that the space H01 contains many functions and it is therefore just as hard to find a function uH01 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=x0<x1<<xN=b be a regular partition of ab. Suppose that the length of Ii=xi1xi is hi=xixi1. We define h=maxi=1,2,,Nhi to be the mesh size. We wish to construct a subspace VhV=H01. Since Vh has finite dimension, we can find one set of basis functions ϕjj=1N1 for Vh, where ϕjVh,j=1,2,,N1 are linearly independent. We remark that Vh is the space spanned by the basis functions i.e.,Vh=vhxvhx=j=1N1cjϕjx. The FE method consists of choosing a basis for the subspace Vh 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. uh must converge to the solution u of the original problem as h0.

It is natural to obtain an approximation uh to u as follows: Find uhVh such that

abuhvdx+abquhvdx=abfvdx,vVh.E5

We call uh the FE approximation of u. We say that (5) is the Galerkin approximation of (4) and the method used to find uhVh is called Galerkin method.

FE approximation using LagrangeP1elements: The simplest finite dimensional space is the piecewise continuous linear function space defined over the triangulation

Vh,01=vhVvhis piecewise continuous linear overabwithvha=vhb=0.

It is easy to show that Vh,01 has a finite dimension even although there are infinite number of elements in Vh,01. The approximation of the FE method is therefore to look for an approximation uh within a small (finite dimensional) subspace Vh,01=vVh1va=vb=0 of H01, consisting of piecewise linear polynomials, where Vh1=vC0abvIiP1Ii.

Let Vh,01 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 ϕii=1N1. The simplest ones are the so-called hat functions satisfying ϕixj=δ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 Vh,01 must vanish at the end points x0=a and xN=b.

The explicit expressions for the hat function ϕix and its derivative ϕix are given by

ϕix=0,axxi1,xxi1hi,xi1xxi,xi+1xhi+1,xixxi+1,0,xi+1xb,,ϕix=0,a<x<xi1,1hi,xi1<x<xi,hi+1,xi<x<xi+1,0,xi+1<x<b,,

for i=1,2,,N1. The FE approximation of (4) thus reads: Find uVh,01, such that

abuhvdx+abquhvdx=abfvdx,vVh,01.E6

We call uh the FE approximation of u. We say that (6) is the Galerkin approximation of (4) and the method used to find uhVh,01 is called Galerkin method.

It can be shown that (6) is equivalent to the N1 equations

abuhϕidx+abquhϕidx=abfϕidx,i=1,2,,N1.E7

Derivation of the discrete system: Since uhVh,01, we can express it as a linear combination of hat functions i.e.,

uh=j=1N1cjϕjx,E8

where cj are real numbers to be determined. We note that the coefficients cj,j=1,2,,N1 are the N1 nodal values of uh to be determined. Note that the index is only from 1 to N1, because of the zero boundary conditions. We remark that uha=uhb=0 and uhxi=ci. So ci is an approximate solution to the exact solution at x=xi.

We can use either the weak/variational form (V), or the minimization form (M), to derive a linear system of equations for the coefficients cj.

Substituting (8) into (7) yields

j=1N1cjabϕiϕjdx+abqϕiϕjdx=abfϕidx,i=1,2,,N1.E9

The problem (7) is now equivalent to the following: Find the real numbers c1,c2,,cN1 that satisfy the linear system (9).

We note that the linear system (9) is equivalent to the system in matrix–vector form

Ac=b,E10

where c=c1c2cN1tRN1 is the unknown vector, A is an N1×N1 matrix, the so-called stiffness matrix when q=0, with entries

aij=abϕiϕj+qϕiϕjdx,i,j=1,2,,N1,E11

and bRN1, the so-called load vector, has entries

bi=abfϕidx,i=1,2,,N1.E12

To obtain the approximate solution we need to solve the linear system for the unknown vector c. We note that aij=aϕiϕj and bi=fϕi, where auv=abuv+quvdx is a bi-linear and fv=abfvdx is a linear form.

2.3.3 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

minvxH01abFv,whereFv=ab12v2+12qv2fvdx.

As before, we look for an approximate solution of the form (8). If we plug this into the functional form above, we get

Fuh=ab12j=1N1cjϕjx2+12qj=1N1cjϕjx2fj=1N1cjϕjxdx,

which is a multi-variable function of c1,c2,,cN1 and can be written as Fuh=Fc1c2cN1. The necessary conditions for a global minimum are Fci=0, i=1,2,,N1. Taking the partial derivatives directly with respect to ci, we get

abϕixj=1N1cjϕjx+qϕixj=1N1cjϕjxfϕixdx=0,i=1,2,,N1.

Exchange the order of integration and the summation, we get

j=1N1cjabϕixϕjx+qϕixϕjxdx=abfϕixdx=0,i=1,2,,N1,

which is exactly the same linear system (9) obtained using the Galerkin method.

2.3.4 Computer implementation

It is straightforward to calculate the entries âi,j=abϕiϕjdx. For ij>1, we have âi,j=0, since ϕi and ϕj lack overlapping support. However, if i=j, then

âi,i=abϕi2dx=xi1xi1hi2dx+xixi+11hi+12dx=1hi+1hi+1,i,j=1,2,,N1.

Furthermore, if j=i+1, then

âi,i+1=abϕiϕi+1'dx=xixi+11hi+11hi+1dx=1hi+1,i,j=1,2,,N2.E13

By symmetry, we also have

âi+1,i=abϕi+1ϕidx=1hi+1,i,j=1,2,,N2.

To obtain a˜i,j=abqϕiϕjdx and bi=abfϕidx, we use the composite trapezoidal rule

abfxdx=i=1Nxi1xifxdx12h1fx0+i=1N1hi+hi+1fxi+hNfxN.

So, we can easily verify that

a˜i,j=abqϕiϕjdxqi2hi+hi+1,i=j0,ij,bi=abfϕidx12hi+hi+1fi,

where qi=qxi and fi=fxi. Thus, the matrix A=âi,j+a˜i,j is tridiagonal and has the form

A=1h1+1h2+q12h1+h21h2001h21h2+1h3+q22h2+h31h3001h301hN1001hN11hN1+1hN+qN12hN1+hN.

Finally, we obtain the following system: c0=cN=0 and

1hici1+1hi+hi+1ci1hi+1ci+1+qihi+hi+12ci=12hi+hi+1fi,i=1,2,,N1,

Remark 2.2Suppose that the partition is uniformi.e.,hi=h=baNfor alli=1,2,,N. Then the stiffness matrixAand the load vectorbhave the form:

A=2h+hq11h001h2h+hq21h001h01hN1001h2h+hqN1,b=hf1f2f3fN1.

Finally, we obtain the following system: c0=cN=0 and

cj1+2cjcj+1h+hqicj=hficj12cj+cj+1h2+qicj=fi,i=1,2,,N1,

which is the same system obtained using the finite difference method, where u is approximated using the second-order midpoint formula uxjuxj12uxj+uxj+1h2. We conclude that the above FE method using the composite trapezoidal rule is equivalent to the finite difference method of order 2.

2.3.5 Existence, uniqueness, and basic a priori error estimate

Lemma 2.1The matrixAwith entriesai,j=abϕiϕjdxis symmetric positive definitei.e.,ai,j=aj,iand

xtAx=i,j=1N1xiai,jxj>0,for all nonzerox=x1xN1tRN1.

Theorem 2.2The linear system (10) obtained using the FE method has a unique solution. Consequently, the FE method solution uh is unique.

Next, we state a general convergence result for the Galerkin method. We first define the following norm and semi-norm: For vH01, we define

v=abv2xdx1/2,v1=v=abvx2dx1/2.

Theorem 2.3Suppose thatqx0, xab. Letube the solution to (4) anduhbe the solution to (6). Then there exists a constantCsuch that

uuhCuvh,vhVh,01,E14

where C is given by C=1+maxxabqx, which is independent of the choice of Vh,01.

Remark 2.3From (14), taking the minimum over allvhVh,01, we getuuhCminvhVh,0uvh. Thus,uuh1CminvhVh,0uvh1, whereC=1+maxxabqx.

Next, we study the convergence of uh to u. Let uH01. Define the piecewise linear interpolant by

πu=j=1NuxjϕjxVh,01,xab.

Since πuVh,01, the estimate (14) gives

uuhCuπu.

This inequality suggest that the error between u and uh is controlled by the interpolation error uπu in the 1-norm.

Theorem 2.4 (A priori error estimate) Suppose thatqx0xab. Letube the solution to (4) anduhbe the solution to (6). Then there exists a constantCsuch that

uuh2Ci=1Nhi2uIi2,

where C is a constant independent of h. Consequently, if h=maxihi, then

uuh2Ch2u2.

Remark 2.4

1. If the partition is not uniform then we obtain the same error estimate with h=maxi=1,2,,Nxixi1.

2. The error is expressed in terms of the exact solution u. If it is expressed in terms of the computed solution uh it is an a posteriori error estimate (this yields a computable error bound).

3. uhu in the v-norm as h=maxihi0. If uuh=0 then uuh is constant, but since u0=uh0 we also have uuh=0 and therefore uh=u.

4. uh is the best approximation within the space Vh,01 with respect to the v-norm.

5. The norm v is referred to as the energy norm and has often a physical meaning.

2.3.6 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: ua=α and ub=β for two real numbers α and β. This BC is also known as strong BC or essential BC.

2. Neumann BCs: ua=α and ub=β for two real numbers α and β. This BC is also known as natural BCs.

3. Robin BCs: ua=αua and ub=βub for two real numbers α and β.

Note that any combination is possible at the two boundary points.

Nonhomogeneous Dirichlet boundary conditions: Let us consider the following two-point BVP: find uC2ab such that

u=fx,xab,ua=α,ub=β,E15

where α and β are given constants and fCab is a given function. In this case, the admissible function space H01=v:v2+v2<va=vb=0 and the FE space Vh,01 defined earlier remain the same. Multiplying (15) by a test function vH01 and integrating by parts gives

abfvdx=abuvdx=ubvb+uava+abuvdx=abuvdx,

since va=vb=0. Hence, the weak or variational form of (15) reads: Given ua=α, ub=β, find uH1=v:v2+v2<, such that

abuvdx=abfvdx,vH01.E16

Let Vh1 and Vh,01, 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=x0<x1<<xN=b be a uniform partition of the interval ab. Moreover let ϕi be the set of hat basis functions of Vh associated with the N+1 nodes xj,j=0,1,,N, such that ϕixj=δij. The FE approximation of (16) thus reads: Find uhVh1 such that uha=α, uhb=β, and

abuhvdx=abfvdx,vVh,01.E17

It can be shown that (17) is equivalent to the N1 equations

abuhϕidx=abfϕidx,i=1,2,,N1.E18

Expanding uh as a linear combination of hat functions

uh=j=0Ncjϕjx=αϕ0x+j=1N1cjϕjx+βϕNx,E19

where the coefficients cj,j=1,2,,N1 are the N1 nodal values of uh to be determined.

Substituting (19) into (18) yields

j=1N1cjabϕiϕjdx=abfϕiαϕ0ϕiβϕNϕidx,i=1,2,,N1,

which is a N1×N1 system of equations for cj. In matrix form we write

Ac=b,E20

where A is a N1×N1 matrix, the so-called stiffness matrix, with entries

ai,j=abϕiϕjdx,i,j=1,2,,N1,E21

c=c1c2cN1t is a N1 vector containing the unknown coefficients cj,j=1,2,,N1, and b is a N1 vector, the so-called load vector, with entries

bi=abfϕiαϕ0ϕiβϕNϕidx,i=1,2,,N1.E22

Computer Implementation: The explicit expression for a hat function ϕix is given by

ϕix=0,axxi1,xxi1hi,xi1<xxi,xi+1xhi+1,xi<xxi+1,0,xi+1<xb,,i=1,2,,N1,ϕ0x=x1xh1,x0<xx1,0,x1<xb,,ϕNx=0,x0<xxN1,xxN1hN,xN1<xb.

For simplicity we assume the partition is uniform so that hi=h for i=1,2,,N. Hence the derivative ϕix is either 1h, 1h, or 0 depending on the interval.

It is straightforward to calculate the entries of the stiffness matrix. For ij>1, we have ai,j=0, since ϕi and ϕj lack overlapping support. However, if i=j, then

ai,j=abϕi2dx=xi1xi1h2dx+xixi+11h2dx=2h,i,j=1,2,,N1,

where we have used that xixi1=xi+1xi=h. Furthermore, if j=i+1, then

ai,i+1=abϕiϕi+1'dx=xixi+11h1hdx=1h,i,j=1,2,,N2.

Changing i to i1 we also have

ai1,i=abϕi1'ϕi'dx=xi1xi1h1hdx=1h,i,j=2,3,,N1.

Thus the stiffness matrix is

A=1h21001210012010012.

The entries bi 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

b1=abfϕ1αϕ0ϕ1βϕNϕ1dx=x0x1fϕ1α1h1hdx+x1x2fϕ1=αh+x0x2fϕ1αh+hfx1,bi=abfϕiαϕ0ϕiβϕNϕidx=xi1xi+1fϕidxhfxi,i=2,,N2,bN1=abfϕN1αϕ0ϕN1βϕNϕN1dx=xN2xN1fϕN1dx+xN1xNfϕN1β1h1hdx=xN2xNfϕN1dx+βhhfxN1+βh.

Assembly: We rewrite (20), (21), (22) as

1h21001210012010012c1c2c3cN1=hfx1+αhhfx2hfx3hfxN1+βh

We note that uha=α=ua and uhb=β=ub. 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. An alternative approach is to use all the basis functions ϕ0ϕN to 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 uC2ab such that

u=fx,xab,ua=α,γub+ub=β,E23

where α,β and γ are given numbers and fCab is a given function. The boundary condition at x=b is called a Robin boundary condition (combination and u and u is prescribed at x=b). In this case, the admissible function space is modified to

H01=v:v2+v2<va=0.

Multiplying (23) by a function vH01 and integrating by parts gives

abfvdx=abuvdx=ubvb+uava+abuvdx
=βγubvb+uava+abuvdx.

Since va=0, we are left with

abuvdx+γubvb=abfvdx+βvb.

Hence, the weak or variational form of (23) reads: Given ua=α, find the approximate solution uH01, such that

abuvdx+γubvb=abfvdx+βvb,vH01.E24

The FE space Vh1 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 uh to the solution u satisfies

abuhvdx+γuhbvb=abfvdx+βvb,vVh1,E25

with uha=α. As before, (25) can be formulated in matrix form.

2.4 Model problem with coefficient and general Robin BCs

Let us consider the following two-point BVP: find uC2ab such that

pxu=fx,xI=ab,paua=κ0uaα,pbub=κ1ubβ,E26

where p=px with pxp0>0, fL2I, κ0,κ10, and α,β are given numbers. Let

V=vC0I:v2+v2<.

Multiplying (26) by a function vV and integrating by parts gives

abfvdx=abpuvdx=abpuvdxpbubvb+pauava=abpuvdxκ1ubβvb+κ0uaαva.

We gather all u-independent terms on the left and obtain

abpuvdxκ1ubvb+κ0uava=abfvdxκ1βvb+κ0αva,vV.

The FE method consists of finding uhVh=vC0abvIiP1Ii such that

abpuhvdxκ1uhbvb+κ0uhava=abfvdxκ1βvb+κ0αva,vVh.E27

Implementation: We need to assemble a stiffness matrix A and a load vector b. Substituting uh=i=0Nciϕi into (27) and taking v=ϕj for j=0,1,,N yields

i=0Nabpϕiϕjdxκ1ϕibϕjb+κ0ϕiaϕja=abfϕjdxκ1βϕjb+κ0αϕja,j=0,1,,N.

which is a N+1×N+1 system of equations for ci. In matrix form we write Ac=b, where c=c0cNt is a N+1 vector containing the unknown coefficients ci,i=0,1,,N, A is a N+1×N+1 matrix with entries

ai,j=abpϕiϕjdxκ1ϕibϕjb+κ0ϕiaϕja,i,j=0,1,,N,

and b is a N+1 vector with entries

bj=abfϕjdxκ1βϕjb+κ0αϕja,j=0,1,,N.

Let for simplification p=1. Then the matrix A and the vector b (when using the trapezoidal rule) are given by

A=κ0+1h11h1001h11h1+1h21h2001h201hN001hN1hNκ1,b=h12f0+κ0αh1+h22f1hN1+hN2fN1hN2fNκ1β.

2.5 The FE method using Lagrange P2 elements

Let a=x0<x1<<xN=b be a regular partition of the interval ab. Suppose that the length of Ii=xi1xi is hi=xixi1. Let Pk=px=j=0kcjxjcjR denotes the vector space of polynomials in one variable and of degree less than or equal to k. The FE method for Lagrange P2 elements involves the discrete space:

Vh2={vxC0[ab]vIiP2Ii,i=1,,N},

and its subspace V0,h2=vVh2va=vb=0. These spaces are composed of continuous, piecewise parabolic functions (polynomials of degree less than or equal to 2). The P2 FE method consists in applying the internal variational approximation approach to these spaces.

Lemma 2.2The spaceVh2is a subspace ofH1abof dimension2N+1. Every functionvhVh2is uniquely defined by its values at the mesh verticesxj,j=0,1,,Nand at the midpointsxj+12=xj+xj+12=xj+hj+12,j=0,1,,N1, wherehj+1=xj+1xj:

vhx=j=0Nvhxjϕjx+j=0N1vhxj+12ϕj+12x,xab,

where ϕjj=0N is the basis of the shape functions ϕj defined as:

ϕjx=ϕxxjhj+1,j=0,1,,N,ϕj+12x=ψxxj+12hj+1,j=0,1,,N1,

with

ϕξ=1+ξ1+2ξ,ξ10,1ξ12ξ,ξ01,0,ξ>1,ψξ=14ξ2,ξ12,0,ξ>12,E28

Figure 3 shows the global shape functions for the space Vh2 and the three quadratic Lagrange P2 shape functions on the reference interval 11.

Remark 2.5Notice that we have:

ϕjxj=δij,ϕjxj+12=0,ϕj+12xj=0,ϕj+12xj+12=δij.

Corollary 2.1The spaceV0,h2is a subspace ofH01abof dimension2N1and every functionvhV0,h2is uniquely defined by its values at the mesh verticesxj,j=1,2,,N1and at the midpointsxj+12,j=0,1,,N1:

vhx=j=1N1vhxjϕjx+j=0N1vhxj+12ϕj+12x,xab,

where ϕjj=0N is the basis of the shape functions ϕj defined as:

ϕjx=ϕxxjhj+1,j=0,1,,N,ϕj+12x=ψxxj+12hj+1,j=0,1,,N1,

with ϕξ and ψξ are defined by (28).

2.5.1 Homogeneous boundary conditions

The variational formulation of the internal approximation of the Dirichlet BVP (3) consists now in finding uhV0,h2, such that:

abuhvdx+abquhvdx=abfvdx,vVh,02.

Here, it is convenient to introduce the notation xj2,j=1,,2N1 for the mesh points and ϕj2,j=1,,2N1 for the basis of V0,h2. Using these notations, we have:

uh=j=12N1cj2ϕj2x,

where cj2=uhxj2uxj2 are the unknowns coefficients. This formulation leads to solve in R2N1 a linear system:

Ac=b,

where c=c12c1cN12tR2N1 is the unknown vector containing the coefficients cj2,j=1,2,,2N1, A is an 2N1×2N1 matrix with entries

aij=abϕi2'ϕj2'+qϕi2ϕj2dx,i,j=1,2,,2N1,

and load vector bR2N1 has entries

bi2=abfϕi2dx,i=1,2,,2N1.

Since the shape functions ϕi have a small support, the matrix A is mostly composed of zeros. However, the main difference with the Lagrange P1 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 ξ11:

x=xj+xj12+xjxj12ξ=xj12+xjxj12ξ,xxj1xj,j=1,2,,N.

Hence, the shape functions can be reduced to only three basic shape functions (Figure 3):

ϕ̂1ξ=ξξ12,ϕ̂0ξ=1ξ1+ξ,ϕ̂1ξ=ξξ+12.

Their respective derivatives are

dϕ̂1ξ=2ξ12,dϕ̂0ξ=2ξ,dϕ̂1ξ=2ξ+12.

This approach consists in considering all computations on an interval Ii=xi1xi on the reference interval 11. Thus, we have:

dϕixdx=dϕixi1/2+xixi12ξdx=2xixi1dϕ̂kξ=2hidϕ̂kξ.

In this case, the elementary contributions of the element Ii to the stiffness matrix and to the mass matrix are given by the 3×3 matrices KIi and MIi:

KIi=Iiϕi1ϕi1ϕi1ϕi12ϕi1ϕiϕi12ϕi1ϕi12ϕi12ϕi12ϕiϕiϕi1ϕiϕi12ϕiϕidx=2hi11ϕ̂1ϕ̂1ϕ̂1ϕ̂0ϕ̂1ϕ̂1ϕ̂0ϕ̂1ϕ̂0ϕ̂0ϕ̂0ϕ̂1ϕ̂1ϕ̂1ϕ̂1ϕ̂0ϕ̂1ϕ̂1=13hi7818168187,
MIi=Iiϕi1ϕi1ϕi1ϕi12ϕi1ϕiϕi12ϕi1ϕi12ϕi12ϕi12ϕiϕiϕi1ϕiϕi12ϕiϕidx=hi211ϕ̂1ϕ̂1ϕ̂1ϕ̂0ϕ̂1ϕ̂1ϕ̂0ϕ̂1ϕ̂0ϕ̂0ϕ̂0ϕ̂1ϕ̂1ϕ̂1ϕ̂1ϕ̂0ϕ̂1ϕ̂1=hi304212162124.

Coefficients of the right-hand sideb: Usually, the function f is only known by its values at the mesh points xi2,i=0,1,,2N and thus, we use the decomposition of f in the basis of shape functions ϕi2,i=0,1,,2N as fx=j=02Nfxj2ϕj2. Each component bi2 of the right-hand side vector is obtained as bi2=k=1Nxk1xkfϕi2dx. Using the previous decomposition of f, we obtain:

bi2=k=1Nxk1xkj=02Nfxj2ϕj2ϕi2dx=j=02Nfxj2k=1Nxk1xkϕi2ϕj2dx.

Thus, the problem is reduced to computing the integrals xk1xkϕi2ϕj2dx. It is easy to see that we obtain expressions very similar to that of the mass matrix. More precisely, the element Ii=xi1xi will contribute to only three components of indices i1, i12 and i as:

bIi=hi304212162124fxi1fxi12fxi..

2.5.2 Nonhomogeneous boundary conditions

Consider the following two-point BVP: find uC2ab such that

u+qxu=fx,xab,ua=α,ub=β,E29

where α and β are given constants and fCab is a given function.

Multiplying (29) by a function vH01=v:v2+v2<va=vb=0 and integrating by parts gives

abfvdx=abu+quvdx=ubvb+uava+abuv+quvdx=abuvdx.

Hence, the weak or variational form of (29) reads: Given ua=α, ub=β, find uH1=v:v2+v2<, such that

abuv+quvdx=abfvdx,vH01.

Let Vh2 and Vh,02, 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=x0<x1<<xN=b of the interval ab.

The FE method scheme consists of finding uhVh2, such that:

abuhvdx+abquhvdx=abfvdx,vVh,02.

Introduce the notation xj2,j=0,1,,2N1,2N for the mesh points and ϕj2,j=0,1,,2N1,2N for the basis of Vh2 and ϕj2,j=1,,2N1 for the basis of V0,h2. Using these notations, we have:

uh=j=02Ncj2ϕj2x,

where cj2=uhxj2uxj2 are the unknowns coefficients. We note that c0=uhx0=α and c2N=uhxN=β. This formulation leads to solve in R2N1 a linear system:

Ac=b,

where c=c12c1cN12tR2N1 is the unknown vector containing the coefficients cj2,j=1,2,,2N1, A is an 2N1×2N1 matrix with entries

aij=abϕi2ϕj2+qϕi2ϕj2dx,i,j=1,2,,2N1,

and the load vector bR2N1 has entries

bi2=abfϕi2dxαabϕi2ϕ0+qϕi2ϕ0dxβabϕi2ϕN+qϕi2ϕNdx,i=1,2,,2N1.

Clearly, the only extra terms are given in the vector with entries

b˜i2=αabϕi2ϕ0+qϕi2ϕ0dxβabϕi2ϕN+qϕi2ϕNdx,i=1,2,,2N1.

Suppose q=0 then for N2, we have

b˜12=αabϕ12ϕ0dxβabϕ12ϕNdx=αx0x1ϕ12ϕ0=8α3h1,b˜1=αabϕ1ϕ0dxβabϕ1ϕNdx=αx0x1ϕ1ϕ0=α3h1,b˜i2=αabϕi2ϕ0dxβabϕi2ϕNdx=0,i=3,,2N3,b˜N1=αabϕN1ϕ0dxβabϕN1ϕNdx=βxN1xNϕN1ϕNdx=β3h1,b˜N12=αabϕN12ϕ0dxβabϕN12ϕNdx=βxN1xNϕN12ϕNdx=8β3h1.

3. The FE for elliptic PDEs

Here, we apply the FE method for two-dimensional elliptic problem: Find u such that

au+bu=fx,xΩ,aun=κgu,on∂Ω,E30

where a>0, b0, κ0, fL2Ω and gC0∂Ω.

3.1 Meshes

Let ΩR2 bounded with ∂Ω assumed to be polygonal. A triangulation Th of Ω is a set of triangles T such that Ω=TThT, and two triangles intersect by either a common triangle edge, or a corner, or nothing. Corners will be referred to as nodes. We let hT=diamT the length or the largest edge.

Let Th have N nodes and M triangles. The data is stored in two matrices. The matrix PR2×N describes the nodes (x1y1,,xNyN) and the matrix KR3×M describes the triangles, i.e., it describes which nodes (numerated from 1 to N) form a triangle T and how it is orientated:

P=x1x2xNy1y2yN,K=n1αn2αnMαn1βn2βnMβn1γn2γnMγ.

This means that triangle Ti is formed by the nodes niα, niβ, and niγ (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.

3.2 Piecewise polynomial spaces

Let T be a triangle with nodes N1=x1y1, N2=x2y2, and N3=x3y3. We define

P1T=vC0Tvxy=c1+c2x+c3yc1c2c3R.

Now let vi=vNi for i=1,2,3. Note that vP1T is determined by vii=13. Given vi we compute ci by

1x1y11x2y21x3y3c1c2c3=v1v2v3.

This is solvable due to

det1x1y11x2y21x3y3=2T0,1x1y11x2y21x3y31=12Tx2y3x3y2x3y1x1y3x1y2x2y1y2y3y3y1y1y2x3x2x1x3x3x1,

where T=12x2y3x3y2x1y3+x3y1+x1y2x2y1, which is ± the area of the triangle T.

Let λjP1T be given by the nodal values λjNi=δij, where δij is the Kronecker symbol. This gives us vxy=α1λ1xy+α2λ2xy+α3λ3xy, where αi=vNi for i=1,2,3. We can compute λixy as follows: Let λixy=ai+bix+ciy. Using λjNi=δij, we get

1x1y11x2y21x3y3a1b1c1=100,1x1y11x2y21x3y3a2b2c2=010,1x1y11x2y21x3y3a3b3c3=001.

Solving the systems, we get

λ1xy=12Tx2y3x3y2+y2y3x+x3x2y,
λ2xy=12Tx3y1x1y3+y3y1x+x1x3y,
λ3xy=12Tx1y2x2y1+y1y2x+x3x1y.

Let Th be a triangulation of Ω, then we let

Vh=vCΩvTP1TTTh.

Functions in Vh are piecewise linear and continuous. We know that vVh is uniquely determined by vNii=12N. We let ϕjNi=δij and let ϕjj=12NVh be a basis for Vh (hat functions), i.e.,

vxy=i=1Nαiϕixy,αi=vNi,i=1,2,,N.

3.3 Interpolation

Given uCT on a single triangle with nodes Ni=xiyi, i=1,2,3, we let

πuxy=i=13uNiϕixy,

in particular πuNi=uNi, i=1,2,,N. We want to estimate the interpolation error uπu. Let

uL2Ω2=Ωux2dxdy,DuL2Ω2=uxL2Ω2+uyL2Ω2,D2uL2Ω2=uxxL2Ω2+2uxyL2Ω2+uyyL2Ω2.

Theorem 3.1Suppose thatuC2T. Then the following hold

uπuL2TChT2D2uL2T,DuπuL2TChTD2uL2T,

where C is a generic constant independent of hT 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=i=1NuNiϕi.

Theorem 3.2Suppose thatuC2Tfor allTTh. Then the following hold

uπuL2Ω2CTThhT4D2uL2T2,DuπuL2Ω2CTThhT2D2uL2T2,

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 Th. Here DuπuL2Ω2=TThDuπuL2T2.

3.4 L2-projection

Let ΩR2. We consider the space L2Ω=vΩv2xydxdy<. Let uL2Ω. We define the L2-projection Ph:L2ΩVh=vC0ΩvTP1TTTh by PhuVh such that

ΩuPhuvhdxdy=0,vhVh.

The problem of finding PhuVh is equivalent to solve the following linear system

ΩuPhuϕidxdy=0,i=1,2,N,

where ϕii=1N is a basis of Vh.

Since PhuVh we can express it as Phu=i=1Nciϕixy, where ciR. Therefore, to find PhuVh we need to find c1,c2,,cNR such that

i=1NciΩϕiϕjdxdy=Ωuϕjdxdy,j=1,2,,N.

The problem can be expressed as a linear system of equations Mc=b, where c=c1c2cNt and the entries of the matrix MRN×N and the vector bRN are given by

mij=Ωϕiϕjdxdy,bj=Ωuϕjdxdy.

In general, we use a quadrature rule to approximate integrals. The general form is

Tfxydxdyj=1nωjfN¯j,

where the ωj ′s denote the weights and the N¯j ′s the quadrature points.

Lemma 3.1The mass matrixMwith entriesmij=Ωϕiϕjdxdyis symmetric and positive definite.

Theorem 3.3For anyuL2ΩtheL2-projectionPhuexists and is unique.

3.5 A priori error estimate

Theorem 3.4LetuL2Ωand letPhube theL2-projection ofu, then

uPhuL2ΩuvhL2Ω,vhVh.

Theorem 3.5Suppose thatuC2ΩwithuC2Tfor allTTh. Then there exists a constantCsuch that

uPhuL2Ω2CTThhT4D2uL2T2.

3.6 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

au+bu=f,inΩ,aun=κgu,on∂Ω,E31

where a>0, b0, κ0, fL2Ω and gC0∂Ω. We seek a weak solution u in

V=H1Ω=vL2Ωvhas a weak derivative andvL2Ω+vL2Ω<.

In order to derive the weak formulation, we multiply (31) with vV, integrate over Ω and use Green’s formula to obtain

Ωfvdxdy=Ωvaudxdy+Ωbuvdxdy=Ωauvdxdy∂Ωvaunds+Ωbuvdxdy=Ωauvdxdy+Ωbuvdxdy+∂Ωκugvds.

We obtain the weak form: Find uV such that

Ωauvdxdy+Ωbuvdxdy+∂Ωκuvds=Ωfvdxdy+∂Ωκgvds,vV.E32

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 uhVh such that

Ωauhvhdxdy+Ωbuhvhdxdy+∂Ωκuhvhds=Ωfvhdxdy+∂Ωκgvhds,vhVh,E33

where Vh=vVvTP1TTTh.

Implementation: Let a=1 and b=g=0. Substituting uh=j=1Ncjϕj into (33) and picking vh=ϕi, we obtain

j=1NcjΩϕjϕidxdy+∂Ωκϕjϕids=Ωfϕidxdy,i=1,2,,N.

This gives us the system A+Rc=b, where c=c1c2cNtRN is the unknown vector and the entries of ARN×N, RRN×N, and bRN are given by

aij=Ωϕjϕidxdy,rij=∂Ωκϕjϕids,bi=Ωfϕidxdy,i,j=1,2,,N.

Assembly of the stiffness matrixA: We can again identify the local contributions that come form a particular triangle T

aijT=Ωϕjϕidxdy,i,j=1,2,3.

where T is an arbitrary triangle with vertices Ni=xiyi and ϕi are the hat functions i.e.,ϕjNi=δij. Let ϕixy=αi+βix+γiy, for i=1,2,3. Then, we compute αi,βi,γi by

1x1y11x2y21x3y3α1β1γ1=100,1x1y11x2y21x3y3α2β2γ2=010,1x1y11x2y21x3y3α3β3γ3=001.

In general we have Bαi=ei for i=1,2,3, where

B=1x1y11x2y21x3y3,αi=αiβiγi,e1=100,e2=010,e3=001.

Furthermore, we obviously have ϕi=βiγit, which gives

aijT=Ωβiβj+γiγjdx=βiβj+γiγjT,i,j=1,2,3.

Assembly of boundary matrixR: Let Γhout denote the set of boundary edges of the triangulation, i.e.Γhout=EE=T∂ΩforTTh. Assume that κ is constant on E. For an edge EΓhout, we define RER2×2 by the entries

rijE=Eκϕjϕids=κ61+δijE,i,j=1,2,

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 TTh

biT=TfϕidxdyT3fNi,i=1,2,,N.

Given A, R and b, we can solve A+Rc=b and write uh=j=1Ncjϕj.

3.7 The Dirichlet problem

Consider the following Dirichlet Problem: Find u such that

Δu=f,inΩ,u=g,on∂Ω,E34

where fL2Ω and gC0∂Ω. We seek a weak solution u in Vg=vVv∂Ω=g. Multiplying (34) by a test function vV0 and integrating over Ω, we get

Ωfvdxdy=ΩvΔudxdy=Ωuvdxdy∂Ωvunds=Ωuvdxdy.

So the weak problem reads: Find uVg such that

Ωuvdxdy=Ωfvdxdy,vV0.

Assume that g is piecewise linear on ∂Ω with respect to the triangulation. Then the FE method in 2D is defined as follows: Find uhVh,g=vVhv∂Ω=g such that

Ωuhvhdxdy=Ωfvhdxdy,vhVh,0.

Assume that we have N nodes and J boundary nodes, then the matrix form of the FE method problem reads:

A0,0A0,gAg,0A0,gc0c1=b0b1,

where A0,0RNJ×NJ, Ag,gRJ×J, A0,gRNJ×J, Ag,0RJ×NJ. Note that c1RJ is known (it contains the values of g in the boundary nodes). We can therefore solve the simplified problem reading: find c0RNJ with A0,0c0=b0A0,gc1.

3.8 The Neumann problem

Consider the following Neumann Problem: Find u such that

Δu=f,inΩ,un=g,on∂Ω,E35

where fL2Ω and gC0∂Ω. Let us try to seek a solution to this problem in the space V=vvL2Ω+vL2Ω<. Multiplying (35) by a test function vV, integrating over Ω, and using Green’s formula, we get

Ωfvdxdy=ΩvΔudxdy=Ωuvdxdy∂Ωvunds=Ωuvdxdy∂Ωvgds.

Thus, the variational formulation reads: find uV such that

Ωuvdxdy∂Ωvgds=Ωfvdxdy,vV.

In order to guarantee solvability, we note that if v=1 then we have

0=Ωu1dxdy=Ωfdxdy+∂Ωgds.

Therefore we need to assume the following compatibility condition

Ωfdxdy+∂Ωgds=0,

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 cR. 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 space

V̂=vVΩvdxdy=0,

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 uV̂ such that

Ωuvdxdy∂Ωvgds=Ωfvdxdy,vV.

Now, the FE method takes the form: find uhV̂hV̂ such that

Ωuhvhdxdy∂Ωvhgds=Ωfvhdxdy,vhV̂h,

where V̂h is the space of all continuous piecewise linear functions with a zero mean.

3.9 Finite elements for mixed Dirichlet-Neumann conditions

Here we describe briefly how Neumann conditions are handled in two-dimensional finite elements. Suppose Ω is a domain in either R2 or R3 and assume that ∂Ω has been partitioned into two disjoint sets: ∂Ω=Γ1Γ2. We consider the following BVP:

κxu=fx,xΩ,u=0,xΓ1,un=0,xΓ2,E36

where fL2Ω. 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 by

V̂=vC2Ω¯:vx=0xΓ1.

Multiplying (36) by a test function vV̂ and integrating over Ω, we get

Ωfvdxdy=Ωvκxudxdy=Ωκxuvdxdy∂Ωκxvunds=ΩκxuvdxdyΓ1κxvundsΓ2κxvunds=Ωκxuvdxdy,

since v=0 on Γ1 and un on Γ1. Thus the weak form of (36) is: Find uV̂ such that

Ωκxuvdxdy=Ωfvdxdy,vV̂.E37

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 of V̂ as follows:

Vh=vCΩ¯:v is linear on Thvz=0for all nodeszΓ1.

A basis for Vh 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.

3.10 The method of shifting the data

3.10.1 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.

The method of shifting the data can be used to transform an inhomogeneous Dirichlet problem to a homogeneous Dirichlet problem. This technique works just as it did for a one-dimensional problem, although in two dimensions it is more difficult to find a function satisfying the boundary conditions. We consider the BVP

Δu=fx,xΩ=0a×0b,ux=gx=g1x,xΓ1,g2x,xΓ2,g3x,xΓ3,g4x,xΓ4,E38

where Γ1, Γ2, Γ3, and Γ4 are, respectively, the bottom, right, top, and left boundary edges of the rectangular domain Ω=0a×0b. We will assume that the boundary data are continuous, so

g10=g40,g1a=g20,g2b=g3a,g30=g4b.

Suppose we find a function w defined on Ω¯ and satisfying wx=gx for all x∂Ω. We then define v=uw and note that

Δv=Δu+Δw=fx+Δw=f̂x,

and vx=uxwx=0 for all x∂Ω. We can then solve

Δv=f̂x,xΩ,vx=0,x∂Ω.E39

Finally, the solution u will be given by u=v+w.

We now describe a method for computing a function w that satisfies the given Dirichlet conditions. We first note that there is a polynomial of the form qxy=c0+c1x+c2y+c3xy, which assumes the desired boundary values at the corners:

q00=g10=g40,qa0=g1a=g20,qab=g2b=g3a,q0b=g30=g4b.

A direct calculation shows that

c0=g10,c1=g1ag10a,c2=g4bg40b,c3=g2b+g10g1ag4bab.

We then define

hx=h1x=g1xg10+g1ag10ax,xΓ1,h2y=g2yg20+g2bg20by,xΓ2,h3x=g3xg30+g3ag30ax,xΓ3,h4y=g4yg40+g4bg40by,xΓ4.

We have thus replaced each gi by a function hi which differs from gi by a linear function, and which has value zero at the two endpoints:

h10=h1a=h20=h2b=h30=h3a=h40=h4b=0.

Finally, we define

wxy=c0+c1x+c2y+c3xy+h1x+h3xh1xby+h4y+h2yh4yax.

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 gi ′s into the hi ′s. The first term in the formula for w undoes this transformation. It is straightforward to verify that w satisfies the desired boundary conditions.

3.10.2 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

Δu=fx,xΩ=0a×0b,nux=gx=g1x,xΓ1,g2x,xΓ2,g3x,xΓ3,g4x,xΓ4,E40

where Γ1, Γ2, Γ3, and Γ4 are, respectively, the bottom, right, top, and left boundary edges of the rectangular domain Ω=0a×0b. We first note that this is equivalent to

uyx=g1x,xΓ1,uxx=g2y,xΓ2,uyx=g3x,xΓ3,uxx=g4y,xΓ4.

We make the following observation: If there is a twice-continuously differentiable function u satisfying the given Neumann conditions, then, since uxy=uyx, we have

uxyx0=g1x,uyx0y=g4y,

which together imply that g10=g40. By similar reasoning, we have all of the following conditions:

g10=g40,g10=g40,g1a=g20,g2b=g3a.E41

We will assume that (41) holds.

We now explain how to compute a function that satisfies the desired Neumann conditions. The method is similar to that used to shift the data in a Dirichlet problem: we will “interpolate” between the Neumann conditions in each dimension and arrange things so that the two interpolations do not interfere with each other. We use the fact that

ψx=αx+α+β2ax2satisfiesψ0=α,ψa=β.E42

The first step is to transform the boundary data glx to a function h1x satisfying h10=h1a=0, and similarly for g2, g3, g4 and h2, h3, h4. Since these derivatives of the boundary data at the corners are (plus or minus) the mixed partial derivatives of the desired function at the corners, it suffices to find a function qxy satisfying the conditions

uxy00=g10,uxya0=g1a,uxy0b=g30,uxyab=g2b.

We can satisfy these conditions with a function of the form qxy=c0xy+c1x2y+c2xy2+c3x2y2. The reader can verify that the necessary coefficients are

c0=g10,c1=g10g1a2a,c2=g30+g102b,c3=g2b+g1ag30g104ab.

If w is to satisfy the desired Neumann conditions, then wq=hi on Γi, i=14, where

h1x=g1x+c0x+c1x2,h2y=g2yc0+2ac1yc2+2ac3y2,
h3x=g3xc0+2bc2xc1+2bc3x2,h4y=g4y+c0y+c2y2.

We can now define wq by the interpolation described by (42):

wxy=qxyh1xy+h3x+h1x2by2h4yx+h2y+h4y2ayx2.

Then w satisfies the original Neumann conditions, as the interested reader can verify directly.

3.11 Eigenvalue problem

Consider the following Eigenvalue Problem: Find λR and u such that

Δu=λu,inΩ,un=0,on∂Ω.E43

In order to derive the weak formulation, we multiply (43) with vV, integrate over Ω and use Green’s formula to obtain

λΩuvdxdy=ΩvΔudxdy=Ωuvdxdy∂Ωvunds=Ωuvdxdy.

We obtain the weak form: Find uV such that

Ωuvdxdy=λΩuvdxdy,vV.E44

The FE method in 2D is defined as follows: Find λhR and uhVh such that

Ωuhvhdxdy=λhΩuhvhdxdy,vhVh,E45

where Vh=vVvTP1TTTh.

Implementation: Substituting uh=j=1Ncjϕj into (45) and picking vh=ϕi, we obtain

j=1NcjΩϕjϕidxdyλhΩϕiϕjdxdy=0,i=1,2,,N.

This leads to an algebraic system of the form Ac=λhMc, i.e. an algebraic eigenvalue problem.

3.12 Error analysis

Consider the following model Problem: Find u such that

Δu=f,inΩ,u=0,on∂Ω.

The weak form: Find uV0 such that

Ωuvdxdy=Ωfvdxdy,vV0.

The FE approximation is defined as follows: Find uhVh,0 such that

Ωuhvhdxdy=Ωfvhdxdy,vhVh,0,

where Vh=vVvTP1TTTh. Expressing uh=j=1Ncjϕj and picking vh=ϕi, we obtain

j=1NcjΩϕjϕidxdy=Ωfϕidxdy,i=1,2,,N.

This leads to system of the form Ac=b, where the entries of ARN×N and bRN are

aij=Ωϕjϕidxdy,bi=Ωfϕidxdy,i,j=1,2,,N.

Theorem 3.6The stiffness matrixAis symmetric and positive definite.

Theorem 3.7 (Galerkin orthogonality)LetuV0denote the weak solution anduhVh,0the corresponding FE method approximation. Then

Ωuuhvhdxdy=0,vhVh,0.

Now, let v2=Ωvvdxdy=Ωv2dxdy be the energy norm on V0.

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)LetuV0denote the weak solution anduhVh,0the corresponding FE method approximation. Then

uuhuvh,vhVh,0.

Theorem 3.9LetuV0denote the weak solution anduhVh,0the corresponding FE method approximation. IfuC2Ω, then there existsCindependent ofhTandusuch that

uuhL2Ω2CTThhT2D2uL2T2.

3.13 The FE method for elliptic problems with a convection term

Consider the following convection-diffusion problem: Find u such that

au+bu+cu=f,inΩ,u=0,on∂Ω.E46

We seek a weak solution u in V0=vVv∂Ω=0. In order to derive the weak formulation, we multiply (46) with vV0, integrate over Ω and use Green’s formula to obtain

Ωfvdxdy=Ωvaudxdy+Ωvbudxdy+Ωcuvdxdy=Ωauvdxdy∂Ωvunds+Ωvbudxdy+Ωcuvdxdy=Ωauvdxdy+Ωvbudxdy+Ωcuvdxdy.

Note that there is no need to apply Green’s formula to Ωvbudxdy. We obtain the weak form: Find uV0 such that

Ωauvdxdy+Ωvbudxdy+Ωcuvdxdy=Ωfvdxdy,vV0.

The FE method in 2D is defined as follows: Find uhVh,0=vVhv∂Ω=0 such that

Ωauhvhdxdy+Ωvhbuhdxdy+Ωcuhvhdxdy=Ωfvhdxdy,vhVh,0,E47

where Vh=vVvTP1TTTh.

Implementation: Substituting uh=j=1Ncjϕj into (47) and picking vh=ϕi, we obtain

j=1NcjΩaϕjϕidxdy+Ωϕibϕjdxdy+Ωcϕiϕjdxdy=Ωfϕidxdy,i=1,2,,N.

This gives us the system A+B+Cc=d, where c=c1cNtRN is the unknown vector and the entries of A,B,CRN×N and dRN are given by

aij=Ωaϕjϕidxdy,bij=Ωϕibϕjdxdy,cij=Ωcϕiϕjdxdy,di=Ωfϕidxdy,

for i,j=1,2,,N. Note that B is not symmetric, i.e.bijbji.

4. The FE method for the heat equation

Consider the following heat/diffusion problem: Find uxt such that

u̇Δu=f,inΩR2,t0T,E48
ut=0,on∂Ωandt0T,E49
ux0=u0x,forxΩandt=0.E50

We seek a weak solution u in V0=vv+v<v∂Ω=0. In order to derive the weak formulation, we multiply (48) with vV0, integrate over Ω and use Green’s formula to obtain, for t0T,

Ωfvdx=Ωu̇vdx+Ωuvdx∂Ωvunds=Ωu̇vdx+Ωuvdx.

The weak form therefore reads: Find utV0 such that for t>0

Ωu̇vdx+Ωuvdx=Ωfvdx,vV0.E51

The semi-discrete FE method in 2D is defined as follows: Find uhtVh,0=vVhv∂Ω=0 such that

Ωu̇hvhdx+Ωuhvhdx=Ωfvhdx,vhVh,0,E52

where Vh=vVvTP1TTTh.

Implementation: Substituting uhxt=j=1Ncjtϕjx into (52) and choosing vh=ϕi, we obtain

j=1NċjΩϕjϕidx+j=1NcjΩϕjϕidx=Ωfϕidx,i=1,2,,N.

This gives us the system of ODEs

Mċt+Atct=bt,t0T,c0=c0,

where c=c1c2cNt=uhN1tuhNNttRN (here Ni denotes the node that belongs to the basis function ϕi) is the unknown vector and the entries of M,ARN×N and bRN are given by

mij=Ωϕiϕjdx,aij=Ωϕjϕidx,bi=Ωfϕidx,i,j=1,2,,N.

Finally, the system of ODEs can be solved with e.g., the backward Euler method as follows: Let 0=t0<t1<<tM=T be a discretization, let km=tmtm1 for m=1,2,,M be the time step size and let cmctm for m=1,2,,M denote corresponding approximations. Then, we can compute cm using

M+kmAmcm=Mcm1+kmbm,m=1,2,,M,

where c0 is obtained from u0x. We can either use c0=c10cN0t=u0N1u0NNt, or we can let c0 to be the L2-projection of u0. We set uh0=j=1Ncj0ϕjx and solve for cj0 using

j=1Ncj0Ωϕjϕidx=Ωu0ϕidx,i=1,2,,N.

Theorem 4.1 (Stability)There hold continuous and discrete stability estimates

utu0+0tfsds,uhmuhm1+kmfmuh0+i=1mkifi.

5. The FE method for the wave equation

Many physical phenomena exhibit wave characteristics. For instance light which is an electromagnetic wave have the ability to disperse and create diffraction patterns, which is typical of waves.

Consider the following wave problem: Find uxt such that

u¨εu=f,inΩR2,t0T,E53
nut=0,on∂Ωandt0T,E54
ux0=u0x,u̇x0=v0x,forxΩandt=0,E55

where f is a given load, ε=εxt is a positive parameter, u0 and v0 are a prescribed initial conditions, and Ω is a bounded domain with boundary ∂Ω and unit outward normal n.

We seek a weak solution u in V=H1Ω=vv+v<. Multiplying the wave Eq. (53) with vV, integrating over Ω, and using Green’s formula, we obtain, for t0T,

Ωfvdx=Ωu¨vdxΩvεudx=Ωu¨vdx+Ωεuvdx∂Ωunds=Ωu¨vdx+Ωεuvdx.

The weak form (variational formulation) therefore reads: Find utV=H1Ω such that for all t>0

Ωu¨vdx+Ωεuvdx=Ωfvdx,vV.E56

Let Vh=vVvTP1TTThV be the space of all continuous piecewise linear functions on a triangle mesh of Ω. The semi-discrete FE method in 2D is defined as follows: Find uhtVh such that

Ωu¨hvhdx+Ωεuhvhdx=Ωfvhdx,vhVh.E57

Implementation: Substituting uhxt=j=1Ncjtϕjx into (57) and choosing vh=ϕi, we obtain

j=1Nc¨jΩϕjϕidx+j=1NcjΩεϕjϕidx=Ωfϕidx,i=1,2,,N.

This gives us the system

Mc¨t+Atct=bt,t0T,E58

where c=c1cNt=uhN1tuhNNttRN (here Ni denotes the node that belongs to the basis function ϕi) is the unknown vector and the entries of the mass and stiffness matrices M,ARN×N and the load vector bRN are given by

mij=Ωϕiϕjdx,aij=Ωεϕjϕidx,bi=Ωfϕidx,i,j=1,2,,N.

Eq. (58) is a semi-discretization of the wave equation in the sense that it does not contain any unknowns with spatial derivatives.

Time discretization: We first transform the system of ODEs into a first-order system. Let dt=ċt, we get the new coupled system

MċtMdt=0,Mḋt+Atct=bt,t0T.

Let w=cdt then the system is equivalent to M̂ẇt+Âtwt=b̂t, t0T, where

M̂=M00M,Â=0MA0,b̂=0b.

Finally, the system of ODEs can be solved with e.g., the backward Euler method as follows: Let 0=t0<t1<<tM=T be a discretization, let km=tmtm1 for m=1,2,,M be the time step size and let wmwtm for m=1,2,,M denote corresponding approximations. Then, we can compute wm using

M̂+kmÂmwm=M̂wm1+kmb̂m,m=1,2,,M,

where w0 is obtained from u0x and v0x.

There are several possible choices of initial data. We can either use w0=w10c2N0t=u0N1u0NNv0N1v0NNt, or we can let w0=w10w20t, where w10 and w20 are the L2-projection of u0 and v0, respectively. We set wh,10=j=1Nwj,10ϕjx and wh,20=j=1Nwj,20ϕjx and solve for wj,10,wj,20 using

j=1Nwj,10Ωϕjϕidx=Ωu0ϕidx,j=1Nwj,20Ωϕjϕidx=Ωv0ϕidx,i=1,2,,N.

We can also use Crank–Nicolson scheme

M̂+km2Âmwm=M̂km2Âm1wm1+km2b̂m1+b̂mgm.

Theorem 5.1 (Conservation of energy)Iff=0, then

u̇htL2Ω2+εuhtL2Ω2=u̇0L2Ω2+εu0L2Ω2.

6. 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.

References

1. 1. Ainsworth M. and Oden J. T. A posteriori Error Estimation in Finite Element Analysis. John Wiley, New York, 2000
2. 2. Brenner S. C. and Scott L. R. The Mathematical Theory of Finite Element Methods, second edition. Springer-Verlag, New York, 2002
3. 3. Ciarlet P. G. The finite element method for elliptic problems. North-Holland Pub. Co., Amsterdam-New York-Oxford, 1978
4. 4. Johnson C. Numerical Solution of Partial Differential Equations by the Finite Element Method. Cambridge University Press, New York, 1987
5. 5. Kaltenbacher M. Numerical simulation of mechatronic sensors and actuators: finite elements for computational multiphysics. Springer, Heidelberg, 2015
6. 6. Larson M. The finite element method: theory, implementation, and applications. Springer, Berlin New York, 2013
7. 7. Oden J. T. and Carey G. F. Finite Elements, Mathematical Aspects. Prentice Hall, Englewood Cliffs, 1983
8. 8. Schwab C. pandhpFinite Element Methods. Oxford University Press, New York, 1998
9. 9. Szabo B. and Babu I. s ka. Finite element analysis. Wiley, New York, 1991
10. 10. Hrennikoff A. Solution of problems of elasticity by the framework method. Journal of Applied Mechanics, 8(4):169–175, 1941
11. 11. Courant R. Variational methods for the solution of problems of equilibrium and vibrations. bulletin of the american mathematical society. 49: 1–23. doi: 10.1090. Technical report, 1943
12. 12. Strang G. and Fix G. J. An analysis of the finite element method. 1973

Written By

Mahboub Baccouch

Submitted: June 9th, 2020 Reviewed: December 8th, 2020 Published: February 17th, 2021