Open access peer-reviewed chapter

A Brief Summary of the Finite Element Method for Differential Equations

By Mahboub Baccouch

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

DOI: 10.5772/intechopen.95423

Downloaded: 246

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.

Advertisement

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 ato 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 xiis called a node or nodal point. The length of the interval (called an element) Ii=xi1xiis 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 ube in the space V. For the model problem (1), the solution space is V=C1ab. We wish to construct a finite dimensional space (subspace) VhVbased 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 P1Iiis the space of polynomials of degree 1on Ii. Functions in Vhare 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 vVhis uniquely determined by its nodal values vxi.

  4. Construct a set of basis functions based on the triangulation. Since Vhhas finite dimension, we can find one set of basis functions. A basis for Vhis ϕjj=0N, where ϕjVhare linearly independent. Then Vh=vhxVvhx=j=0Ncjϕjxis 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.

Figure 1.

A continuous piecewise linear functionv.

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=1and ϕixj=0for 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 uby a continuous piecewise linear function uhx. The FE method consists of finding uhVhsuch 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,,cNare unknowns, and choosing v=ϕi,i=1,2,,Nto 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=0and uhxi=j=0Ncjϕjxi=ciϕixi=cifor i=1,2,,N. Thus, we get the following linear system

Figure 2.

A typical hat functionϕion a mesh. Also shown is the half hat functionsϕ0andϕN.

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+fband using ϕixi1=ϕixi+1=0and ϕ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=0for i=2,,Nand

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ϕjxand 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,,cNusing 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ϕidxand β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,jas

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 uC2absuch that

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

where u:Ω¯=abRis the sought solution, qx0is 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 vin 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 ato 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 uH01that satisfies (4). We note that a solution uto (4) is less regular than the solution u(3). Indeed, (4) has only uwhereas (3) contains u. Furthermore, we can easily verify the following:

  1. If uis strong solution (i.e.,solves (3)) then uis also weak solution (i.e.,solves (4)).

  2. Conversely, if uis 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 H01contains many functions and it is therefore just as hard to find a function uH01which 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=bbe a regular partition of ab. Suppose that the length of Ii=xi1xiis hi=xixi1. We define h=maxi=1,2,,Nhito be the mesh size. We wish to construct a subspace VhV=H01. Since Vhhas finite dimension, we can find one set of basis functions ϕjj=1N1for Vh, where ϕjVh,j=1,2,,N1are linearly independent. We remark that Vhis 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 Vhthat satisfies the following properties

  1. The matrix Amust 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. uhmust converge to the solution uof the original problem as h0.

It is natural to obtain an approximation uhto uas follows: Find uhVhsuch that

abuhvdx+abquhvdx=abfvdx,vVh.E5

We call uhthe FE approximation of u. We say that (5) is the Galerkin approximation of (4) and the method used to find uhVhis 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,01has 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 uhwithin a small (finite dimensional) subspace Vh,01=vVh1va=vb=0of H01, consisting of piecewise linear polynomials, where Vh1=vC0abvIiP1Ii.

Let Vh,01be the space of all continuous piecewise linear functions, which vanish at the end points aand b. There are many types of basis functions ϕii=1N1. The simplest ones are the so-called hat functions satisfying ϕixj=δij, where δijis the Kronecker symbol. Note especially that there is no need to construct hat functions ϕ0and ϕNsince any function of Vh,01must vanish at the end points x0=aand xN=b.

The explicit expressions for the hat function ϕixand its derivative ϕixare 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 uhthe FE approximation of u. We say that (6) is the Galerkin approximation of (4) and the method used to find uhVh,01is called Galerkin method.

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

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 cjare real numbers to be determined. We note that the coefficients cj,j=1,2,,N1are the N1nodal values of uhto be determined. Note that the index is only from 1 to N1, because of the zero boundary conditions. We remark that uha=uhb=0and uhxi=ci. So ciis 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,,cN1that 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=c1c2cN1tRN1is the unknown vector, Ais an N1×N1matrix, 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ϕjand bi=fϕi, where auv=abuv+quvdxis a bi-linear and fv=abfvdxis 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,,cN1and 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 âi,j=abϕiϕjdx. For ij>1, we have âi,j=0, since ϕiand ϕjlack overlapping support. However, if i=j, then

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

Furthermore, if j=i+1, then

â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

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

To obtain a˜i,j=abqϕiϕjdxand 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=qxiand fi=fxi. Thus, the matrix A=âi,j+a˜i,jis tridiagonal and has the form

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

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

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=0and

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 uis 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 priorierror 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 solutionuhis 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 Cis 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 uhto 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 uand uhis controlled by the interpolation error uπuin 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 Cis 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 uhit is an a posteriorierror estimate (this yields a computable error bound).

  3. uhuin the v-norm as h=maxihi0. If uuh=0then uuhis constant, but since u0=uh0we also have uuh=0and therefore uh=u.

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

  5. The norm vis 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=αuaand ub=βubfor 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 uC2absuch that

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

where αand βare given constants and fCabis a given function. In this case, the admissible function space H01=v:v2+v2<va=vb=0and the FE space Vh,01defined earlier remain the same. Multiplying (15) by a test function vH01and 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 Vh1and 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 aand b. We also let a=x0<x1<<xN=bbe a uniform partition of the interval ab. Moreover let ϕibe the set of hat basis functions of Vhassociated with the N+1nodes xj,j=0,1,,N, such that ϕixj=δij. The FE approximation of (16) thus reads: Find uhVh1such that uha=α, uhb=β, and

abuhvdx=abfvdx,vVh,01.E17

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

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

Expanding uhas a linear combination of hat functions

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

where the coefficients cj,j=1,2,,N1are the N1nodal values of uhto 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×N1system of equations for cj. In matrix form we write

Ac=b,E20

where Ais a N1×N1matrix, the so-called stiffness matrix, with entries

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

c=c1c2cN1tis a N1vector containing the unknown coefficients cj,j=1,2,,N1, and bis a N1vector, 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 ϕixis 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=hfor i=1,2,,N. Hence the derivative ϕixis 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 ϕiand ϕjlack 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 ito i1we also have

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

Thus the stiffness matrix is

A=1h21001210012010012.

The entries biof the load vector must often be evaluated using quadrature, since they involve the function fwhich 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=α=uaand uhb=β=ub. Therefore, we see that the system matrix Aremains the same, and only the first and last entries of the load vector bneed to be modified because of the definition of the basis functions ϕ0ϕN. An alternative approach is to use all the basis functions ϕ0ϕNto form a larger system of equation, i.e., and N+1×N+1system. 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 Aexcept 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 uC2absuch that

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

where α,βand γare given numbers and fCabis a given function. The boundary condition at x=bis called a Robin boundary condition (combination and uand uis 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 vH01and 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 Vh1is 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 uhto the solution usatisfies

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 uC2absuch that

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

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

V=vC0I:v2+v2<.

Multiplying (26) by a function vVand 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=vC0abvIiP1Iisuch that

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

Implementation: We need to assemble a stiffness matrix Aand a load vector b. Substituting uh=i=0Nciϕiinto (27) and taking v=ϕjfor j=0,1,,Nyields

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+1system of equations for ci. In matrix form we write Ac=b, where c=c0cNtis a N+1vector containing the unknown coefficients ci,i=0,1,,N, Ais a N+1×N+1matrix with entries

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

and bis a N+1vector with entries

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

Let for simplification p=1. Then the matrix Aand 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 P2elements

Let a=x0<x1<<xN=bbe a regular partition of the interval ab. Suppose that the length of Ii=xi1xiis hi=xixi1. Let Pk=px=j=0kcjxjcjRdenotes the vector space of polynomials in one variable and of degree less than or equal to k. The FE method for Lagrange P2elements 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 P2FE 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=0Nis the basis of the shape functions ϕjdefined 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 Vh2and the three quadratic Lagrange P2shape functions on the reference interval 11.

Figure 3.

(left) global shape functions for the spaceVh2. (right) the three quadratic LagrangeP2shape functions on the reference interval11.

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=0Nis the basis of the shape functions ϕjdefined 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,,2N1for the mesh points and ϕj2,j=1,,2N1for the basis of V0,h2. Using these notations, we have:

uh=j=12N1cj2ϕj2x,

where cj2=uhxj2uxj2are the unknowns coefficients. This formulation leads to solve in R2N1a linear system:

Ac=b,

where c=c12c1cN12tR2N1is the unknown vector containing the coefficients cj2,j=1,2,,2N1, Ais an 2N1×2N1matrix with entries

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

and load vector bR2N1has entries

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

Since the shape functions ϕihave a small support, the matrix Ais mostly composed of zeros. However, the main difference with the Lagrange P1FE method, the matrix Ais no longer a tridiagonal matrix.

Computer Implementation: The coefficients of the matrix Acan 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=xi1xion 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 Iito the stiffness matrix and to the mass matrix are given by the 3×3matrices KIiand 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 fis only known by its values at the mesh points xi2,i=0,1,,2Nand thus, we use the decomposition of fin the basis of shape functions ϕi2,i=0,1,,2Nas fx=j=02Nfxj2ϕj2. Each component bi2of 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=xi1xiwill contribute to only three components of indices i1, i12and ias:

bIi=hi304212162124fxi1fxi12fxi..

2.5.2 Nonhomogeneous boundary conditions

Consider the following two-point BVP: find uC2absuch that

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

where αand βare given constants and fCabis a given function.

Multiplying (29) by a function vH01=v:v2+v2<va=vb=0and 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 Vh2and 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 aand b, on a uniform partition a=x0<x1<<xN=bof 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,2Nfor the mesh points and ϕj2,j=0,1,,2N1,2Nfor the basis of Vh2and ϕj2,j=1,,2N1for the basis of V0,h2. Using these notations, we have:

uh=j=02Ncj2ϕj2x,

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

Ac=b,

where c=c12c1cN12tR2N1is the unknown vector containing the coefficients cj2,j=1,2,,2N1, Ais an 2N1×2N1matrix with entries

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

and the load vector bR2N1has 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=0then 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.
Advertisement

3. The FE for elliptic PDEs

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

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

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

3.1 Meshes

Let ΩR2bounded with ∂Ωassumed to be polygonal. A triangulation Thof Ωis a set of triangles Tsuch 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=diamTthe length or the largest edge.

Let Thhave Nnodes and Mtriangles. The data is stored in two matrices. The matrix PR2×Ndescribes the nodes (x1y1,,xNyN)and the matrix KR3×Mdescribes the triangles, i.e.,it describes which nodes (numerated from 1to N) form a triangle Tand how it is orientated:

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

This means that triangle Tiis 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 Tbe a triangle with nodes N1=x1y1, N2=x2y2, and N3=x3y3. We define

P1T=vC0Tvxy=c1+c2x+c3yc1c2c3R.

Now let vi=vNifor i=1,2,3. Note that vP1Tis determined by vii=13. Given viwe compute ciby

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 λjP1Tbe given by the nodal values λjNi=δij, where δijis the Kronecker symbol. This gives us vxy=α1λ1xy+α2λ2xy+α3λ3xy, where αi=vNifor i=1,2,3.We can compute λixyas 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 Thbe a triangulation of Ω, then we let

Vh=vCΩvTP1TTTh.

Functions in Vhare piecewise linear and continuous. We know that vVhis uniquely determined by vNii=12N. We let ϕjNi=δijand let ϕjj=12NVhbe a basis for Vh(hat functions), i.e.,

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

3.3 Interpolation

Given uCTon 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 Cis a generic constant independent of hTand 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 Cis a generic constant independent of hand 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ΩvTP1TTThby PhuVhsuch that

ΩuPhuvhdxdy=0,vhVh.

The problem of finding PhuVhis equivalent to solve the following linear system

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

where ϕii=1Nis a basis of Vh.

Since PhuVhwe can express it as Phu=i=1Nciϕixy, where ciR. Therefore, to find PhuVhwe need to find c1,c2,,cNRsuch 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=c1c2cNtand the entries of the matrix MRN×Nand the vector bRNare 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 usuch that

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

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

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 uVsuch 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 uhVhsuch that

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

where Vh=vVvTP1TTTh.

Implementation: Let a=1and b=g=0. Substituting uh=j=1Ncjϕjinto (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=c1c2cNtRNis the unknown vector and the entries of ARN×N, RRN×N, and bRNare 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 Tis an arbitrary triangle with vertices Ni=xiyiand ϕiare the hat functions i.e.,ϕjNi=δij. Let ϕixy=αi+βix+γiy, for i=1,2,3. Then, we compute αi,βi,γiby

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

In general we have Bαi=eifor 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 Γhoutdenote 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×2by the entries

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

where Eis the length of Eand δijis 1 for i=jand 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, Rand b, we can solve A+Rc=band write uh=j=1Ncjϕj.

3.7 The Dirichlet problem

Consider the following Dirichlet Problem: Find usuch that

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

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

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

So the weak problem reads: Find uVgsuch that

Ωuvdxdy=Ωfvdxdy,vV0.

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

Ωuhvhdxdy=Ωfvhdxdy,vhVh,0.

Assume that we have Nnodes and Jboundary 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 c1RJis known (it contains the values of gin the boundary nodes). We can therefore solve the simplified problem reading: find c0RNJwith A0,0c0=b0A0,gc1.

3.8 The Neumann problem

Consider the following Neumann Problem: Find usuch 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 uVsuch that

Ωuvdxdy∂Ωvgds=Ωfvdxdy,vV.

In order to guarantee solvability, we note that if v=1then 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 uexists, it is only determined up to a constant, since u+cis a solution if uis 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̂his 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 R2or R3and 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=0on Γ1and unon Γ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 Γ1includes its “endpoints”). We can then choose the approximating subspace of V̂as follows:

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

A basis for Vhis 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 Γ4are, 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 wdefined on Ω¯and satisfying wx=gxfor all x∂Ω. We then define v=uwand note that

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

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

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

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

We now describe a method for computing a function wthat 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 giby a function hiwhich differs from giby 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 Γ1and Γ3, while the third term interpolates between the boundary values on Γ2and Γ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 wundoes this transformation. It is straightforward to verify that wsatisfies 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 Γ4are, 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 usatisfying 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 glxto a function h1xsatisfying h10=h1a=0, and similarly for g2, g3, g4and 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 qxysatisfying 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 wis to satisfy the desired Neumann conditions, then wq=hion Γi, i=14, where

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

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

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

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

3.11 Eigenvalue problem

Consider the following Eigenvalue Problem: Find λRand usuch 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 uVsuch that

Ωuvdxdy=λΩuvdxdy,vV.E44

The FE method in 2D is defined as follows: Find λhRand uhVhsuch that

Ωuhvhdxdy=λhΩuhvhdxdy,vhVh,E45

where Vh=vVvTP1TTTh.

Implementation: Substituting uh=j=1Ncjϕjinto (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 usuch that

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

The weak form: Find uV0such that

Ωuvdxdy=Ωfvdxdy,vV0.

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

Ωuhvhdxdy=Ωfvhdxdy,vhVh,0,

where Vh=vVvTP1TTTh. Expressing uh=j=1Ncjϕjand 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×Nand bRNare

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=Ωv2dxdybe the energy norm on V0.

There are two different kinds of error estimates, a prioriestimates, where the error is bounded in terms of the exact solution, and a posteriorierror 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 usuch that

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

We seek a weak solution uin 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 uV0such that

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

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

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

where Vh=vVvTP1