A Brief Summary of the Finite Element Method for Differential Equations

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

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^{\prime }=f\left(x\right),x\in \left[a,b\right],u\left(a\right)=u_0.$$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 $v\left(x\right)\in V_0=\left\{v\in L^2\left[a,b\right]:{\parallel v\parallel }^2+{\parallel v\prime \parallel }^2<\infty ,v\left(a\right)=0\right\}$, where ${\parallel v\parallel }^2={\int }_a^b{v}^2\left(x\right)\mathit{dx}$, integrating from $a$ to $b$, using integration by parts, and applying $v\left(a\right)=0$, to get

   ${\int }_a^b\mathit{\text{fvdx}}={\int }_a^b{u}^{\prime }\mathit{vdx}=-{\int }_a^{b}u{v}^{\prime }\mathit{dx}+u\left(b\right)v\left(b\right)-u\left(a\right)v\left(a\right)=-{\int }_a^{b}u{v}^{\prime }\mathit{dx}+u\left(b\right)v\left(b\right).$

2. Generate a triangulation (also called a mesh) of the computational domain $\left[a,b\right]$. For a one-dimensional problem, a mesh is a set of points in the interval $\left[a,b\right]$, say, $a=x_0\le x_1\le \cdots \le x_N=b$. The point $x_i$ is called a node or nodal point. The length of the interval (called an element) $I_i=\left[x_{i-1},x_i\right]$ is $h_i=x_i-x_{i-1}$. Let $h=\underset{1\le i\le N}{max}{h}_i$ (called a mesh size that measures how fine the partition is). If the mesh is uniformly distributed, then $x_i=a+ih$, $i=0,1,\dots ,N$, where $h=\frac{b-a}{N}$.

3. Define a finite dimensional space over the triangulation: Let the solution $u$ be in the space $V$. For the model problem (1), the solution space is $V=C^1\left[a,b\right]$. We wish to construct a finite dimensional space (subspace) $V_h\subset V$ based on the mesh. When the FE space is a subspace of the solution space, the method is called conforming. It is known that in this case, the FE solution converges to the true solution provided the FE space approximates the given space in some sense [3]. Different finite dimensional spaces will generate different FE solutions.

   Define the FE space as the set of all continuous piecewise linear polynomials ${\left.V_h=\{v:v\right|}_{I_i}\in P^1\left(I_i\right),i=1,2,\dots ,N,v\left(a\right)=0\}$, where $P^1\left(I_i\right)$ is the space of polynomials of degree $\le 1$ on $I_i$. Functions in $V_h$ are linear on each $I_i$, and continuous on the whole interval $\left[a,b\right]$. An example of such a function is shown in Figure 1.

   We remark that any function $v\in V_h$ is uniquely determined by its nodal values $v\left(x_i\right)$.

4. Construct a set of basis functions based on the triangulation. Since $V_h$ has finite dimension, we can find one set of basis functions. A basis for $V_h$ is ${\left\{{\varphi }_j\right\}}_{j=0}^N$, where ${\varphi }_j\in V_h$ are linearly independent. Then $V_h=\left\{v_h\left(x\right)\in V,v_h\left(x\right)={\sum }_{j=0}^N{c}_j{\varphi }_j\left(x\right)\right\}$ is the space spanned by the basis functions ${\left\{{\varphi }_i\right\}}_{i=0}^N$. The simplest finite dimensional space is the piecewise continuous linear function space defined over the triangulation.

Figure 1.

A continuous piecewise linear function $v$.

$$V_h=\left\{v_h\left(x\right)\in V,v_h\left(x\right)\text{is piecewise continuous linear over}\left[\mathrm{a},\mathrm{b}\right]\text{with}{v}_h\left(a\right)=0\right\}.$$

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 ${\varphi }_i\left(x_i\right)=1$ and ${\varphi }_i\left(x_j\right)=0$ for $i\ne j$. The analytic form is (see Figure 2)

$${\varphi }_0\left(x\right)=\left\{\begin{array}{ll}\frac{x_1-x}{h},& x\in I_1,\\ 0,& \text{else},\end{array}\right.,{\varphi }_N\left(x\right)=\left\{\begin{array}{ll}\frac{x-x_{N-1}}{h},& x\in I_N,\\ 0,& \text{else},\end{array}\right.,{\varphi }_i\left(x\right)=\left\{\begin{array}{ll}\frac{x-x_{i-1}}{h},& x\in I_i,\\ \frac{x_{i+1}-x}{h},& x\in I_{i+1},\\ 0,& \text{else}.\end{array}\right.$$

1. Approximate the exact solution $u$ by a continuous piecewise linear function $u_h\left(x\right)$. The FE method consists of finding $u_h\in V_h$ such that

$$-{\int }_a^b{u}_h{v}^{\prime }\mathit{dx}+u_h\left(b\right)v\left(b\right)={\int }_a^b\mathit{\text{fvdx}},\forall v\in V_h.$$

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 $u_h\left(x\right)={\sum }_{j=0}^N{c}_j{\varphi }_j\left(x\right)$, where $c_0,c_1,\dots ,c_N$ are unknowns, and choosing $v={\varphi }_i,i=1,2,\dots ,N$ to get

$$-\sum _{j=0}^N{c}_j{\int }_a^b{\varphi }_j{\varphi }_i^{\prime }\mathit{dx}+c_N{\varphi }_i\left(b\right)={\int }_a^{b}f{\varphi }_i\mathit{dx},i=1,2,\dots ,N,$$

since $u_h\left(b\right)=c_N$. Note that using the hat functions, we have $u_h\left(x_0\right)=0$ and $u_h\left(x_i\right)={\sum }_{j=0}^N{c}_j{\varphi }_j\left(x_i\right)=c_i{\varphi }_i\left(x_i\right)=c_i$ for $i=1,2,\dots ,N$. Thus, we get the following linear system

Figure 2.

A typical hat function ${\varphi }_i$ on a mesh. Also shown is the half hat functions ${\varphi }_0$ and ${\varphi }_N$.

$$-\sum _{j=1}^N{c}_j{\int }_a^b{\varphi }_j{\varphi }_i^{\prime }\mathit{dx}+c_N{\varphi }_i\left(b\right)={\int }_a^b{\varphi }_0{\varphi }_i^{\prime }\mathit{dx},i=1,2,\dots ,N.$$

Finally, we solve the linear system for $c_1,\dots ,c_N$. We note that for $i=1,2,\dots ,N-1$, we have

$${\int }_a^b{\varphi }_i{\varphi }_i^{\prime }\mathit{dx}={\int }_{x_{i-1}}^{x_{i+1}}{\varphi }_i{\varphi }_i^{\prime }\mathit{dx}=\frac{1}{h_i}{\int }_{x_{i-1}}^{x_i}\left(\frac{x-x_{i-1}}{h_i}\right)\mathit{dx}-\frac{1}{h_i}{\int }_{x_i}^{x_{i+1}}\left(\frac{x_{i+1}-x}{h_i}\right)\mathit{dx}=0.$$

However, for $i=N$, we have

$${\int }_a^b{\varphi }_N{\varphi }_N\text{'}\mathit{dx}={\int }_{x_{N-1}}^{x_N}{\varphi }_N{\varphi }_N\text{'}\mathit{dx}={\int }_{x_{N-1}}^{x_N}\left(\frac{x-x_{N-1}}{h_N}\right)\left(\frac{x-x_{N-1}}{h_N}\right)\mathit{dx}=\frac{1}{h_N}{\int }_{x_{N-1}}^{x_N}\frac{x-x_{N-1}}{h_N}\mathit{dx}=\frac{1}{2}.$$

Similarly, for $i=1,2,\dots ,N$, we have

$${\int }_a^b{\varphi }_{i-1}{\varphi }_i^{\prime }\mathit{dx}={\int }_{x_{i-1}}^{x_i}{\varphi }_{i-1}{\varphi }_i^{\prime }\mathit{dx}={\int }_{x_{i-1}}^{x_i}\left(\frac{x_i-x}{h_i}\right)\left(\frac{x-x_{i-1}}{h_i}\right)\mathit{dx}=\frac{1}{h_i}{\int }_{x_{i-1}}^{x_i}\frac{x_i-x}{h_i}\mathit{dx}=\frac{1}{2},$$

$${\int }_a^b{\varphi }_{i+1}{\varphi }_i^{\prime }\mathit{dx}={\int }_{x_i}^{x_{i+1}}{\varphi }_{i+1}{\varphi }_i^{\prime }\mathit{dx}={\int }_{x_i}^{x_{i+1}}\left(\frac{x-x_i}{h_{i+1}}\right)\left(\frac{x_{i+1}-x}{h_{i+1}}\right)\mathit{dx}=-\frac{1}{h_{i+1}}{\int }_{x_i}^{x_{i+1}}\frac{x-x_i}{h_{i+1}}\mathit{dx}=-\frac{1}{2}.$$

We next calculate ${\int }_a^{b}f{\varphi }_i\mathit{dx}$. 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 ${\int }_a^{b}f\left(x\right)\mathit{dx}\approx \frac{b-a}{2}\left(f\left(a\right)+f\left(b\right)\right)$ and using ${\varphi }_i\left(x_{i-1}\right)={\varphi }_i\left(x_{i+1}\right)=0$ and ${\varphi }_i\left(x_i\right)=1$, we get

$${\int }_a^{b}f{\varphi }_i\mathit{dx}={\int }_{x_{i-1}}^{x_i}f{\varphi }_i\mathit{dx}+{\int }_{x_i}^{x_{i+1}}f{\varphi }_i\mathit{dx}\approx \frac{h_i+h_{i+1}}{2}f\left(x_i\right),i=1,2,\dots ,N-1,$$

$${\int }_a^{b}f\left(x\right){\varphi }_N\mathit{dx}={\int }_{x_{N-1}}^{x_N}f\left(x\right){\varphi }_N\mathit{dx}\approx \frac{h_N}{2}\left(f\left(x_{N-1}\right){\varphi }_N\left(x_{N-1}\right)+f\left(x_N\right){\varphi }_N\left(x_N\right)\right)=\frac{h_N}{2}f\left(x_N\right).$$

Thus, we obtain the following linear system of equations

$$\left[\begin{array}{ccccc}0& \frac{1}{2}& 0& \cdots & 0\\ -\frac{1}{2}& 0& \frac{1}{2}& \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & -\frac{1}{2}& 0& \frac{1}{2}\\ 0& \cdots & 0& -\frac{1}{2}& \frac{1}{2}\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_{N-1}\\ c_N\\ \end{array}\right]=\left[\begin{array}{c}\frac{h_1+h_2}{2}f\left(x_1\right)\\ \frac{h_2+h_3}{2}f\left(x_2\right)\\ ⋮\\ \frac{h_{N-1}+h_N}{2}f\left(x_{N-1}\right)\\ \frac{h_N}{2}f\left(x_N\right)\end{array}\right].$$

The determinant of the above matrix is $\frac{1}{2^N}$. Thus, the system has a unique solution $c_1,c_2,\dots ,c_N$.

Remark 2.1Suppose that$u\left(a\right)=u_0$, then we let$u_h\left(x\right)={\sum }_{j=0}^N{c}_j{\varphi }_j\left(x\right)$. Since$u_0=u_h\left(x_0\right)={\sum }_{j=1}^N{c}_j{\varphi }_j\left(x_0\right)=c_0{\varphi }_0\left(x_0\right)=c_0$, we only need to find$c_1,c_2,\dots ,c_N$. Choosing$v={\varphi }_i,i=1,2,\dots ,N$, we get the following linear system

$$-\sum _{j=1}^N{c}_j{\int }_a^b{\varphi }_j{\varphi }_i^{\prime }\mathit{dx}+c_N{\varphi }_i\left(b\right)={\int }_a^{b}f{\varphi }_i\mathit{dx}+u_0{\int }_a^b{\varphi }_0{\varphi }_i^{\prime }\mathit{dx},i=1,2,\dots ,N.$$

Finally, we solve the linear system for $c_1,\dots ,c_N$. We note that ${\int }_a^b{\varphi }_0{\varphi }_i^{\prime }\mathit{dx}=0$ for $i=2,\dots ,N$ and

$${\int }_a^b{\varphi }_0{\varphi }_1\text{'}\mathit{dx}={\int }_{x_0}^{x_1}\left(\frac{x_1-x}{h_1}\right){\left(\frac{x-x_0}{h_1}\right)}^{\prime }\mathit{dx}=\frac{1}{h_1}{\int }_{x_0}^{x_1}\frac{x_1-x}{h_1}\mathit{dx}=\frac{1}{2}.$$

Following the same steps used for the case $u\left(a\right)=0$, we obtain the following linear system of equations

$$\left[\begin{array}{ccccc}0& \frac{1}{2}& 0& \cdots & 0\\ -\frac{1}{2}& 0& \frac{1}{2}& \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & -\frac{1}{2}& 0& \frac{1}{2}\\ 0& \cdots & 0& -\frac{1}{2}& \frac{1}{2}\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_{N-1}\\ c_N\end{array}\right]=\left[\begin{array}{c}\frac{h_1+h_2}{2}f\left(x_1\right)+\frac{u_0}{2}\\ \frac{h_2+h_3}{2}f\left(x_2\right)\\ ⋮\\ \frac{h_{N-1}+h_N}{2}f\left(x_{N-1}\right)\\ \frac{h_N}{2}f\left(x_N\right)\end{array}\right].$$

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^{\prime }=f\left(x,u\right),x\in \left[a,b\right],u\left(a\right)=u_0.$$E2

The FE method consists of finding ${\left.u_h\in V_h=\{v:v\right|}_{I_i}\in P^1\left(I_i\right),i=1,2,\dots ,N,v\left(a\right)=0\}$, such that

$$u_h\left(b\right)v\left(b\right)-{\int }_a^b{u}_h{v}^{\prime }\mathit{dx}={\int }_a^{b}f\left(x,u_h\right)\mathit{vdx},\forall v\in V_h.$$

Writing $u_h\left(x\right)={\sum }_{j=0}^N{c}_j{\varphi }_j\left(x\right)$ and choosing $v={\varphi }_i,i=1,2,\dots ,N$, we get

$$c_N\left({\varphi }_i-{\int }_a^b{\varphi }_N{\varphi }_i^{\prime }\mathit{dx}\right)-\sum _{j=0}^{N-1}{c}_j{\int }_a^b{\varphi }_j{\varphi }_i^{\prime }\mathit{dx}-{\int }_a^{b}f\left(x,\sum _{j=0}^N{c}_j{\varphi }_j\right){\varphi }_i\mathit{dx}=0,i=1,2,\dots ,N,$$

where $u_h\left(x_0\right)=c_0=u_0$. Finally, we solve the nonlinear system for $c_1,c_2,\dots ,c_N$ using $e.g.,$ Newton’s method for systems of nonlinear equations. The system can be written as $F_i\left(c_1,c_2,\dots ,c_N\right)=0,i=1,2,\dots ,N$, where

$$F_i=c_N\left({\varphi }_i-{\int }_a^b{\varphi }_N{\varphi }_i^{\prime }\mathit{dx}\right)-\sum _{j=0}^{N-1}{c}_j{\int }_a^b{\varphi }_j{\varphi }_i^{\prime }\mathit{dx}-{\int }_a^{b}f\left(x,\sum _{j=0}^N{c}_j{\varphi }_j\right){\varphi }_i\mathit{dx},i=1,2,\dots ,N.$$

Let ${\alpha }_i={\sum }_{j=0}^N{c}_j{\int }_a^b{\varphi }_j{\varphi }_i^{\prime }\mathit{dx}$ and ${\beta }_i={\int }_a^{b}f\left(x,{\sum }_{j=0}^N{c}_j{\varphi }_j\right){\varphi }_i\mathit{dx}$. Then, for $i=1,2,\dots ,N-1$,

$$\begin{array}{cc}{\alpha }_i& =c_{i-1}{\int }_{x_{i-1}}^{x_i}{\varphi }_{i-1}{\varphi }_i^{\prime }\mathit{dx}+c_i\left({\int }_{x_{i-1}}^{x_i}{\varphi }_i{\varphi }_i^{\prime }\mathit{dx}+{\int }_{x_i}^{x_{i+1}}{\varphi }_i{\varphi }_i^{\prime }\mathit{dx}\right)+c_{i+1}{\int }_{x_i}^{x_{i+1}}{\varphi }_{i+1}{\varphi }_i^{\prime }\mathit{dx}\\ & =c_{i-1}{\int }_{x_{i-1}}^{x_i}\frac{x_i-x}{h_i^2}\mathit{dx}+c_i\left({\int }_{x_{i-1}}^{x_i}\frac{x-x_{i-1}}{h_i^2}\mathit{dx}-{\int }_{x_i}^{x_{i+1}}\frac{x_{i+1}-x}{h_{i+1}^2}\mathit{dx}\right)-c_{i+1}{\int }_{x_i}^{x_{i+1}}\frac{x-x_i}{h_{i+1}^2}\mathit{dx}\\ & =\frac{1}{2}{c}_{i-1}+c_i\left(\frac{1}{2}-\frac{1}{2}\right)-\frac{1}{2}{c}_{i+1}=\frac{1}{2}{c}_{i-1}-\frac{1}{2}{c}_{i+1},\\ {\alpha }_N& =c_{N-1}{\int }_{x_{N-1}}^{x_N}{\varphi }_{N-1}{\varphi }_N\text{'}\mathit{dx}+c_N{\int }_{x_{N-1}}^{x_N}{\varphi }_N{\varphi }_N\text{'}\mathit{dx}\\ & =c_{N-1}{\int }_{x_{N-1}}^{x_N}\frac{x_N-x}{h_N^2}\mathit{dx}+c_N{\int }_{x_{N-1}}^{x_N}\frac{x-x_{N-1}}{h_N^2}\mathit{dx}=\frac{1}{2}{c}_{N-1}+\frac{1}{2}{c}_N.\end{array}$$

Similarly,

$${\beta }_i={\int }_{x_{i-1}}^{x_{i+1}}f\left(x,\sum _{j=0}^N{c}_j{\varphi }_j\right){\varphi }_i\mathit{dx}={\int }_{x_{i-1}}^{x_i}f\left(x,\sum _{j=0}^N{c}_j{\varphi }_j\right){\varphi }_i\mathit{dx}+{\int }_{x_i}^{x_{i+1}}f\left(x,\sum _{j=0}^N{c}_j{\varphi }_j\right){\varphi }_i\mathit{dx}.$$

Using Simpson’s Rule ${\int }_a^{b}f\left(x\right)\mathit{dx}\approx \frac{b-a}{6}\left(f\left(a\right)+4f\left(\frac{a+b}{2}\right)+f\left(b\right)\right)$, and using ${\varphi }_i\left(x_{i-1}\right)={\varphi }_i\left(x_{i+1}\right)=0$, ${\varphi }_i\left(x_i\right)=1$, ${\sum }_{j=0}^N{c}_j{\varphi }_j\left(x_{i-1}+\frac{h_i}{2}\right)=\frac{c_{i-1}+c_i}{2}$, ${\varphi }_i\left(x_{i-1}+\frac{h_i}{2}\right)=\frac{1}{2}$, ${\sum }_{j=0}^N{c}_j{\varphi }_j\left(x_i\right)=c_i$, we have, for $i=1,2,\dots ,N-1$,

$${\beta }_i\approx \frac{h_i}{3}f\left(x_{i-1}+\frac{h_i}{2},\frac{c_{i-1}+c_i}{2}\right)+\frac{h_i+h_{i+1}}{6}f\left(x_i,c_i\right)+\frac{h_{i+1}}{3}f\left(x_i+\frac{h_{i+1}}{2},\frac{c_i+c_{i+1}}{2}\right).$$

However, for $i=N$, we have

$${\beta }_N\approx \frac{h_N}{6}\left(2f\left(x_{N-1}+\frac{h}{2},\frac{c_{N-1}+c_N}{2}\right)+f\left(x_N,c_N\right)\right).$$

Next, we compute the Jacobian matrix with entries

$$J_{i,j}=\frac{\partial F_i}{\partial c_j}={\int }_a^b{\varphi }_j{\varphi }_i^{\prime }\mathit{dx}-{\int }_a^b{f}_u\left(x,\sum _{j=0}^N{c}_j{\varphi }_j\right){\varphi }_j{\varphi }_i\mathit{dx}=a_{i,j}-b_{i,j},i=1,2,\dots ,N.$$

We already computed the entries $a_{i,j}$ as

$$\begin{gathered} a_{i,i-1}={\int }_a^b{\varphi }_{i-1}{\varphi }_i^{\prime }\mathit{dx}=\frac{1}{2},a_{i,i}={\int }_a^b{\varphi }_i{\varphi }_i^{\prime }\mathit{dx}=0,i=1,2,\dots ,N-1, \\ {a}_{N,N}={\int }_a^b{\varphi }_N{\varphi }_N\text{'}\mathit{dx}=\frac{1}{2},a_{i,i+1}={\int }_a^b{\varphi }_{i+1}{\varphi }_i^{\prime }\mathit{dx}=-\frac{1}{2}. \end{gathered}$$

Using Simpson’s Rule, we get

$$\begin{array}{cc}{b}_{i,i-1}& ={\int }_{x_{i-1}}^{x_i}{\varphi }_{i-1}{\varphi }_i{f}_u\left(x,\sum _{j=0}^N{c}_j{\varphi }_j\right)\mathit{dx}\approx \frac{h_i}{6}{f}_u\left(x_{i-1}+\frac{h_i}{2},\frac{c_{i-1}+c_i}{2}\right),\\ b_{i,i+1}& ={\int }_{x_i}^{x_{i+1}}{\varphi }_{i+1}{\varphi }_i{f}_u\left(x,\sum _{j=0}^N{c}_j{\varphi }_j\right)\mathit{dx}\approx \frac{h_{i+1}}{6}{f}_u\left(x_i+\frac{h_{i+1}}{2},\frac{c_i+c_{i+1}}{2}\right),\\ b_{i,i}& ={\int }_{x_{i-1}}^{x_i}{\varphi }_i^2{f}_u\left(x,\sum _{j=0}^N{c}_j{\varphi }_j\right)\mathit{dx}+{\int }_{x_i}^{x_{i+1}}{\varphi }_i^2{f}_u\left(x,\sum _{j=0}^N{c}_j{\varphi }_j\right)\mathit{dx}\\ & \approx \frac{h_i}{6}{f}_u\left(x_{i-1}+\frac{h_i}{2},\frac{c_{i-1}+c_i}{2}\right)+\frac{h_i+h_{i+1}}{6}{f}_u\left(x_i,c_i\right)+\frac{h_{i+1}}{6}{f}_u\left(x_i+\frac{h_{i+1}}{2},\frac{c_i+c_{i+1}}{2}\right),\\ b_{N,N}& ={\int }_{x_{N-1}}^{x_N}{\varphi }_N^2{f}_u\left(x,\sum _{j=0}^N{c}_j{\varphi }_j\right)\mathit{dx}\approx \frac{h_N}{6}\left(f_u\left(x_{N-1}+\frac{h}{2},\frac{c_{N-1}+c_N}{2}\right)+f_u\left(x_N,c_N\right)\right).\end{array}$$

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 $u\in C^2\left[a,b\right]$ such that

$$-u^{″}+q\left(x\right)u=f\left(x\right),x\in \mathrm{\Omega }=\left(a,b\right),u\left(a\right)=u\left(b\right)=0,$$E3

where $u:\overline{\mathrm{\Omega }}=\left[a,b\right]\to \mathbb{R}$ is the sought solution, $q\left(x\right)\ge 0$ is a continuous function on $\left[a,b\right]$, and $f\in L^2\left[a,b\right]$. Under these assumptions, (3) has a unique solution $u\in C^2\left[a,b\right]$. For general $q\left(x\right)$, 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: ${\int }_a^b{u}^{\prime }{v}^{\prime }\mathit{dx}+{\int }_a^b\mathit{\text{quvdx}}={\int }_a^b\mathit{\text{fvdx}}$, for any test function $v$ in the Sobolev space $H_0^1\left[a,b\right]=\left\{v\in L^2\left[a,b\right]:{‖v‖}^2+{‖v^{\prime }‖}^2<\infty ,v\left(a\right)=v\left(b\right)=0\right\}$, where ${‖v‖}^2={\int }_a^b{v}^2\left(x\right)\mathit{dx}$. 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: ${min}_{v\left(x\right)\in H_0^1\left[a,b\right]}{\int }_a^b\left(\frac{1}{2}{\left(v^{\prime }\right)}^2+\frac{1}{2}{\mathit{qv}}^2-\mathit{fv}\right)\mathit{dx}$. 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 that$u^{″}$exists and continuous on$\left[a,b\right]$. 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 $H_0^1=\left\{v\in L^2\left(a,b\right):{‖v‖}^2+{‖v^{\prime }‖}^2<\infty ,v\left(a\right)=v\left(b\right)=0\right\}$. Multiplying (3) by a test function $v\in H_0^1$, integrating from $a$ to $b$, and using integration by parts, we get

$${\int }_a^b\mathit{\text{fvdx}}={\int }_a^b-u^{″}\mathit{vdx}+{\int }_a^b\mathit{\text{quvdx}}={\int }_a^b{u}^{\prime }{v}^{\prime }\mathit{dx}+{\int }_a^b\mathit{\text{quvdx}},$$

since $v\left(a\right)=v\left(b\right)=0$. Hence, the weak (or variational) form of (3) reads: Find $u\in H_0^1$, such that

$${\int }_a^b{u}^{\prime }{v}^{\prime }\mathit{dx}+{\int }_a^b\mathit{\text{quvdx}}={\int }_a^b\mathit{\text{fvdx}},\forall v\in H_0^1.$$E4

We want to find $u\in H_0^1$ that satisfies (4). We note that a solution $u$ to (4) is less regular than the solution $u$ (3). Indeed, (4) has only $u^{\prime }$ 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 $u\in C^2\left[a,b\right]$, 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‖}_{\mathrm{\Omega }}$, where $\mathrm{\Omega }=\left[a,b\right]$.

The FE formulation: The FE method is based on the variational form (4). We note that the space $H_0^1$ contains many functions and it is therefore just as hard to find a function $u\in H_0^1$ 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=x_0<x_1<\cdots <x_N=b$ be a regular partition of $\left[a,b\right]$. Suppose that the length of $I_i=\left[x_{i-1},x_i\right]$ is $h_i=x_i-x_{i-1}$. We define $h=\underset{i=1,2,\dots ,N}{max}{h}_i$ to be the mesh size. We wish to construct a subspace $V_h\subset V=H_0^1$. Since $V_h$ has finite dimension, we can find one set of basis functions ${\left\{{\varphi }_j\right\}}_{j=1}^{N-1}$ for $V_h$, where ${\varphi }_j\in V_h,j=1,2,\dots ,N-1$ are linearly independent. We remark that $V_h$ is the space spanned by the basis functions $i.e.,$$V_h=\left\{v_h\left(x\right),v_h\left(x\right)={\sum }_{j=1}^{N-1}{c}_j{\varphi }_j\left(x\right)\right\}$. The FE method consists of choosing a basis for the subspace $V_h$ that satisfies the following properties

1. The matrix $A$ must be sparse ($e.g.$ traditional or banded matrix). In this case, iterative methods for solving linear systems can be adapted to obtain an efficient solution.

2. $u_h$ must converge to the solution $u$ of the original problem as $h\to 0$.

It is natural to obtain an approximation $u_h$ to $u$ as follows: Find $u_h\in V_h$ such that

$${\int }_a^b{u}_h^{\prime }{v}^{\prime }\mathit{dx}+{\int }_a^b{\mathit{qu}}_h\mathit{vdx}={\int }_a^b\mathit{\text{fvdx}},\forall v\in V_h.$$E5

We call $u_h$ the FE approximation of $u$. We say that (5) is the Galerkin approximation of (4) and the method used to find $u_h\in V_h$ is called Galerkin method.

FE approximation using Lagrange${\mathbb{P}}_1$elements: The simplest finite dimensional space is the piecewise continuous linear function space defined over the triangulation

$$V_{h,0}^1=\left\{v_h\in V,v_h\text{is piecewise continuous linear over}\left[\mathrm{a},\mathrm{b}\right]\text{with}{v}_h\left(a\right)=v_h\left(b\right)=0\right\}.$$

It is easy to show that $V_{h,0}^1$ has a finite dimension even although there are infinite number of elements in $V_{h,0}^1$. The approximation of the FE method is therefore to look for an approximation $u_h$ within a small (finite dimensional) subspace $V_{h,0}^1=\left\{v\in V_h^1|v\left(a\right)=v\left(b\right)=0\right\}$ of $H_0^1$, consisting of piecewise linear polynomials, where $V_h^1=\left\{v\in C^0\left[a,b\right]{\left|v\right|}_{I_i}\in P^1\left(I_i\right)\right\}$.

Let $V_{h,0}^1$ 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 ${\left\{{\varphi }_i\right\}}_{i=1}^{N-1}$. The simplest ones are the so-called hat functions satisfying ${\varphi }_i\left(x_j\right)={\delta }_{\mathit{ij}}$, where ${\delta }_{\mathit{ij}}$ is the Kronecker symbol. Note especially that there is no need to construct hat functions ${\varphi }_0$ and ${\varphi }_N$ since any function of $V_{h,0}^1$ must vanish at the end points $x_0=a$ and $x_N=b$.

The explicit expressions for the hat function ${\varphi }_i\left(x\right)$ and its derivative ${\varphi }_i^{\prime }\left(x\right)$ are given by

$${\varphi }_i\left(x\right)=\left\{\begin{array}{ll}0,& a\le x\le x_{i-1},\\ \frac{x-x_{i-1}}{h_i},& x_{i-1}\le x\le x_i,\\ \frac{x_{i+1}-x}{h_{i+1}},& x_i\le x\le x_{i+1},\\ 0,& x_{i+1}\le x\le b,\end{array}\right.,{\varphi }_i^{\prime }\left(x\right)=\left\{\begin{array}{ll}0,& a<x<x_{i-1},\\ \frac{1}{h_i},& x_{i-1}<x<x_i,\\ -\frac{}{h_{i+1}},& x_i<x<x_{i+1},\\ 0,& x_{i+1}<x<b,\end{array}\right.,$$

for $i=1,2,\dots ,N-1$. The FE approximation of (4) thus reads: Find $u\in V_{h,0}^1$, such that

$${\int }_a^b{u}_h^{\prime }{v}^{\prime }\mathit{dx}+{\int }_a^b{\mathit{qu}}_h\mathit{vdx}={\int }_a^b\mathit{\text{fvdx}},\forall v\in V_{h,0}^1.$$E6

We call $u_h$ the FE approximation of $u$. We say that (6) is the Galerkin approximation of (4) and the method used to find $u_h\in V_{h,0}^1$ is called Galerkin method.

It can be shown that (6) is equivalent to the $N-1$ equations

$${\int }_a^b{u}_h^{\prime }{\varphi }_i^{\prime }\mathit{dx}+{\int }_a^b{\mathit{qu}}_h{\varphi }_i\mathit{dx}={\int }_a^{b}f{\varphi }_i\mathit{dx},i=1,2,\dots ,N-1.$$E7

Derivation of the discrete system: Since $u_h\in V_{h,0}^1$, we can express it as a linear combination of hat functions $i.e.,$

$$u_h=\sum _{j=1}^{N-1}{c}_j{\varphi }_j\left(x\right),$$E8

where $c_j$ are real numbers to be determined. We note that the coefficients $c_j,j=1,2,\dots ,N-1$ are the $N-1$ nodal values of $u_h$ to be determined. Note that the index is only from 1 to $N-1$, because of the zero boundary conditions. We remark that $u_h\left(a\right)=u_h\left(b\right)=0$ and $u_h\left(x_i\right)=c_i$. So $c_i$ is an approximate solution to the exact solution at $x=x_i$.

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

Substituting (8) into (7) yields

$$\sum _{j=1}^{N-1}{c}_j\left({\int }_a^b{\varphi }_i^{\prime }{\varphi }_j^{\prime }\mathit{dx}+{\int }_a^{b}q{\varphi }_i{\varphi }_j\mathit{dx}\right)={\int }_a^{b}f{\varphi }_i\mathit{dx},i=1,2,\dots ,N-1.$$E9

The problem (7) is now equivalent to the following: Find the real numbers $c_1,c_2,\dots ,c_{N-1}$ that satisfy the linear system (9).

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

$$A\mathbf{c}=\mathbf{b},$$E10

where $\mathbf{c}={\left[c_1,c_2,\dots ,c_{N-1}\right]}^t\in {\mathbb{R}}^{N-1}$ is the unknown vector, $A$ is an $\left(N-1\right)\times \left(N-1\right)$ matrix, the so-called stiffness matrix when $q=0$, with entries

$$a_{\mathit{ij}}={\int }_a^b\left({\varphi }_i^{\prime }{\varphi }_j^{\prime }+q{\varphi }_i{\varphi }_j\right)\mathit{dx},i,j=1,2,\dots ,N-1,$$E11

and $\mathbf{b}\in {\mathbb{R}}^{N-1}$, the so-called load vector, has entries

$$b_i={\int }_a^{b}f{\varphi }_i\mathit{dx},i=1,2,\dots ,N-1.$$E12

To obtain the approximate solution we need to solve the linear system for the unknown vector $\mathbf{c}$. We note that $a_{\mathit{ij}}=a\left({\varphi }_i,{\varphi }_j\right)$ and $b_i=\left(f,{\varphi }_i\right)$, where $a\left(u,v\right)={\int }_a^b\left(u^{\prime }{v}^{\prime }+\mathit{quv}\right)\mathit{dx}$ is a bi-linear and $\left(f,v\right)={\int }_a^b\mathit{\text{fvdx}}$ 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

$$\underset{v\left(x\right)\in H_0^1\left[a,b\right]}{min}F\left(v\right),\text{where}F\left(v\right)={\int }_a^b\left(\frac{1}{2}{\left(v^{\prime }\right)}^2+\frac{1}{2}{\mathit{qv}}^2-\mathit{fv}\right)\mathit{dx}.$$

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

$$F\left(u_h\right)={\int }_a^b\left(\frac{1}{2}{\left(\sum _{j=1}^{N-1}{c}_j{\varphi }_j^{\prime }\left(x\right)\right)}^2+\frac{1}{2}q{\left(\sum _{j=1}^{N-1}{c}_j{\varphi }_j\left(x\right)\right)}^2-f\left(\sum _{j=1}^{N-1}{c}_j{\varphi }_j\left(x\right)\right)\right)\mathit{dx},$$

which is a multi-variable function of $c_1,c_2,\dots ,c_{N-1}$ and can be written as $F\left(u_h\right)=F\left(c_1,c_2,\dots ,c_{N-1}\right)$. The necessary conditions for a global minimum are $\frac{\partial F}{\partial c_i}=0$, $i=1,2,\dots ,N-1$. Taking the partial derivatives directly with respect to $c_i$, we get

$${\int }_a^b\left({\varphi }_i^{\prime }\left(x\right)\sum _{j=1}^{N-1}{c}_j{\varphi }_j^{\prime }\left(x\right)+q{\varphi }_i\left(x\right)\sum _{j=1}^{N-1}{c}_j{\varphi }_j\left(x\right)-f{\varphi }_i\left(x\right)\right)\mathit{dx}=0,i=1,2,\dots ,N-1.$$

Exchange the order of integration and the summation, we get

$$\sum _{j=1}^{N-1}{c}_j{\int }_a^b\left({\varphi }_i^{\prime }\left(x\right){\varphi }_j^{\prime }\left(x\right)+q{\varphi }_i\left(x\right){\varphi }_j\left(x\right)\right)\mathit{dx}={\int }_a^{b}f{\varphi }_i\left(x\right)\mathit{dx}=0,i=1,2,\dots ,N-1,$$

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 ${\widehat{a}}_{i,j}={\int }_a^b{\varphi }_i^{\prime }{\varphi }_j^{\prime }\mathit{dx}$. For $\mid i-j\mid >1$, we have ${\widehat{a}}_{i,j}=0$, since ${\varphi }_i$ and ${\varphi }_j$ lack overlapping support. However, if $i=j$, then

$${\widehat{a}}_{i,i}={\int }_a^b{\left({\varphi }_i^{\prime }\right)}^2\mathit{dx}={\int }_{x_{i-1}}^{x_i}{\left(\frac{1}{h_i}\right)}^2\mathit{dx}+{\int }_{x_i}^{x_{i+1}}{\left(-\frac{1}{h_{i+1}}\right)}^2\mathit{dx}=\frac{1}{h_i}+\frac{1}{h_{i+1}},i,j=1,2,\dots ,N-1.$$

Furthermore, if $j=i+1$, then

$${\widehat{a}}_{i,i+1}={\int }_a^b{\varphi }_i^{\prime }{\varphi }_{i+1}^{\text{'}}\mathit{dx}={\int }_{x_i}^{x_{i+1}}\left(-\frac{1}{h_{i+1}}\right)\left(\frac{1}{h_{i+1}}\right)\mathit{dx}=-\frac{1}{h_{i+1}},i,j=1,2,\dots ,N-2.$$E13

By symmetry, we also have

$${\widehat{a}}_{i+1,i}={\int }_a^b{\varphi }_{i+1}^{\prime }{\varphi }_i^{\prime }\mathit{dx}=-\frac{1}{h_{i+1}},i,j=1,2,\dots ,N-2.$$

To obtain ${\tilde{a}}_{i,j}={\int }_a^{b}q{\varphi }_i{\varphi }_j\mathit{dx}$ and $b_i={\int }_a^{b}f{\varphi }_i\mathit{dx}$, we use the composite trapezoidal rule

$${\int }_a^{b}f\left(x\right)\mathit{dx}=\sum _{i=1}^N{\int }_{x_{i-1}}^{x_i}f\left(x\right)\mathit{dx}\approx \frac{1}{2}\left[h_{1}f\left(x_0\right)+\sum _{i=1}^{N-1}\left(h_i+h_{i+1}\right)f\left(x_i\right)+h_{N}f\left(x_N\right)\right].$$

So, we can easily verify that

$${\tilde{a}}_{i,j}={\int }_a^{b}q{\varphi }_i{\varphi }_j\mathit{dx}\approx \left\{\begin{array}{ll}\frac{q_i}{2}\left(h_i+h_{i+1}\right),& i=j\\ 0,& i\ne j\end{array}\right.,b_i={\int }_a^{b}f{\varphi }_i\mathit{dx}\approx \frac{1}{2}\left(h_i+h_{i+1}\right)f_i,$$

where $q_i=q\left(x_i\right)$ and $f_i=f\left(x_i\right)$. Thus, the matrix $A=\left({\widehat{a}}_{i,j}+{\tilde{a}}_{i,j}\right)$ is tridiagonal and has the form

$$A=\left[\begin{array}{ccccc}\frac{1}{h_1}+\frac{1}{h_2}+\frac{q_1}{2}\left(h_1+h_2\right)& -\frac{1}{h_2}& 0& \cdots & 0\\ -\frac{1}{h_2}& \frac{1}{h_2}+\frac{1}{h_3}+\frac{q_2}{2}\left(h_2+h_3\right)& -\frac{1}{h_3}& \ddots & 0\\ 0& -\frac{1}{h_3}& \ddots & \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & -\frac{1}{h_{N-1}}\\ 0& \cdots & 0& -\frac{1}{h_{N-1}}& \frac{1}{h_{N-1}}+\frac{1}{h_N}+\frac{q_{N-1}}{2}\left(h_{N-1}+h_N\right)\end{array}\right].$$

Finally, we obtain the following system: $c_0=c_N=0$ and

$$-\frac{1}{h_i}{c}_{i-1}+\frac{1}{h_i+h_{i+1}}{c}_i-\frac{1}{h_{i+1}}{c}_{i+1}+\frac{q_i\left(h_i+h_{i+1}\right)}{2}{c}_i=\frac{1}{2}\left(h_i+h_{i+1}\right)f_i,i=1,2,\dots ,N-1,$$

Remark 2.2Suppose that the partition is uniform$i.e.,$$h_i=h=\frac{b-a}{N}$for all$i=1,2,\dots ,N$. Then the stiffness matrix$A$and the load vector$\mathbf{b}$have the form:

$$A=\left[\begin{array}{ccccc}\frac{2}{h}+{\mathit{hq}}_1& -\frac{1}{h}& 0& \cdots & 0\\ -\frac{1}{h}& \frac{2}{h}+{\mathit{hq}}_2& -\frac{1}{h}& \ddots & 0\\ 0& -\frac{1}{h}& \ddots & \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & -\frac{1}{h_{N-1}}\\ 0& \cdots & 0& -\frac{1}{h}& \frac{2}{h}+{\mathit{hq}}_{N-1}\end{array}\right],\mathbf{b}=h\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\\ ⋮\\ f_{N-1}\end{array}\right].$$

Finally, we obtain the following system: $c_0=c_N=0$ and

$$\frac{-c_{j-1}+2c_j-c_{j+1}}{h}+{\mathit{hq}}_i{c}_j={\mathit{hf}}_i\Rightarrow -\frac{c_{j-1}-2c_j+c_{j+1}}{h^2}+q_i{c}_j=f_i,i=1,2,\dots ,N-1,$$

which is the same system obtained using the finite difference method, where $u^{″}$ is approximated using the second-order midpoint formula $u^{″}\left(x_j\right)\approx \frac{u\left(x_{j-1}\right)-2u\left(x_j\right)+u\left(x_{j+1}\right)}{h^2}$. 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 matrix$A$with entries$a_{i,j}={\int }_a^b{\varphi }_i^{\prime }{\varphi }_j^{\prime }\mathit{dx}$is symmetric positive definite$i.e.,$$a_{i,j}=a_{j,i}$and

$${\mathbf{x}}^{t}A\mathbf{x}=\sum _{i,j=1}^{N-1}{x}_i{a}_{i,j}{x}_j>0,\text{for all nonzero}\mathbf{x}={\left[x_1,\dots ,x_{N-1}\right]}^t\in {\mathbb{R}}^{N-1}.$$

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

Next, we state a general convergence result for the Galerkin method. We first define the following norm and semi-norm: For $v\in H_0^1$, we define

$$‖v‖={\left({\int }_a^b{v}^2\left(x\right)\mathit{dx}\right)}^{1/2},{\left|v\right|}_1=‖v^{\prime }‖={\left({\int }_a^b{\left(v^{\prime }\left(x\right)\right)}^2\mathit{dx}\right)}^{1/2}.$$

Theorem 2.3Suppose that$q\left(x\right)\ge 0$, $\forall x\in \left[a,b\right]$. Let$u$be the solution to (4) and$u_h$be the solution to (6). Then there exists a constant$C$such that

$$‖{\left(u-u_h\right)}^{\prime }‖\le C‖{\left(u-v_h\right)}^{\prime }‖,\forall v_h\in V_{h,0}^1,$$E14

where $C$ is given by $C=1+{max}_{x\in \left[a,b\right]}\mid q\left(x\right)\mid$, which is independent of the choice of $V_{h,0}^1$.

Remark 2.3From (14), taking the minimum over all$v_h\in V_{h,0}^1$, we get$‖{\left(u-u_h\right)}^{\prime }‖\le C\underset{v_h\in V_{h,0}}{min}‖{\left(u-v_h\right)}^{\prime }‖$. Thus,${\left|u-u_h\right|}_1\le C{min}_{v_h\in V_{h,0}}{\left|u-v_h\right|}_1$, where$C=1+{max}_{x\in \left[a,b\right]}\mid q\left(x\right)\mid$.

Next, we study the convergence of $u_h$ to $u$. Let $u\in H_0^1$. Define the piecewise linear interpolant by

$$\mathit{\pi u}=\sum _{j=1}^{N}u\left(x_j\right){\varphi }_j\left(x\right)\in V_{h,0}^1,x\in \left[a,b\right].$$

Since $\mathit{\pi u}\in V_{h,0}^1$, the estimate (14) gives

$$‖{\left(u-u_h\right)}^{\prime }‖\le C‖{\left(u-\mathit{\pi u}\right)}^{\prime }‖.$$

This inequality suggest that the error between $u$ and $u_h$ is controlled by the interpolation error $u-\mathit{\pi u}$ in the ${\left|\cdot \right|}_1$-norm.

Theorem 2.4 (A priori error estimate) Suppose that$q\left(x\right)\ge 0$$\forall x\in \left[a,b\right]$. Let$u$be the solution to (4) and$u_h$be the solution to (6). Then there exists a constant$C$such that

$${‖{\left(u-u_h\right)}^{\prime }‖}^2\le C\sum _{i=1}^N{h}_i^2{‖u^{″}‖}_{I_i}^2,$$

where $C$ is a constant independent of $h$. Consequently, if $h={max}_i{h}_i$, then

$${‖{\left(u-u_h\right)}^{\prime }‖}^2\le {\mathit{Ch}}^2{‖u^{″}‖}^2.$$

Remark 2.4

1. If the partition is not uniform then we obtain the same error estimate with $h={max}_{i=1,2,\dots ,N}\left(x_i-x_{i-1}\right)$.

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

3. $u_h\to u$ in the $‖v^{\prime }‖$-norm as $h={\mathit{max}}_i\left(h_i\right)\to 0$. If $‖{\left(u-u_h\right)}^{\prime }‖=0$ then $u-u_h$ is constant, but since $u\left(0\right)=u_h\left(0\right)$ we also have $u-u_h=0$ and therefore $u_h=u$.

4. $u_h$ is the best approximation within the space $V_{h,0}^1$ with respect to the $‖v^{\prime }‖$-norm.

5. The norm $‖v^{\prime }‖$ 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: $u\left(a\right)=\alpha$ and $u\left(b\right)=\beta$ for two real numbers $\alpha$ and $\beta$. This BC is also known as strong BC or essential BC.

2. Neumann BCs: $u^{\prime }\left(a\right)=\alpha$ and $u^{\prime }\left(b\right)=\beta$ for two real numbers $\alpha$ and $\beta$. This BC is also known as natural BCs.

3. Robin BCs: $u^{\prime }\left(a\right)=\mathit{\alpha u}\left(a\right)$ and $u^{\prime }\left(b\right)=\mathit{\beta u}\left(b\right)$ for two real numbers $\alpha$ and $\beta$.

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

Nonhomogeneous Dirichlet boundary conditions: Let us consider the following two-point BVP: find $u\in C^2\left(a,b\right)$ such that

$$-u^{″}=f\left(x\right),x\in \left(a,b\right),u\left(a\right)=\alpha ,u\left(b\right)=\beta ,$$E15

where $\alpha$ and $\beta$ are given constants and $f\in C\left(a,b\right)$ is a given function. In this case, the admissible function space $H_0^1=\left\{v:{‖v‖}^2+{‖v^{\prime }‖}^2<\infty ,v\left(a\right)=v\left(b\right)=0\right\}$ and the FE space $V_{h,0}^1$ defined earlier remain the same. Multiplying (15) by a test function $v\in H_0^1$ and integrating by parts gives

$${\int }_a^b\mathit{\text{fvdx}}={\int }_a^b-u^{″}\mathit{vdx}=-u^{\prime }\left(b\right)v\left(b\right)+u^{\prime }\left(a\right)v\left(a\right)+{\int }_a^b{u}^{\prime }{v}^{\prime }\mathit{dx}={\int }_a^b{u}^{\prime }{v}^{\prime }\mathit{dx},$$

since $v\left(a\right)=v\left(b\right)=0$. Hence, the weak or variational form of (15) reads: Given $u\left(a\right)=\alpha$, $u\left(b\right)=\beta$, find $u\in H^1=\left\{v:{‖v‖}^2+{‖v^{\prime }‖}^2<\infty \right\}$, such that

$${\int }_a^b{u}^{\prime }{v}^{\prime }\mathit{dx}={\int }_a^b\mathit{\text{fvdx}},\forall v\in H_0^1.$$E16

Let $V_h^1$ and $V_{h,0}^1$, 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=x_0<x_1<\cdots <x_N=b$ be a uniform partition of the interval $\left[a,b\right]$. Moreover let $\left\{{\varphi }_i\right\}$ be the set of hat basis functions of $V_h$ associated with the $N+1$ nodes $x_j,j=0,1,\dots ,N$, such that ${\varphi }_i\left(x_j\right)={\delta }_{\mathit{ij}}$. The FE approximation of (16) thus reads: Find $u_h\in V_h^1$ such that $u_h\left(a\right)=\alpha$, $u_h\left(b\right)=\beta$, and

$${\int }_a^b{u}_h^{\prime }{v}^{\prime }\mathit{dx}={\int }_a^b\mathit{\text{fvdx}},\forall v\in V_{h,0}^1.$$E17

It can be shown that (17) is equivalent to the $N-1$ equations

$${\int }_a^b{u}_h^{\prime }{\varphi }_i^{\prime }\mathit{dx}={\int }_a^{b}f{\varphi }_i\mathit{dx},i=1,2,\dots ,N-1.$$E18

Expanding $u_h$ as a linear combination of hat functions

$$u_h=\sum _{j=0}^N{c}_j{\varphi }_j\left(x\right)={\mathit{\alpha \varphi }}_0\left(x\right)+\sum _{j=1}^{N-1}{c}_j{\varphi }_j\left(x\right)+{\mathit{\beta \varphi }}_N\left(x\right),$$E19

where the coefficients $c_j,j=1,2,\dots ,N-1$ are the $N-1$ nodal values of $u_h$ to be determined.

Substituting (19) into (18) yields

$$\sum _{j=1}^{N-1}{c}_j\left({\int }_a^b{\varphi }_i^{\prime }{\varphi }_j^{\prime }\mathit{dx}\right)={\int }_a^b\left(f{\varphi }_i-{\mathit{\alpha \varphi }}_0^{\prime }{\varphi }_i^{\prime }-{\mathit{\beta \varphi }}_N^{\prime }{\varphi }_i^{\prime }\right)\mathit{dx},i=1,2,\dots ,N-1,$$

which is a $\left(N-1\right)\times \left(N-1\right)$ system of equations for $c_j$. In matrix form we write

$$A\mathbf{c}=\mathbf{b},$$E20

where $A$ is a $\left(N-1\right)\times \left(N-1\right)$ matrix, the so-called stiffness matrix, with entries

$$a_{i,j}={\int }_a^b{\varphi }_i^{\prime }{\varphi }_j^{\prime }\mathit{dx},i,j=1,2,\dots ,N-1,$$E21

$\mathbf{c}={\left[c_1,c_2,\dots ,c_{N-1}\right]}^t$ is a $\left(N-1\right)$ vector containing the unknown coefficients $c_j,j=1,2,\dots ,N-1$, and $\mathbf{b}$ is a $\left(N-1\right)$ vector, the so-called load vector, with entries

$$b_i={\int }_a^b\left(f{\varphi }_i-{\mathit{\alpha \varphi }}_0^{\prime }{\varphi }_i^{\prime }-{\mathit{\beta \varphi }}_N^{\prime }{\varphi }_i^{\prime }\right)\mathit{dx},i=1,2,\dots ,N-1.$$E22

Computer Implementation: The explicit expression for a hat function ${\varphi }_i\left(x\right)$ is given by

$$\begin{array}{l}{\varphi }_i\left(x\right)=\left\{\begin{array}{ll}0,& a\le x\le x_{i-1},\\ \frac{x-x_{i-1}}{h_i},& x_{i-1}<x\le x_i,\\ \frac{x_{i+1}-x}{h_{i+1}},& x_i<x\le x_{i+1},\\ 0,& x_{i+1}<x\le b,\end{array}\right.,i=1,2,\dots ,N-1,\\ {\varphi }_0\left(x\right)=\left\{\begin{array}{ll}\frac{x_1-x}{h_1},& x_0<x\le x_1,\\ 0,& x_1<x\le b,\end{array}\right.,{\varphi }_N\left(x\right)=\left\{\begin{array}{ll}0,& x_0<x\le x_{N-1},\\ \frac{x-x_{N-1}}{h_N},& x_{N-1}<x\le b.\end{array}\right.\end{array}$$

For simplicity we assume the partition is uniform so that $h_i=h$ for $i=1,2,\dots ,N$. Hence the derivative ${\varphi }_i^{\prime }\left(x\right)$ is either $-\frac{1}{h}$, $\frac{1}{h}$, or 0 depending on the interval.

It is straightforward to calculate the entries of the stiffness matrix. For $\mid i-j\mid >1$, we have $a_{i,j}=0$, since ${\varphi }_i$ and ${\varphi }_j$ lack overlapping support. However, if $i=j$, then

$$a_{i,j}={\int }_a^b{\left({\varphi }_i^{\prime }\right)}^2\mathit{dx}={\int }_{x_{i-1}}^{x_i}{\left(\frac{1}{h}\right)}^2\mathit{dx}+{\int }_{x_i}^{x_{i+1}}{\left(-\frac{1}{h}\right)}^2\mathit{dx}=\frac{2}{h},i,j=1,2,\dots ,N-1,$$

where we have used that $x_i-x_{i-1}=x_{i+1}-x_i=h$. Furthermore, if $j=i+1$, then

$$a_{i,i+1}={\int }_a^b{\varphi }_i^{\prime }{\varphi }_{i+1}^{\text{'}}\mathit{dx}={\int }_{x_i}^{x_{i+1}}\left(-\frac{1}{h}\right)\left(\frac{1}{h}\right)\mathit{dx}=-\frac{1}{h},i,j=1,2,\dots ,N-2.$$

Changing $i$ to $i-1$ we also have

$$a_{i-1,i}={\int }_a^b{\varphi }_{i-1}^{\text{'}}{\varphi }_i\mathit{dx}={\int }_{x_{i-1}}^{x_i}\left(\frac{1}{h}\right)\left(-\frac{1}{h}\right)\mathit{dx}=-\frac{1}{h},i,j=2,3,\dots ,N-1.$$

Thus the stiffness matrix is

$$A=\frac{1}{h}\left[\begin{array}{ccccc}2& -1& 0& \cdots & 0\\ -1& 2& -1& \ddots & 0\\ 0& -1& 2& \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & -1\\ 0& \cdots & 0& -1& 2\end{array}\right].$$

The entries $b_i$ of the load vector must often be evaluated using quadrature, since they involve the function $f$ which can be hard to integrate analytically. For example, using the trapezoidal rule one obtains the approximate load vector entries

$$\begin{array}{cc}{b}_1& ={\int }_a^b\left(f{\varphi }_1-{\mathit{\alpha \varphi }}_0^{\prime }{\varphi }_1^{\prime }-{\mathit{\beta \varphi }}_N^{\prime }{\varphi }_1^{\prime }\right)\mathit{dx}={\int }_{x_0}^{x_1}\left(f{\varphi }_1-\alpha \left(-\frac{1}{h}\right)\left(\frac{1}{h}\right)\right)\mathit{dx}+{\int }_{x_1}^{x_2}f{\varphi }_1\\ & =\frac{\alpha }{h}+{\int }_{x_0}^{x_2}f{\varphi }_1\approx \frac{\alpha }{h}+\mathit{hf}\left(x_1\right),\\ b_i& ={\int }_a^b\left(f{\varphi }_i-{\mathit{\alpha \varphi }}_0^{\prime }{\varphi }_i^{\prime }-{\mathit{\beta \varphi }}_N^{\prime }{\varphi }_i^{\prime }\right)\mathit{dx}={\int }_{x_{i-1}}^{x_{i+1}}f{\varphi }_i\mathit{dx}\approx \mathit{hf}\left(x_i\right),i=2,\dots ,N-2,\\ b_{N-1}& ={\int }_a^b\left(f{\varphi }_{N-1}^{\prime }-{\mathit{\alpha \varphi }}_0^{\prime }{\varphi }_{N-1}-{\mathit{\beta \varphi }}_N^{\prime }{\varphi }_{N-1}^{\prime }\right)\mathit{dx}\\ & ={\int }_{x_{N-2}}^{x_{N-1}}f{\varphi }_{N-1}\mathit{dx}+{\int }_{x_{N-1}}^{x_N}\left(f{\varphi }_{N-1}-\beta \left(\frac{1}{h}\right)\left(-\frac{1}{h}\right)\right)\mathit{dx}={\int }_{x_{N-2}}^{x_N}f{\varphi }_{N-1}\mathit{dx}+\frac{\beta }{h}\approx \mathit{hf}\left(x_{N-1}\right)+\frac{\beta }{h}.\end{array}$$

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

$$\frac{1}{h}\left[\begin{array}{ccccc}2& -1& 0& \cdots & 0\\ -1& 2& -1& \ddots & 0\\ 0& -1& 2& \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & -1\\ 0& \cdots & 0& -1& 2\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\\ ⋮\\ c_{N-1}\end{array}\right]=\left[\begin{array}{c}\mathit{hf}\left(x_1\right)+\frac{\alpha }{h}\\ \mathit{hf}\left(x_2\right)\\ \mathit{hf}\left(x_3\right)\\ ⋮\\ \mathit{hf}\left(x_{N-1}\right)+\frac{\beta }{h}\end{array}\right]$$

We note that $u_h\left(a\right)=\alpha =u\left(a\right)$ and $u_h\left(b\right)=\beta =u\left(b\right)$. Therefore, we see that the system matrix $A$ remains the same, and only the first and last entries of the load vector $\mathbf{b}$ need to be modified because of the definition of the basis functions $\left\{{\varphi }_0,\dots ,{\varphi }_N\right\}$. An alternative approach is to use all the basis functions $\left\{{\varphi }_0,\dots ,{\varphi }_N\right\}$ to form a larger system of equation, $i.e.$, and $\left(N+1\right)\times \left(N+1\right)$ system. The procedure for inserting the boundary conditions into the system equation is: enter zeros in the first and $\left(N+1\right)$-th rows of the system matrix $A$ except for unity in the main diagonal positions of these two rows, and enter $\alpha$ and $\beta$ in the first and $\left(N+1\right)$-th rows of the vector $\mathbf{b}$, respectively.

General boundary conditions: Let us consider the following two-point BVP: find $u\in C^2\left(a,b\right)$ such that

$$-u^{″}=f\left(x\right),x\in \left[a,b\right],u\left(a\right)=\alpha ,\mathit{\gamma u}\left(b\right)+u^{\prime }\left(b\right)=\beta ,$$E23

where $\alpha ,\beta$ and $\gamma$ are given numbers and $f\in C\left(a,b\right)$ is a given function. The boundary condition at $x=b$ is called a Robin boundary condition (combination and $u$ and $u^{\prime }$ is prescribed at $x=b$). In this case, the admissible function space is modified to

$$H_0^1=\left\{v:{‖v‖}^2+{‖v^{\prime }‖}^2<\infty ,v\left(a\right)=0\right\}.$$

Multiplying (23) by a function $v\in H_0^1$ and integrating by parts gives

$${\int }_a^b\mathit{\text{fvdx}}={\int }_a^b-u^{″}\mathit{vdx}=-u^{\prime }\left(b\right)v\left(b\right)+u^{\prime }\left(a\right)v\left(a\right)+{\int }_a^b{u}^{\prime }{v}^{\prime }\mathit{dx}$$

$$=-\left(\beta -\mathit{\gamma u}\left(b\right)\right)v\left(b\right)+u^{\prime }\left(a\right)v\left(a\right)+{\int }_a^b{u}^{\prime }{v}^{\prime }\mathit{dx}.$$

Since $v\left(a\right)=0$, we are left with

$${\int }_a^b{u}^{\prime }{v}^{\prime }\mathit{dx}+\mathit{\gamma u}\left(b\right)v\left(b\right)={\int }_a^b\mathit{\text{fvdx}}+\mathit{\beta v}\left(b\right).$$

Hence, the weak or variational form of (23) reads: Given $u\left(a\right)=\alpha$, find the approximate solution $u\in H_0^1$, such that

$${\int }_a^b{u}^{\prime }{v}^{\prime }\mathit{dx}+\mathit{\gamma u}\left(b\right)v\left(b\right)={\int }_a^b\mathit{\text{fvdx}}+\mathit{\beta v}\left(b\right),\forall v\in H_0^1.$$E24

The FE space $V_h^1$ is now the set of all continuous piecewise linear functions which vanish at the end point $a$. The FE approximation of (24) thus reads: Find the piecewise linear approximation $u_h$ to the solution $u$ satisfies

$${\int }_a^b{u}_h^{\prime }{v}^{\prime }\mathit{dx}+\gamma u_h\left(b\right)v\left(b\right)={\int }_a^b\mathit{\text{fvdx}}+\mathit{\beta v}\left(b\right),\forall v\in V_h^1,$$E25

with $u_h\left(a\right)=\alpha$. 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 $u\in C^2\left(a,b\right)$ such that

$$\begin{array}{l}-{\left(p\left(x\right)u^{\prime }\right)}^{\prime }=f\left(x\right),x\in I=\left[a,b\right],p\left(a\right)u^{\prime }\left(a\right)={\kappa }_0\left(u\left(a\right)-\alpha \right),\\ p\left(b\right)u^{\prime }\left(b\right)={\kappa }_1\left(u\left(b\right)-\beta \right),\end{array}$$E26

where $p=p\left(x\right)$ with $p\left(x\right)\ge p_0>0$, $f\in L^2\left(I\right)$, ${\kappa }_0,{\kappa }_1\ge 0$, and $\alpha ,\beta$ are given numbers. Let

$$V=\left\{v\in C^0\left(I\right):{‖v‖}^2+{‖v^{\prime }‖}^2<\infty \right\}.$$

Multiplying (26) by a function $v\in V$ and integrating by parts gives

$$\begin{array}{cc}{\int }_a^b\mathit{\text{fvdx}}& ={\int }_a^b-{\left(p{u}^{\prime }\right)}^{\prime }\mathit{vdx}={\int }_a^{b}p{u}^{\prime }{v}^{\prime }\mathit{dx}-p\left(b\right)u^{\prime }\left(b\right)v\left(b\right)+p\left(a\right)u^{\prime }\left(a\right)v\left(a\right)\\ & ={\int }_a^{b}p{u}^{\prime }{v}^{\prime }\mathit{dx}-{\kappa }_1\left(u\left(b\right)-\beta \right)v\left(b\right)+{\kappa }_0\left(u\left(a\right)-\alpha \right)v\left(a\right).\end{array}$$

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

$${\int }_a^{b}p{u}^{\prime }{v}^{\prime }\mathit{dx}-{\kappa }_{1}u\left(b\right)v\left(b\right)+{\kappa }_{0}u\left(a\right)v\left(a\right)={\int }_a^b\mathit{\text{fvdx}}-{\kappa }_1\mathit{\beta v}\left(b\right)+{\kappa }_0\mathit{\alpha v}\left(a\right),\forall v\in V.$$

The FE method consists of finding $u_h\in V_h=\left\{v\in C^0\left(a,b\right){\left|v\right|}_{I_i}\in P^1\left(I_i\right)\right\}$ such that

$${\int }_a^b{\mathit{pu}}_h^{\prime }{v}^{\prime }\mathit{dx}-{\kappa }_1{u}_h\left(b\right)v\left(b\right)+{\kappa }_0{u}_h\left(a\right)v\left(a\right)={\int }_a^b\mathit{\text{fvdx}}-{\kappa }_1\mathit{\beta v}\left(b\right)+{\kappa }_0\mathit{\alpha v}\left(a\right),\forall v\in V_h.$$E27

Implementation: We need to assemble a stiffness matrix $A$ and a load vector $b$. Substituting $u_h={\sum }_{i=0}^N{c}_i{\varphi }_i$ into (27) and taking $v={\varphi }_j$ for $j=0,1,\dots ,N$ yields

$$\begin{gathered} \sum _{i=0}^N{\int }_a^{b}p{\varphi }_i^{\prime }{\varphi }_j^{\prime }\mathit{dx}-{\kappa }_1{\varphi }_i\left(b\right){\varphi }_j\left(b\right)+{\kappa }_0{\varphi }_i\left(a\right){\varphi }_j\left(a\right)={\int }_a^{b}f{\varphi }_j\mathit{dx}-{\kappa }_1{\mathit{\beta \varphi }}_j\left(b\right)+{\kappa }_0{\mathit{\alpha \varphi }}_j\left(a\right), \\ \forall j=0,1,\dots ,N. \end{gathered}$$

which is a $\left(N+1\right)\times \left(N+1\right)$ system of equations for $c_i$. In matrix form we write $A\mathbf{c}=\mathbf{b}$, where $\mathbf{c}={\left[c_0,\dots ,c_N\right]}^t$ is a $\left(N+1\right)$ vector containing the unknown coefficients $c_i,i=0,1,\dots ,N$, $A$ is a $\left(N+1\right)\times \left(N+1\right)$ matrix with entries

$$a_{i,j}={\int }_a^{b}p{\varphi }_i^{\prime }{\varphi }_j^{\prime }\mathit{dx}-{\kappa }_1{\varphi }_i\left(b\right){\varphi }_j\left(b\right)+{\kappa }_0{\varphi }_i\left(a\right){\varphi }_j\left(a\right),i,j=0,1,\dots ,N,$$

and $\mathbf{b}$ is a $\left(N+1\right)$ vector with entries

$$b_j={\int }_a^{b}f{\varphi }_j\mathit{dx}-{\kappa }_1{\mathit{\beta \varphi }}_j\left(b\right)+{\kappa }_0{\mathit{\alpha \varphi }}_j\left(a\right),j=0,1,\dots ,N.$$

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

$$A=\left[\begin{array}{ccccc}{\kappa }_0+\frac{1}{h_1}& -\frac{1}{h_1}& 0& \cdots & 0\\ -\frac{1}{h_1}& \frac{1}{h_1}+\frac{1}{h_2}& -\frac{1}{h_2}& \ddots & 0\\ 0& -\frac{1}{h_2}& \ddots & \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & -\frac{1}{h_N}\\ 0& \cdots & 0& -\frac{1}{h_N}& \frac{1}{h_N}-{\kappa }_1\end{array}\right],\mathbf{b}=\left[\begin{array}{c}\frac{h_1}{2}{f}_0+{\kappa }_0\alpha \\ \frac{h_1+h_2}{2}{f}_1\\ ⋮\\ \frac{h_{N-1}+h_N}{2}{f}_{N-1}\\ \frac{h_N}{2}{f}_N-{\kappa }_1\beta \end{array}\right].$$

2.5 The FE method using Lagrange ${\mathbb{P}}_2$ elements

Let $a=x_0<x_1<\cdots <x_N=b$ be a regular partition of the interval $\left[a,b\right]$. Suppose that the length of $I_i=\left[x_{i-1},x_i\right]$ is $h_i=x_i-x_{i-1}$. Let $P^k=\left\{p\left(x\right)={\sum }_{j=0}^k{c}_j{x}^j,c_j\in \mathbb{R}\right\}$ denotes the vector space of polynomials in one variable and of degree less than or equal to $k$. The FE method for Lagrange $P^2$ elements involves the discrete space:

$${\left.V_h^2=\{v\left(x\right)\in C^0[a,b],v\right|}_{I_i}\in P^2\left(I_i\right),i=1,\dots ,N\},$$

and its subspace $V_{0,h}^2=\left\{v\in V_h^2|v\left(a\right)=v\left(b\right)=0\right\}.$ These spaces are composed of continuous, piecewise parabolic functions (polynomials of degree less than or equal to 2). The $P^2$ FE method consists in applying the internal variational approximation approach to these spaces.

Lemma 2.2The space$V_h^2$is a subspace of$H^1\left[a,b\right]$of dimension$2N+1$. Every function$v_h\in V_h^2$is uniquely defined by its values at the mesh vertices$x_j,j=0,1,\dots ,N$and at the midpoints$x_{j+\frac{1}{2}}=\frac{x_j+x_{j+1}}{2}=x_j+\frac{h_{j+1}}{2},j=0,1,\dots ,N-1$, where$h_{j+1}=x_{j+1}-x_j$:

$$v_h\left(x\right)=\sum _{j=0}^N{v}_h\left(x_j\right){\varphi }_j\left(x\right)+\sum _{j=0}^{N-1}{v}_h\left(x_{j+\frac{1}{2}}\right){\varphi }_{j+\frac{1}{2}}\left(x\right),\forall x\in \left[a,b\right],$$

where ${\left\{{\varphi }_j\right\}}_{j=0}^N$ is the basis of the shape functions ${\varphi }_j$ defined as:

$${\varphi }_j\left(x\right)=\varphi \left(\frac{x-x_j}{h_{j+1}}\right),j=0,1,\dots ,N,{\varphi }_{j+\frac{1}{2}}\left(x\right)=\psi \left(\frac{x-x_{j+\frac{1}{2}}}{h_{j+1}}\right),j=0,1,\dots ,N-1,$$

with

$$\varphi \left(\xi \right)=\left\{\begin{array}{ll}\left(1+\xi \right)\left(1+2\xi \right),& \xi \in \left[-1,0\right],\\ \left(1-\xi \right)\left(1-2\xi \right),& \xi \in \left[0,1\right],\\ 0,& \mid \xi \mid >1,\end{array}\right.\psi \left(\xi \right)=\left\{\begin{array}{ll}1-4{\xi }^2,& \mid \xi \mid \le \frac{1}{2},\\ 0,& \mid \xi \mid >\frac{1}{2},\end{array}\right.$$E28

Figure 3 shows the global shape functions for the space $V_h^2$ and the three quadratic Lagrange $P^2$ shape functions on the reference interval $\left[-1,1\right]$.

Figure 3.

(left) global shape functions for the space $V_h^2$. (right) the three quadratic Lagrange $P^2$ shape functions on the reference interval $\left[-1,1\right]$.

Remark 2.5Notice that we have:

$${\varphi }_j\left(x_j\right)={\delta }_{\mathit{ij}},{\varphi }_j\left(x_{j+\frac{1}{2}}\right)=0,{\varphi }_{j+\frac{1}{2}}\left(x_j\right)=0,{\varphi }_{j+\frac{1}{2}}\left(x_{j+\frac{1}{2}}\right)={\delta }_{\mathit{ij}}.$$

Corollary 2.1The space$V_{0,h}^2$is a subspace of$H_0^1\left[a,b\right]$of dimension$2N-1$and every function$v_h\in V_{0,h}^2$is uniquely defined by its values at the mesh vertices$x_j,j=1,2,\dots ,N-1$and at the midpoints$x_{j+\frac{1}{2}},j=0,1,\dots ,N-1$:

$$v_h\left(x\right)=\sum _{j=1}^{N-1}{v}_h\left(x_j\right){\varphi }_j\left(x\right)+\sum _{j=0}^{N-1}{v}_h\left(x_{j+\frac{1}{2}}\right){\varphi }_{j+\frac{1}{2}}\left(x\right),\forall x\in \left[a,b\right],$$

where ${\left\{{\varphi }_j\right\}}_{j=0}^N$ is the basis of the shape functions ${\varphi }_j$ defined as:

$${\varphi }_j\left(x\right)=\varphi \left(\frac{x-x_j}{h_{j+1}}\right),j=0,1,\dots ,N,{\varphi }_{j+\frac{1}{2}}\left(x\right)=\psi \left(\frac{x-x_{j+\frac{1}{2}}}{h_{j+1}}\right),j=0,1,\dots ,N-1,$$

with $\varphi \left(\xi \right)$ and $\psi \left(\xi \right)$ 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 $u_h\in V_{0,h}^2$, such that:

$${\int }_a^b{u}_h^{\prime }{v}^{\prime }\mathit{dx}+{\int }_a^b{\mathit{qu}}_h\mathit{vdx}={\int }_a^b\mathit{\text{fvdx}},\forall v\in V_{h,0}^2.$$

Here, it is convenient to introduce the notation $x_{\frac{j}{2}},j=1,\dots ,2N-1$ for the mesh points and ${\varphi }_{\frac{j}{2}},j=1,\dots ,2N-1$ for the basis of $V_{0,h}^2$. Using these notations, we have:

$$u_h=\sum _{j=1}^{2N-1}{c}_{\frac{j}{2}}{\varphi }_{\frac{j}{2}}\left(x\right),$$

where $c_{\frac{j}{2}}=u_h\left(x_{\frac{j}{2}}\right)\approx u\left(x_{\frac{j}{2}}\right)$ are the unknowns coefficients. This formulation leads to solve in ${\mathbb{R}}^{2N-1}$ a linear system:

$$A\mathbf{c}=\mathbf{b},$$

where $\mathbf{c}={\left[c_{\frac{1}{2}},c_1,\dots ,c_{N-\frac{1}{2}}\right]}^t\in {\mathbb{R}}^{2N-1}$ is the unknown vector containing the coefficients $c_{\frac{j}{2}},j=1,2,\dots ,2N-1$, $A$ is an $\left(2N-1\right)\times \left(2N-1\right)$ matrix with entries

$$a_{\mathit{ij}}={\int }_a^b\left({\varphi }_{\frac{i}{2}}\text{'}{\varphi }_{\frac{j}{2}}\text{'}+q{\varphi }_{\frac{i}{2}}{\varphi }_{\frac{j}{2}}\right)\mathit{dx},i,j=1,2,\dots ,2N-1,$$

and load vector $\mathbf{b}\in {\mathbb{R}}^{2N-1}$ has entries

$$b_{\frac{i}{2}}={\int }_a^{b}f{\varphi }_{\frac{i}{2}}\mathit{dx},i=1,2,\dots ,2N-1.$$

Since the shape functions ${\varphi }_i$ have a small support, the matrix $A$ is mostly composed of zeros. However, the main difference with the Lagrange $P^1$ FE method, the matrix $A$ is no longer a tridiagonal matrix.

Computer Implementation: The coefficients of the matrix $A$ can be computed more easily by considering the following change of variables, for $\xi \in \left[-1,1\right]$:

$$x=\frac{x_j+x_{j-1}}{2}+\frac{x_j-x_{j-1}}{2}\xi =x_{j-\frac{1}{2}}+\frac{x_j-x_{j-1}}{2}\xi ,\forall x\in \left[x_{j-1},x_j\right],j=1,2,\dots ,N.$$

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

$${\widehat{\varphi }}_{-1}\left(\xi \right)=\frac{\xi \left(\xi -1\right)}{2},{\widehat{\varphi }}_0\left(\xi \right)=\left(1-\xi \right)\left(1+\xi \right),{\widehat{\varphi }}_1\left(\xi \right)=\frac{\xi \left(\xi +1\right)}{2}.$$

Their respective derivatives are

$$\frac{d{\widehat{\varphi }}_{-1}\left(\xi \right)}{\mathit{d\xi }}=\frac{2\xi -1}{2},\frac{d{\widehat{\varphi }}_0\left(\xi \right)}{\mathit{d\xi }}=-2\xi ,\frac{d{\widehat{\varphi }}_1\left(\xi \right)}{\mathit{d\xi }}=\frac{2\xi +1}{2}.$$

This approach consists in considering all computations on an interval $I_i=\left[x_{i-1},x_i\right]$ on the reference interval $\left[-1,1\right]$. Thus, we have:

$$\frac{d{\varphi }_i\left(x\right)}{\mathit{dx}}=\frac{d{\varphi }_i\left(x_{i-1/2}+\frac{x_i-x_{i-1}}{2}\xi \right)}{\mathit{d\xi }}\frac{\mathit{d\xi }}{\mathit{dx}}=\frac{2}{x_i-x_{i-1}}\frac{d{\widehat{\varphi }}_k\left(\xi \right)}{\mathit{d\xi }}=\frac{2}{h_i}\frac{d{\widehat{\varphi }}_k\left(\xi \right)}{\mathit{d\xi }}.$$

In this case, the elementary contributions of the element $I_i$ to the stiffness matrix and to the mass matrix are given by the $3\times 3$ matrices $K^{I_i}$ and $M^{I_i}$:

$$K^{I_i}={\int }_{I_i}\left[\begin{array}{ccc}{\varphi }_{i-1}^{\prime }{\varphi }_{i-1}^{\prime }& {\varphi }_{i-1}^{\prime }{\varphi }_{i-\frac{1}{2}}^{\prime }& {\varphi }_{i-1}^{\prime }{\varphi }_i^{\prime }\\ {\varphi }_{i-\frac{1}{2}}^{\prime }{\varphi }_{i-1}^{\prime }& {\varphi }_{i-\frac{1}{2}}^{\prime }{\varphi }_{i-\frac{1}{2}}^{\prime }& {\varphi }_{i-\frac{1}{2}}^{\prime }{\varphi }_i^{\prime }\\ {\varphi }_i^{\prime }{\varphi }_{i-1}^{\prime }& {\varphi }_i^{\prime }{\varphi }_{i-\frac{1}{2}}^{\prime }& {\varphi }_i^{\prime }{\varphi }_i^{\prime }\end{array}\right]\mathit{dx}=\frac{2}{h_i}{\int }_{-1}^1\left[\begin{array}{ccc}{\widehat{\varphi }}_{-1}^{\prime }{\widehat{\varphi }}_{-1}^{\prime }& {\widehat{\varphi }}_{-1}^{\prime }{\widehat{\varphi }}_0^{\prime }& {\widehat{\varphi }}_{-1}^{\prime }{\widehat{\varphi }}_1^{\prime }\\ {\widehat{\varphi }}_0^{\prime }{\widehat{\varphi }}_{-1}^{\prime }& {\widehat{\varphi }}_0^{\prime }{\widehat{\varphi }}_0^{\prime }& {\widehat{\varphi }}_0^{\prime }{\widehat{\varphi }}_1^{\prime }\\ {\widehat{\varphi }}_1^{\prime }{\widehat{\varphi }}_{-1}^{\prime }& {\widehat{\varphi }}_1^{\prime }{\widehat{\varphi }}_0^{\prime }& {\widehat{\varphi }}_1^{\prime }{\widehat{\varphi }}_1^{\prime }\end{array}\right]\mathit{d\xi }=\frac{1}{3h_i}\left[\begin{array}{ccc}7& -8& 1\\ -8& 16& -8\\ 1& -8& 7\end{array}\right],$$

$$M^{I_i}={\int }_{I_i}\left[\begin{array}{ccc}{\varphi }_{i-1}{\varphi }_{i-1}& {\varphi }_{i-1}{\varphi }_{i-\frac{1}{2}}& {\varphi }_{i-1}{\varphi }_i\\ {\varphi }_{i-\frac{1}{2}}{\varphi }_{i-1}& {\varphi }_{i-\frac{1}{2}}{\varphi }_{i-\frac{1}{2}}& {\varphi }_{i-\frac{1}{2}}{\varphi }_i\\ {\varphi }_i{\varphi }_{i-1}& {\varphi }_i{\varphi }_{i-\frac{1}{2}}& {\varphi }_i{\varphi }_i\\ & & \end{array}\right]\mathit{dx}=\frac{h_i}{2}{\int }_{-1}^1\left[\begin{array}{ccc}{\widehat{\varphi }}_{-1}{\widehat{\varphi }}_{-1}& {\widehat{\varphi }}_{-1}{\widehat{\varphi }}_0& {\widehat{\varphi }}_{-1}{\widehat{\varphi }}_1\\ {\widehat{\varphi }}_0{\widehat{\varphi }}_{-1}& {\widehat{\varphi }}_0{\widehat{\varphi }}_0& {\widehat{\varphi }}_0{\widehat{\varphi }}_1\\ {\widehat{\varphi }}_1{\widehat{\varphi }}_{-1}& {\widehat{\varphi }}_{-1}{\widehat{\varphi }}_0& {\widehat{\varphi }}_1{\widehat{\varphi }}_1\\ & & \end{array}\right]\mathit{d\xi }=\frac{h_i}{30}\left[\begin{array}{ccc}4& 2& -1\\ 2& 16& 2\\ -1& 2& 4\end{array}\right].$$

Coefficients of the right-hand side$\mathbf{b}$: Usually, the function $f$ is only known by its values at the mesh points $x_{\frac{i}{2}},i=0,1,\dots ,2N$ and thus, we use the decomposition of $f$ in the basis of shape functions ${\varphi }_{\frac{i}{2}},i=0,1,\dots ,2N$ as $f\left(x\right)={\sum }_{j=0}^{2N}f\left(x_{\frac{j}{2}}\right){\varphi }_{\frac{j}{2}}$. Each component $b_{\frac{i}{2}}$ of the right-hand side vector is obtained as $b_{\frac{i}{2}}={\sum }_{k=1}^N{\int }_{x_{k-1}}^{x_k}f{\varphi }_{\frac{i}{2}}\mathit{dx}$. Using the previous decomposition of $f$, we obtain:

$$b_{\frac{i}{2}}=\sum _{k=1}^N{\int }_{x_{k-1}}^{x_k}\sum _{j=0}^{2N}f\left(x_{\frac{j}{2}}\right){\varphi }_{\frac{j}{2}}{\varphi }_{\frac{i}{2}}\mathit{dx}=\sum _{j=0}^{2N}f\left(x_{\frac{j}{2}}\right)\left(\sum _{k=1}^N{\int }_{x_{k-1}}^{x_k}{\varphi }_{\frac{i}{2}}{\varphi }_{\frac{j}{2}}\mathit{dx}\right).$$

Thus, the problem is reduced to computing the integrals ${\int }_{x_{k-1}}^{x_k}{\varphi }_{\frac{i}{2}}{\varphi }_{\frac{j}{2}}\mathit{dx}$. It is easy to see that we obtain expressions very similar to that of the mass matrix. More precisely, the element $I_i=\left[x_{i-1},x_i\right]$ will contribute to only three components of indices $i-1$, $i-\frac{1}{2}$ and $i$ as:

$${\mathbf{b}}^{I_i}=\frac{h_i}{30}\left[\begin{array}{ccc}4& 2& -1\\ 2& 16& 2\\ -1& 2& 4\end{array}\right]\left[\begin{array}{c}f\left(x_{i-1}\right)\\ f\left(x_{i-\frac{1}{2}}\right)\\ f\left(x_i\right).\end{array}\right].$$

2.5.2 Nonhomogeneous boundary conditions

Consider the following two-point BVP: find $u\in C^2\left(a,b\right)$ such that

$$-u^{″}+q\left(x\right)u=f\left(x\right),x\in \left[a,b\right],u\left(a\right)=\alpha ,u\left(b\right)=\beta ,$$E29

where $\alpha$ and $\beta$ are given constants and $f\in C\left(a,b\right)$ is a given function.

Multiplying (29) by a function $v\in H_0^1=\left\{v:{‖v‖}^2+{‖v^{\prime }‖}^2<\infty ,v\left(a\right)=v\left(b\right)=0\right\}$ and integrating by parts gives

$${\int }_a^b\mathit{\text{fvdx}}={\int }_a^b\left(-u^{″}+\mathit{qu}\right)\mathit{vdx}=-u^{\prime }\left(b\right)v\left(b\right)+u^{\prime }\left(a\right)v\left(a\right)+{\int }_a^b\left(u^{\prime }{v}^{\prime }+\mathit{quv}\right)\mathit{dx}={\int }_a^b{u}^{\prime }{v}^{\prime }\mathit{dx}.$$

Hence, the weak or variational form of (29) reads: Given $u\left(a\right)=\alpha$, $u\left(b\right)=\beta$, find $u\in H^1=\left\{v:{‖v‖}^2+{‖v^{\prime }‖}^2<\infty \right\}$, such that

$${\int }_a^b\left(u^{\prime }{v}^{\prime }+\mathit{quv}\right)\mathit{dx}={\int }_a^b\mathit{\text{fvdx}},\forall v\in H_0^1.$$

Let $V_h^2$ and $V_{h,0}^2$, 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=x_0<x_1<\cdots <x_N=b$ of the interval $\left[a,b\right]$.

The FE method scheme consists of finding $u_h\in V_h^2$, such that:

$${\int }_a^b{u}_h^{\prime }{v}^{\prime }\mathit{dx}+{\int }_a^b{\mathit{qu}}_h\mathit{vdx}={\int }_a^b\mathit{\text{fvdx}},\forall v\in V_{h,0}^2.$$

Introduce the notation $x_{\frac{j}{2}},j=0,1,\dots ,2N-1,2N$ for the mesh points and ${\varphi }_{\frac{j}{2}},j=0,1,\dots ,2N-1,2N$ for the basis of $V_h^2$ and ${\varphi }_{\frac{j}{2}},j=1,\dots ,2N-1$ for the basis of $V_{0,h}^2$. Using these notations, we have:

$$u_h=\sum _{j=0}^{2N}{c}_{\frac{j}{2}}{\varphi }_{\frac{j}{2}}\left(x\right),$$

where $c_{\frac{j}{2}}=u_h\left(x_{\frac{j}{2}}\right)\approx u\left(x_{\frac{j}{2}}\right)$ are the unknowns coefficients. We note that $c_0=u_h\left(x_0\right)=\alpha$ and $c_{2N}=u_h\left(x_N\right)=\beta$. This formulation leads to solve in ${\mathbb{R}}^{2N-1}$ a linear system:

$$A\mathbf{c}=\mathbf{b},$$

where $\mathbf{c}={\left[c_{\frac{1}{2}},c_1,\dots ,c_{N-\frac{1}{2}}\right]}^t\in {\mathbb{R}}^{2N-1}$ is the unknown vector containing the coefficients $c_{\frac{j}{2}},j=1,2,\dots ,2N-1$, $A$ is an $\left(2N-1\right)\times \left(2N-1\right)$ matrix with entries

$$a_{\mathit{ij}}={\int }_a^b\left({\varphi }_{\frac{i}{2}}\prime {\varphi }_{\frac{j}{2}}\prime +q{\varphi }_{\frac{i}{2}}{\varphi }_{\frac{j}{2}}\right)\mathit{dx},i,j=1,2,\dots ,2N-1,$$

and the load vector $\mathbf{b}\in {\mathbb{R}}^{2N-1}$ has entries

$$b_{\frac{i}{2}}={\int }_a^{b}f{\varphi }_{\frac{i}{2}}\mathit{dx}-\alpha {\int }_a^b\left({\varphi }_{\frac{i}{2}}\prime {\varphi }_0^{\prime }+q{\varphi }_{\frac{i}{2}}{\varphi }_0\right)\mathit{dx}-\beta {\int }_a^b\left({\varphi }_{\frac{i}{2}}\prime {\varphi }_N\prime +q{\varphi }_{\frac{i}{2}}{\varphi }_N\right)\mathit{dx},i=1,2,\dots ,2N-1.$$

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

$${\tilde{b}}_{\frac{i}{2}}=-\alpha {\int }_a^b\left({\varphi }_{\frac{i}{2}}\prime {\varphi }_0^{\prime }+q{\varphi }_{\frac{i}{2}}{\varphi }_0\right)\mathit{dx}-\beta {\int }_a^b\left({\varphi }_{\frac{i}{2}}\prime {\varphi }_N\prime +q{\varphi }_{\frac{i}{2}}{\varphi }_N\right)\mathit{dx},i=1,2,\dots ,2N-1.$$

Suppose $q=0$ then for $N\ge 2$, we have

$$\begin{array}{ll}{\tilde{b}}_{\frac{1}{2}}& =-\alpha {\int }_a^b{\varphi }_{\frac{1}{2}}^{\prime }{\varphi }_0^{\prime }\mathit{dx}-\beta {\int }_a^b{\varphi }_{\frac{1}{2}}^{\prime }{\varphi }_N^{\prime }\mathit{dx}=-\alpha {\int }_{x_0}^{x_1}{\varphi }_{\frac{1}{2}}^{\prime }{\varphi }_0^{\prime }=\frac{8\alpha }{3h_1},\\ {\tilde{b}}_1& =-\alpha {\int }_a^b{\varphi }_1^{\prime }{\varphi }_0^{\prime }\mathit{dx}-\beta {\int }_a^b{\varphi }_1^{\prime }{\varphi }_N^{\prime }\mathit{dx}=-\alpha {\int }_{x_0}^{x_1}{\varphi }_1^{\prime }{\varphi }_0^{\prime }=-\frac{\alpha }{3h_1},\\ {\tilde{b}}_{\frac{i}{2}}& =-\alpha {\int }_a^b{\varphi }_{\frac{i}{2}}^{\prime }{\varphi }_0^{\prime }\mathit{dx}-\beta {\int }_a^b{\varphi }_{\frac{i}{2}}^{\prime }{\varphi }_N^{\prime }\mathit{dx}=0,i=3,\dots ,2N-3,\\ {\tilde{b}}_{N-1}& =-\alpha {\int }_a^b{\varphi }_{N-1}^{\prime }{\varphi }_0^{\prime }\mathit{dx}-\beta {\int }_a^b{\varphi }_{N-1}^{\prime }{\varphi }_N^{\prime }\mathit{dx}=-\beta {\int }_{x_{N-1}}^{x_N}{\varphi }_{N-1}^{\prime }{\varphi }_N^{\prime }\mathit{dx}=-\frac{\beta }{3h_1},\\ {\tilde{b}}_{N-\frac{1}{2}}& =-\alpha {\int }_a^b{\varphi }_{N-\frac{1}{2}}^{\prime }{\varphi }_0^{\prime }\mathit{dx}-\beta {\int }_a^b{\varphi }_{N-\frac{1}{2}}^{\prime }{\varphi }_N^{\prime }\mathit{dx}=-\beta {\int }_{x_{N-1}}^{x_N}{\varphi }_{N-\frac{1}{2}}^{\prime }{\varphi }_N^{\prime }\mathit{dx}=\frac{8\beta }{3h_1}.\end{array}$$

3. The FE for elliptic PDEs

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

$$-\nabla \cdot \left(a\nabla u\right)+\mathit{bu}=f\left(\mathbf{x}\right),\mathbf{x}\in \mathrm{\Omega },a\nabla u\cdot \mathbf{n}=\kappa \left(g-u\right),\mathrm{on}\mathrm{\partial \Omega },$$E30

where $a>0$, $b\ge 0$, $\kappa \ge 0$, $f\in L^2\left(\mathrm{\Omega }\right)$ and $g\in C^0\left(\mathrm{\partial \Omega }\right)$.

3.1 Meshes

Let $\mathrm{\Omega }\subset {\mathbb{R}}^2$ bounded with $\mathrm{\partial \Omega }$ assumed to be polygonal. A triangulation ${\mathcal{T}}_h$ of $\mathrm{\Omega }$ is a set of triangles $T$ such that $\mathrm{\Omega }={\bigcup }_{T\in {\mathcal{T}}_h}T$, and two triangles intersect by either a common triangle edge, or a corner, or nothing. Corners will be referred to as nodes. We let $h_T=\mathit{\text{diam}}\left(T\right)$ the length or the largest edge.

Let ${\mathcal{T}}_h$ have $N$ nodes and $M$ triangles. The data is stored in two matrices. The matrix $P\in {\mathbb{R}}^{2\times N}$ describes the nodes ($\left(x_1,y_1\right),\dots ,\left(x_N,y_N\right))$ and the matrix $K\in {\mathbb{R}}^{3\times 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=\left[\begin{array}{cccc}{x}_1& x_2& \cdots & x_N\\ y_1& y_2& \cdots & y_N\end{array}\right],K=\left[\begin{array}{cccc}{n}_1^{\alpha }& n_2^{\alpha }& \cdots & n_M^{\alpha }\\ n_1^{\beta }& n_2^{\beta }& \cdots & n_M^{\beta }\\ n_1^{\gamma }& n_2^{\gamma }& \cdots & n_M^{\gamma }\end{array}\right].$$

This means that triangle $T_i$ is formed by the nodes $n_i^{\alpha }$, $n_i^{\beta }$, and $n_i^{\gamma }$ (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 $N_1=\left(x_1,y_1\right)$, $N_2=\left(x_2,y_2\right)$, and $N_3=\left(x_3,y_3\right)$. We define

$$P^1\left(T\right)=\left\{v\in C^0\left(T\right)|v\left(x,y\right)=c_1+c_{2}x+c_{3}y,c_1,c_2,c_3\in \mathbb{R}\right\}.$$

Now let $v_i=v\left(N_i\right)$ for $i=\mathrm{1,2,3}$. Note that $v\in P^1\left(T\right)$ is determined by ${\left\{v_i\right\}}_{i=1}^3$. Given $v_i$ we compute $c_i$ by

$$\left[\begin{array}{ccc}1& x_1& y_1\\ 1& x_2& y_2\\ 1& x_3& y_3\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right]=\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right].$$

This is solvable due to

$$det\left[\begin{array}{ccc}1& x_1& y_1\\ 1& x_2& y_2\\ 1& x_3& y_3\end{array}\right]=2\mid T\mid \ne 0,{\left[\begin{array}{ccc}1& x_1& y_1\\ 1& x_2& y_2\\ 1& x_3& y_3\end{array}\right]}^{-1}=\frac{1}{2\mid T\mid }\left[\begin{array}{ccc}{x}_2{y}_3-x_3{y}_2& x_3{y}_1-x_1{y}_3& x_1{y}_2-x_2{y}_1\\ y_2-y_3& y_3-y_1& y_1-y_2\\ x_3-x_2& x_1-x_3& x_3-x_1\end{array}\right],$$

where $\mid T\mid =\frac{1}{2}\left(x_2{y}_3-x_3{y}_2-x_1{y}_3+x_3{y}_1+x_1{y}_2-x_2{y}_1\right)$, which is $\pm$ the area of the triangle $T$.

Let ${\lambda }_j\in P^1\left(T\right)$ be given by the nodal values ${\lambda }_j\left(N_i\right)={\delta }_{\mathit{ij}}$, where ${\delta }_{\mathit{ij}}$ is the Kronecker symbol. This gives us $v\left(x,y\right)={\alpha }_1{\lambda }_1\left(x,y\right)+{\alpha }_2{\lambda }_2\left(x,y\right)+{\alpha }_3{\lambda }_3\left(x,y\right)$, where ${\alpha }_i=v\left(N_i\right)$ for $i=\mathrm{1,2,3}.$ We can compute ${\lambda }_i\left(x,y\right)$ as follows: Let ${\lambda }_i\left(x,y\right)=a_i+b_{i}x+c_{i}y$. Using ${\lambda }_j\left(N_i\right)={\delta }_{\mathit{ij}}$, we get

$$\left[\begin{array}{ccc}1& x_1& y_1\\ 1& x_2& y_2\\ 1& x_3& y_3\end{array}\right]\left[\begin{array}{c}{a}_1\\ b_1\\ c_1\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right],\left[\begin{array}{ccc}1& x_1& y_1\\ 1& x_2& y_2\\ 1& x_3& y_3\end{array}\right]\left[\begin{array}{c}{a}_2\\ b_2\\ c_2\end{array}\right]=\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right],\left[\begin{array}{ccc}1& x_1& y_1\\ 1& x_2& y_2\\ 1& x_3& y_3\end{array}\right]\left[\begin{array}{c}{a}_3\\ b_3\\ c_3\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right].$$

Solving the systems, we get

$${\lambda }_1\left(x,y\right)=\frac{1}{2\mid T\mid }\left(x_2{y}_3-x_3{y}_2+\left(y_2-y_3\right)x+\left(x_3-x_2\right)y\right),$$

$${\lambda }_2\left(x,y\right)=\frac{1}{2\mid T\mid }\left(x_3{y}_1-x_1{y}_3+\left(y_3-y_1\right)x+\left(x_1-x_3\right)y\right),$$

$${\lambda }_3\left(x,y\right)=\frac{1}{2\mid T\mid }\left(x_1{y}_2-x_2{y}_1+\left(y_1-y_2\right)x+\left(x_3-x_1\right)y\right).$$

Let ${\mathcal{T}}_h$ be a triangulation of $\mathrm{\Omega }$, then we let

$$V_h=\left\{v\in C\left(\mathrm{\Omega }\right){\left|v\right|}_T\in P^1\left(T\right),\forall T\in {\mathcal{T}}_h\right\}.$$

Functions in $V_h$ are piecewise linear and continuous. We know that $v\in V_h$ is uniquely determined by $\left\{v\left(N_i\right),i=1,2,\dots ,N\right\}$. We let ${\varphi }_j\left(N_i\right)={\delta }_{\mathit{ij}}$ and let $\left\{{\varphi }_j,j=1,2,\dots ,N\right\}\subset V_h$ be a basis for $V_h$ (hat functions), $i.e.,$

$$v\left(x,y\right)=\sum _{i=1}^N{\alpha }_i{\varphi }_i\left(x,y\right),{\alpha }_i=v\left(N_i\right),i=1,2,\dots ,N.$$

3.3 Interpolation

Given $u\in C\left(T\right)$ on a single triangle with nodes $N_i=\left(x_i,y_i\right)$, $i=\mathrm{1,2,3}$, we let

$$\mathit{\pi u}\left(x,y\right)=\sum _{i=1}^{3}u\left(N_i\right){\varphi }_i\left(x,y\right),$$

in particular $\mathit{\pi u}\left(N_i\right)=u\left(N_i\right)$, $i=1,2,\dots ,N$. We want to estimate the interpolation error $u-\mathit{\pi u}$. Let

$$\begin{array}{cc}{‖u‖}_{L^2\left(\mathrm{\Omega }\right)}^2& ={\int }_{\mathrm{\Omega }}{\left|u\left(x\right)\right|}^2\mathit{\text{dxdy}},{‖\mathit{Du}‖}_{L^2\left(\mathrm{\Omega }\right)}^2={‖u_x‖}_{L^2\left(\mathrm{\Omega }\right)}^2+{‖u_y‖}_{L^2\left(\mathrm{\Omega }\right)}^2,\\ {‖D^{2}u‖}_{L^2\left(\mathrm{\Omega }\right)}^2& ={‖u_{\mathit{xx}}‖}_{L^2\left(\mathrm{\Omega }\right)}^2+2{‖u_{\mathit{xy}}‖}_{L^2\left(\mathrm{\Omega }\right)}^2+{‖u_{\mathit{yy}}‖}_{L^2\left(\mathrm{\Omega }\right)}^2.\end{array}$$

Theorem 3.1Suppose that$u\in C^2\left(T\right)$. Then the following hold

$${‖u-\mathit{\pi u}‖}_{L^2\left(T\right)}\le {\mathit{Ch}}_T^2{‖D^{2}u‖}_{L^2\left(T\right)},{‖D\left(u-\mathit{\pi u}\right)‖}_{L^2\left(T\right)}\le {\mathit{Ch}}_T{‖D^{2}u‖}_{L^2\left(T\right)},$$

where $C$ is a generic constant independent of $h_T$ and $u$, but it depends on the ratio between smallest and largest interior angle of the triangle $T$.

Now, we consider the piecewise continuous interpolant $\mathit{\pi u}={\sum }_{i=1}^{N}u\left(N_i\right){\varphi }_i$.

Theorem 3.2Suppose that$u\in C^2\left(T\right)$for all$T\in {\mathcal{T}}_h$. Then the following hold

$${‖u-\mathit{\pi u}‖}_{L^2\left(\mathrm{\Omega }\right)}^2\le C\sum _{T\in {\mathcal{T}}_h}{h}_T^4{‖D^{2}u‖}_{L^2\left(T\right)}^2,{‖D\left(u-\mathit{\pi u}\right)‖}_{L^2\left(\mathrm{\Omega }\right)}^2\le C\sum _{T\in {\mathcal{T}}_h}{h}_T^2{‖D^{2}u‖}_{L^2\left(T\right)}^2,$$

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 ${\mathcal{T}}_h$. Here ${‖D\left(u-\mathit{\pi u}\right)‖}_{L^2\left(\mathrm{\Omega }\right)}^2={\sum }_{T\in {\mathcal{T}}_h}{‖D\left(u-\mathit{\pi u}\right)‖}_{L^2\left(T\right)}^2$.

3.4 $L^2$-projection

Let $\mathrm{\Omega }\subset {\mathbb{R}}^2$. We consider the space $L^2\left(\mathrm{\Omega }\right)=\left\{v|{\int }_{\mathrm{\Omega }}{v}^2\left(x,y\right)\mathit{\text{dxdy}}<\infty \right\}$. Let $u\in L^2\left(\mathrm{\Omega }\right).$ We define the $L^2$-projection $P_h:L^2\left(\mathrm{\Omega }\right)\to V_h=\left\{v\in C^0\left(\mathrm{\Omega }\right){\left|v\right|}_T\in P^1\left(T\right),\forall T\in {\mathcal{T}}_h\right\}$ by $P_{h}u\in V_h$ such that

$${\int }_{\mathrm{\Omega }}\left(u-P_{h}u\right)v_h\mathit{\text{dxdy}}=0,\forall v_h\in V_h.$$

The problem of finding $P_{h}u\in V_h$ is equivalent to solve the following linear system

$${\int }_{\mathrm{\Omega }}\left(u-P_{h}u\right){\varphi }_i\mathit{\text{dxdy}}=0,i=1,2\dots ,N,$$

where ${\left\{{\varphi }_i\right\}}_{i=1}^N$ is a basis of $V_h$.

Since $P_{h}u\in V_h$ we can express it as $P_{h}u={\sum }_{i=1}^N{c}_i{\varphi }_i\left(x,y\right)$, where $c_i\in \mathbb{R}$. Therefore, to find $P_{h}u\in V_h$ we need to find $c_1,c_2,\dots ,c_N\in \mathbb{R}$ such that

$$\sum _{i=1}^N{c}_i{\int }_{\mathrm{\Omega }}{\varphi }_i{\varphi }_j\mathit{\text{dxdy}}={\int }_{\mathrm{\Omega }}u{\varphi }_j\mathit{\text{dxdy}},j=1,2,\dots ,N.$$

The problem can be expressed as a linear system of equations $M\mathbf{c}=\mathbf{b}$, where $\mathbf{c}={\left[c_1,c_2,\dots ,c_N\right]}^t$ and the entries of the matrix $M\in {\mathbb{R}}^{N\times N}$ and the vector $\mathbf{b}\in {\mathbb{R}}^N$ are given by

$$m_{\mathit{ij}}={\int }_{\mathrm{\Omega }}{\varphi }_i{\varphi }_j\mathit{\text{dxdy}},b_j={\int }_{\mathrm{\Omega }}u{\varphi }_j\mathit{\text{dxdy}}.$$

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

$${\int }_{T}f\left(x,y\right)\mathit{\text{dxdy}}\approx \sum _{j=1}^n{\omega }_{j}f\left({\overline{N}}_j\right),$$

where the ${\omega }_j$ ′s denote the weights and the $\left({\overline{N}}_j\right)$ ′s the quadrature points.

Lemma 3.1The mass matrix$M$with entries$m_{\mathit{ij}}={\int }_{\mathrm{\Omega }}{\varphi }_i{\varphi }_j\mathit{\text{dxdy}}$is symmetric and positive definite.

Theorem 3.3For any$u\in L^2\left(\mathrm{\Omega }\right)$the$L^2$-projection$P_{h}u$exists and is unique.

3.5 A priori error estimate

Theorem 3.4Let$u\in L^2\left(\mathrm{\Omega }\right)$and let$P_{h}u$be the$L^2$-projection of$u$, then

$${‖u-P_{h}u‖}_{L^2\left(\mathrm{\Omega }\right)}\le {‖u-v_h‖}_{L^2\left(\mathrm{\Omega }\right)},\forall v_h\in V_h.$$

Theorem 3.5Suppose that$u\in C^2\left(\mathrm{\Omega }\right)$with$u\in C^2\left(T\right)$for all$T\in {\mathcal{T}}_h$. Then there exists a constant$C$such that

$${‖u-P_{h}u‖}_{L^2\left(\mathrm{\Omega }\right)}^2\le C\sum _{T\in {\mathcal{T}}_h}{h}_T^4{‖D^{2}u‖}_{L^2\left(T\right)}^2.$$

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

$$-\nabla \cdot \left(a\nabla u\right)+\mathit{bu}=f,\text{in}\mathrm{\Omega },a\nabla u\cdot \mathbf{n}=\kappa \left(g-u\right),\mathrm{on}\mathrm{\partial \Omega },$$E31

where $a>0$, $b\ge 0$, $\kappa \ge 0$, $f\in L^2\left(\mathrm{\Omega }\right)$ and $g\in C^0\left(\mathrm{\partial \Omega }\right)$. We seek a weak solution $u$ in

$$V=H^1\left(\mathrm{\Omega }\right)=\left\{v\in L^2\left(\mathrm{\Omega }\right)|v\text{has a weak derivative and}{‖v‖}_{L^2\left(\mathrm{\Omega }\right)}+{‖\nabla v‖}_{L^2\left(\mathrm{\Omega }\right)}<\infty \right\}.$$

In order to derive the weak formulation, we multiply (31) with $v\in V$, integrate over $\mathrm{\Omega }$ and use Green’s formula to obtain

$$\begin{array}{cc}{\int }_{\mathrm{\Omega }}\mathit{\text{fvdxdy}}& =-{\int }_{\mathrm{\Omega }}v\nabla \cdot \left(a\nabla u\right)\mathit{\text{dxdy}}+{\int }_{\mathrm{\Omega }}\mathit{\text{buvdxdy}}\\ & ={\int }_{\mathrm{\Omega }}a\nabla u\cdot \nabla \mathit{\text{vdxdy}}-{\int }_{\mathrm{\partial \Omega }}v\left(a\nabla u\right)\cdot \mathbf{n}\mathit{ds}+{\int }_{\mathrm{\Omega }}\mathit{\text{buvdxdy}}\\ & ={\int }_{\mathrm{\Omega }}a\nabla u\cdot \nabla \mathit{\text{vdxdy}}+{\int }_{\mathrm{\Omega }}\mathit{\text{buvdxdy}}+{\int }_{\mathrm{\partial \Omega }}\kappa \left(u-g\right)\mathit{vds}.\end{array}$$

We obtain the weak form: Find $u\in V$ such that

$${\int }_{\mathrm{\Omega }}a\nabla u\cdot \nabla \mathit{\text{vdxdy}}+{\int }_{\mathrm{\Omega }}\mathit{\text{buvdxdy}}+{\int }_{\mathrm{\partial \Omega }}\mathit{\text{κuvds}}={\int }_{\mathrm{\Omega }}\mathit{\text{fvdxdy}}+{\int }_{\mathrm{\partial \Omega }}\mathit{\text{κgvds}},v\in V.$$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 $u_h\in V_h$ such that

$${\int }_{\mathrm{\Omega }}a\nabla u_h\cdot \nabla v_h\mathit{\text{dxdy}}+{\int }_{\mathrm{\Omega }}{\mathit{bu}}_h{v}_h\mathit{\text{dxdy}}+{\int }_{\mathrm{\partial \Omega }}\kappa u_h{v}_h\mathit{ds}={\int }_{\mathrm{\Omega }}{\mathit{fv}}_h\mathit{\text{dxdy}}+{\int }_{\mathrm{\partial \Omega }}\kappa {\mathit{gv}}_h\mathit{ds},v_h\in V_h,$$E33

where $V_h=\left\{v\in V{\left|v\right|}_T\in P^1\left(T\right),\forall T\in {\mathcal{T}}_h\right\}$.

Implementation: Let $a=1$ and $b=g=0$. Substituting $u_h={\sum }_{j=1}^N{c}_j{\varphi }_j$ into (33) and picking $v_h={\varphi }_i$, we obtain

$$\sum _{j=1}^N{c}_j\left({\int }_{\mathrm{\Omega }}\nabla {\varphi }_j\cdot \nabla {\varphi }_i\mathit{\text{dxdy}}+{\int }_{\mathrm{\partial \Omega }}{\mathit{\kappa \varphi }}_j{\varphi }_i\mathit{ds}\right)={\int }_{\mathrm{\Omega }}f{\varphi }_i\mathit{\text{dxdy}},i=1,2,\dots ,N.$$

This gives us the system $\left(A+R\right)\mathbf{c}=\mathbf{b},$ where $\mathbf{c}={\left[c_1,c_2,\dots ,c_N\right]}^t\in {\mathbb{R}}^N$ is the unknown vector and the entries of $A\in {\mathbb{R}}^{N\times N}$, $R\in {\mathbb{R}}^{N\times N}$, and $\mathbf{b}\in {\mathbb{R}}^N$ are given by

$$a_{\mathit{ij}}={\int }_{\mathrm{\Omega }}\nabla {\varphi }_j\cdot \nabla {\varphi }_i\mathit{\text{dxdy}},r_{\mathit{ij}}={\int }_{\mathrm{\partial \Omega }}{\mathit{\kappa \varphi }}_j{\varphi }_i\mathit{ds},b_i={\int }_{\mathrm{\Omega }}f{\varphi }_i\mathit{\text{dxdy}},i,j=1,2,\dots ,N.$$

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

$$a_{\mathit{ij}}^T={\int }_{\mathrm{\Omega }}\nabla {\varphi }_j\cdot \nabla {\varphi }_i\mathit{\text{dxdy}},i,j=\mathrm{1,2,3}.$$

where $T$ is an arbitrary triangle with vertices $N_i=\left(x_i,y_i\right)$ and ${\varphi }_i$ are the hat functions $i.e.,$${\varphi }_j\left(N_i\right)={\delta }_{\mathit{ij}}$. Let ${\varphi }_i\left(x,y\right)={\alpha }_i+{\beta }_{i}x+{\gamma }_{i}y$, for $i=\mathrm{1,2,3}$. Then, we compute ${\alpha }_i,{\beta }_i,{\gamma }_i$ by

$$\left[\begin{array}{ccc}1& x_1& y_1\\ 1& x_2& y_2\\ 1& x_3& y_3\end{array}\right]\left[\begin{array}{c}{\alpha }_1\\ {\beta }_1\\ {\gamma }_1\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right],\left[\begin{array}{ccc}1& x_1& y_1\\ 1& x_2& y_2\\ 1& x_3& y_3\end{array}\right]\left[\begin{array}{c}{\alpha }_2\\ {\beta }_2\\ {\gamma }_2\end{array}\right]=\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right],\left[\begin{array}{ccc}1& x_1& y_1\\ 1& x_2& y_2\\ 1& x_3& y_3\end{array}\right]\left[\begin{array}{c}{\alpha }_3\\ {\beta }_3\\ {\gamma }_3\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right].$$

In general we have $B{\alpha }_i={\mathbf{e}}_i$ for $i=\mathrm{1,2,3}$, where

$$B=\left[\begin{array}{ccc}1& x_1& y_1\\ 1& x_2& y_2\\ 1& x_3& y_3\end{array}\right],{\alpha }_i=\left[\begin{array}{c}{\alpha }_i\\ {\beta }_i\\ {\gamma }_i\end{array}\right],{\mathbf{e}}_1=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right],{\mathbf{e}}_2=\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right],{\mathbf{e}}_3=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right].$$

Furthermore, we obviously have $\nabla {\varphi }_i={\left[{\beta }_i,{\gamma }_i\right]}^t$, which gives

$$a_{\mathit{ij}}^T={\int }_{\mathrm{\Omega }}\left({\beta }_i{\beta }_j+{\gamma }_i{\gamma }_j\right)\mathit{dx}=\left({\beta }_i{\beta }_j+{\gamma }_i{\gamma }_j\right)\mid T\mid ,i,j=\mathrm{1,2,3}.$$

Assembly of boundary matrix$R$: Let ${\mathrm{\Gamma }}_h^{\mathit{\text{out}}}$ denote the set of boundary edges of the triangulation, $i.e.$${\mathrm{\Gamma }}_h^{\mathit{\text{out}}}=\left\{E|E=T\cap \mathrm{\partial \Omega },\mathit{\text{for}}T\in {\mathcal{T}}_h\right\}$. Assume that $\kappa$ is constant on $E$. For an edge $E\in {\mathrm{\Gamma }}_h^{\mathit{\text{out}}}$, we define $R^E\in {\mathbb{R}}^{2\times 2}$ by the entries

$$r_{\mathit{ij}}^E={\int }_E{\mathit{\kappa \varphi }}_j{\varphi }_i\mathit{ds}=\frac{\kappa }{6}\left(1+{\delta }_{\mathit{ij}}\right)\mid E\mid ,i,j=1,2,$$

where $\mid E\mid$ is the length of $E$ and ${\delta }_{\mathit{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 $T\in {\mathcal{T}}_h$

$$b_i^T={\int }_{T}f{\varphi }_i\mathit{\text{dxdy}}\approx \frac{\mid T\mid }{3}f\left(N_i\right),i=1,2,\dots ,N.$$

Given $A$, $R$ and $\mathbf{b}$, we can solve $\left(A+R\right)\mathbf{c}=\mathbf{b}$ and write $u_h={\sum }_{j=1}^N{c}_j{\varphi }_j$.

3.7 The Dirichlet problem

Consider the following Dirichlet Problem: Find $u$ such that

$$-\mathrm{\Delta }u=f,\text{in}\mathrm{\Omega },u=g,\mathrm{on}\mathrm{\partial \Omega },$$E34

where $f\in L^2\left(\mathrm{\Omega }\right)$ and $g\in C^0\left(\mathrm{\partial \Omega }\right)$. We seek a weak solution $u$ in $V_g=\left\{v\in V{\left|v\right|}_{\mathrm{\partial \Omega }}=g\right\}.$ Multiplying (34) by a test function $v\in V_0$ and integrating over $\mathrm{\Omega }$, we get

$${\int }_{\mathrm{\Omega }}\mathit{\text{fvdxdy}}=-{\int }_{\mathrm{\Omega }}v\mathrm{\Delta }\mathit{\text{udxdy}}={\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{\text{vdxdy}}-{\int }_{\mathrm{\partial \Omega }}v\nabla u\cdot \mathbf{n}\mathit{ds}={\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{\text{vdxdy}}.$$

So the weak problem reads: Find $u\in V_g$ such that

$${\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{\text{vdxdy}}={\int }_{\mathrm{\Omega }}\mathit{\text{fvdxdy}},v\in V_0.$$

Assume that $g$ is piecewise linear on $\mathrm{\partial \Omega }$ with respect to the triangulation. Then the FE method in 2D is defined as follows: Find $u_h\in V_{h,g}=\left\{v\in V_h{\left|v\right|}_{\mathrm{\partial \Omega }}=g\right\}$ such that

$${\int }_{\mathrm{\Omega }}\nabla u_h\cdot \nabla v_h\mathit{\text{dxdy}}={\int }_{\mathrm{\Omega }}{\mathit{fv}}_h\mathit{\text{dxdy}},v_h\in V_{h,0}.$$

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

$$\left[\begin{array}{cc}{A}_{0,0}& A_{0,g}\\ A_{g,0}& A_{0,g}\end{array}\right]\left[\begin{array}{c}{\mathbf{c}}_0\\ {\mathbf{c}}_1\end{array}\right]=\left[\begin{array}{c}{\mathbf{b}}_0\\ {\mathbf{b}}_1\end{array}\right],$$

where $A_{0,0}\in {\mathbb{R}}^{\left(N-J\right)\times \left(N-J\right)}$, $A_{g,g}\in {\mathbb{R}}^{J\times J}$, $A_{0,g}\in {\mathbb{R}}^{\left(N-J\right)\times J}$, $A_{g,0}\in {\mathbb{R}}^{J\times \left(N-J\right)}$. Note that ${\mathbf{c}}_1\in {\mathbb{R}}^J$ is known (it contains the values of $g$ in the boundary nodes). We can therefore solve the simplified problem reading: find ${\mathbf{c}}_0\in {\mathbb{R}}^{N-J}$ with $A_{0,0}{\mathbf{c}}_0={\mathbf{b}}_0-A_{0,g}{\mathbf{c}}_1$.

3.8 The Neumann problem

Consider the following Neumann Problem: Find $u$ such that

$$-\mathrm{\Delta }u=f,\text{in}\mathrm{\Omega },\nabla u\cdot \mathbf{n}=g,\mathrm{on}\mathrm{\partial \Omega },$$E35

where $f\in L^2\left(\mathrm{\Omega }\right)$ and $g\in C^0\left(\mathrm{\partial \Omega }\right)$. Let us try to seek a solution to this problem in the space $V=\left\{v|{‖v‖}_{L^2\left(\mathrm{\Omega }\right)}+{‖\nabla v‖}_{L^2\left(\mathrm{\Omega }\right)}<\infty \right\}$. Multiplying (35) by a test function $v\in V$, integrating over $\mathrm{\Omega }$, and using Green’s formula, we get

$${\int }_{\mathrm{\Omega }}\mathit{\text{fvdxdy}}=-{\int }_{\mathrm{\Omega }}v\mathrm{\Delta }\mathit{\text{udxdy}}={\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{\text{vdxdy}}-{\int }_{\mathrm{\partial \Omega }}v\nabla u\cdot \mathbf{n}\mathit{ds}={\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{\text{vdxdy}}-{\int }_{\mathrm{\partial \Omega }}\mathit{\text{vgds}}.$$

Thus, the variational formulation reads: find $u\in V$ such that

$${\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{\text{vdxdy}}-{\int }_{\mathrm{\partial \Omega }}\mathit{\text{vgds}}={\int }_{\mathrm{\Omega }}\mathit{\text{fvdxdy}},\forall v\in V.$$

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

$$0={\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla 1\mathit{\text{dxdy}}={\int }_{\mathrm{\Omega }}\mathit{\text{fdxdy}}+{\int }_{\mathrm{\partial \Omega }}\mathit{gds}.$$

Therefore we need to assume the following compatibility condition

$${\int }_{\mathrm{\Omega }}\mathit{\text{fdxdy}}+{\int }_{\mathrm{\partial \Omega }}\mathit{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 $c\in \mathbb{R}$. To fix this constant and obtain a unique solution a common trick is to impose the additional constraint ${\int }_{\mathrm{\Omega }}\mathit{\text{udxdy}}=0$. We therefore define the weak solution space

$$\widehat{V}=\left\{v\in V|{\int }_{\mathrm{\Omega }}\mathit{\text{vdxdy}}=0\right\},$$

which contains only functions with a zero mean value. This is a called a quotient space. This space guarantees a unique weak solution (with weak formulation as usual with test functions in $V$). So the weak problem reads: Find $u\in \widehat{V}$ such that

$${\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{\text{vdxdy}}-{\int }_{\mathrm{\partial \Omega }}\mathit{\text{vgds}}={\int }_{\mathrm{\Omega }}\mathit{\text{fvdxdy}},\forall v\in V.$$

Now, the FE method takes the form: find $u_h\in {\widehat{V}}_h\subset \widehat{V}$ such that

$${\int }_{\mathrm{\Omega }}\nabla u_h\cdot \nabla v_h\mathit{\text{dxdy}}-{\int }_{\mathrm{\partial \Omega }}{v}_h\mathit{gds}={\int }_{\mathrm{\Omega }}{\mathit{fv}}_h\mathit{\text{dxdy}},\forall v_h\in {\widehat{V}}_h,$$

where ${\widehat{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 $\mathrm{\Omega }$ is a domain in either ${\mathbb{R}}^2$ or ${\mathbb{R}}^3$ and assume that $\mathrm{\partial \Omega }$ has been partitioned into two disjoint sets: $\mathrm{\partial \Omega }={\mathrm{\Gamma }}_1\cup {\mathrm{\Gamma }}_2$. We consider the following BVP:

$$-\nabla \cdot \left(\kappa \left(\mathbf{x}\right)\nabla u\right)=f\left(\mathbf{x}\right),\mathbf{x}\in \mathrm{\Omega },u=0,\mathbf{x}\in {\mathrm{\Gamma }}_1,\nabla u\cdot \mathbf{n}=0,\mathbf{x}\in {\mathrm{\Gamma }}_2,$$E36

where $f\in L^2\left(\mathrm{\Omega }\right)$. 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

$$\widehat{V}=\left\{v\in C^2\left(\overline{\mathrm{\Omega }}\right):v\left(\mathbf{x}\right)=0,\mathbf{x}\in {\mathrm{\Gamma }}_1\right\}.$$

Multiplying (36) by a test function $v\in \widehat{V}$ and integrating over $\mathrm{\Omega }$, we get

$$\begin{array}{cc}{\int }_{\mathrm{\Omega }}\mathit{\text{fvdxdy}}& =-{\int }_{\mathrm{\Omega }}v\nabla \cdot \left(\kappa \left(\mathbf{x}\right)\nabla u\right)\mathit{\text{dxdy}}={\int }_{\mathrm{\Omega }}\kappa \left(\mathbf{x}\right)\nabla u\cdot \nabla \mathit{\text{vdxdy}}-{\int }_{\mathrm{\partial \Omega }}\kappa \left(\mathbf{x}\right)v\nabla u\cdot \mathbf{n}\mathit{ds}\\ & ={\int }_{\mathrm{\Omega }}\kappa \left(\mathbf{x}\right)\nabla u\cdot \nabla \mathit{\text{vdxdy}}-{\int }_{{\mathrm{\Gamma }}_1}\kappa \left(\mathbf{x}\right)v\nabla u\cdot \mathbf{n}\mathit{ds}-{\int }_{{\mathrm{\Gamma }}_2}\kappa \left(\mathbf{x}\right)v\nabla u\cdot \mathbf{n}\mathit{ds}\\ & ={\int }_{\mathrm{\Omega }}\kappa \left(\mathbf{x}\right)\nabla u\cdot \nabla \mathit{\text{vdxdy}},\end{array}$$

since $v=0$ on ${\mathrm{\Gamma }}_1$ and $\nabla u\cdot \mathbf{n}$ on ${\mathrm{\Gamma }}_1$. Thus the weak form of (36) is: Find $u\in \widehat{V}$ such that

$${\int }_{\mathrm{\Omega }}\kappa \left(\mathbf{x}\right)\nabla u\cdot \nabla \mathit{\text{vdxdy}}={\int }_{\mathrm{\Omega }}\mathit{\text{fvdxdy}},v\in \widehat{V}.$$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 $\widehat{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 ${\mathrm{\Gamma }}_1$, that is, that ${\mathrm{\Gamma }}_1$ includes its “endpoints”). We can then choose the approximating subspace of $\widehat{V}$ as follows:

$$V_h=\left\{v\in C\left(\overline{\mathrm{\Omega }}\right):\text{v is linear on}{\mathcal{T}}_{\mathrm{h}},v\left(\mathbf{z}\right)=0\text{for all nodes}\mathbf{z}\in {\mathrm{\Gamma }}_1\right\}.$$

A basis for $V_h$ is formed by including all basis functions corresponding to interior boundary nodes that do not belong to ${\mathrm{\Gamma }}_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

$$-\mathrm{\Delta }u=f\left(\mathbf{x}\right),\mathbf{x}\in \mathrm{\Omega }=\left(0,a\right)\times \left(0,b\right),u\left(\mathbf{x}\right)=g\left(\mathbf{x}\right)=\left\{\begin{array}{ll}{g}_1\left(\mathbf{x}\right),& \mathbf{x}\in {\mathrm{\Gamma }}_1,\\ g_2\left(\mathbf{x}\right),& \mathbf{x}\in {\mathrm{\Gamma }}_2,\\ g_3\left(\mathbf{x}\right),& \mathbf{x}\in {\mathrm{\Gamma }}_3,\\ g_4\left(\mathbf{x}\right),& \mathbf{x}\in {\mathrm{\Gamma }}_4,\end{array}\right.$$E38

where ${\mathrm{\Gamma }}_1$, ${\mathrm{\Gamma }}_2$, ${\mathrm{\Gamma }}_3$, and ${\mathrm{\Gamma }}_4$ are, respectively, the bottom, right, top, and left boundary edges of the rectangular domain $\mathrm{\Omega }=\left(0,a\right)\times \left(0,b\right)$. We will assume that the boundary data are continuous, so

$$g_1\left(0\right)=g_4\left(0\right),g_1\left(a\right)=g_2\left(0\right),g_2\left(b\right)=g_3\left(a\right),g_3\left(0\right)=g_4\left(b\right).$$

Suppose we find a function $w$ defined on $\overline{\mathrm{\Omega }}$ and satisfying $w\left(\mathbf{x}\right)=g\left(\mathbf{x}\right)$ for all $\mathbf{x}\in \mathrm{\partial \Omega }$. We then define $v=u-w$ and note that

$$-\mathrm{\Delta }v=-\mathrm{\Delta }u+\mathrm{\Delta }w=f\left(\mathbf{x}\right)+\mathrm{\Delta }w=\widehat{f}\left(\mathbf{x}\right),$$

and $v\left(\mathbf{x}\right)=u\left(\mathbf{x}\right)-w\left(\mathbf{x}\right)=0$ for all $\mathbf{x}\in \mathrm{\partial \Omega }$. We can then solve

$$-\mathrm{\Delta }v=\widehat{f}\left(\mathbf{x}\right),\mathbf{x}\in \mathrm{\Omega },v\left(\mathbf{x}\right)=0,\mathbf{x}\in \mathrm{\partial \Omega }.$$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 $q\left(x,y\right)=c_0+c_{1}x+c_{2}y+c_3\mathit{xy}$, which assumes the desired boundary values at the corners:

$$q\left(0,0\right)=g_1\left(0\right)=g_4\left(0\right),q\left(a,0\right)=g_1\left(a\right)=g_2\left(0\right),q\left(a,b\right)=g_2\left(b\right)=g_3\left(a\right),q\left(0,b\right)=g_3\left(0\right)=g_4\left(b\right).$$

A direct calculation shows that

$$\begin{gathered} c_0=g_1\left(0\right),c_1=\frac{g_1\left(a\right)-g_1\left(0\right)}{a},c_2=\frac{g_4\left(b\right)-g_4\left(0\right)}{b}, \\ {c}_3=\frac{g_2\left(b\right)+g_1\left(0\right)-g_1\left(a\right)-g_4\left(b\right)}{\mathit{ab}}. \end{gathered}$$

We then define

$$h\left(\mathbf{x}\right)=\left\{\begin{array}{ll}{h}_1\left(x\right)=g_1\left(x\right)-\left(g_1\left(0\right)+\frac{g_1\left(a\right)-g_1\left(0\right)}{a}x\right),& \mathbf{x}\in {\mathrm{\Gamma }}_1,\\ h_2\left(y\right)=g_2\left(y\right)-\left(g_2\left(0\right)+\frac{g_2\left(b\right)-g_2\left(0\right)}{b}y\right),& \mathbf{x}\in {\mathrm{\Gamma }}_2,\\ h_3\left(x\right)=g_3\left(x\right)-\left(g_3\left(0\right)+\frac{g_3\left(a\right)-g_3\left(0\right)}{a}x\right),& \mathbf{x}\in {\mathrm{\Gamma }}_3,\\ h_4\left(y\right)=g_4\left(y\right)-\left(g_4\left(0\right)+\frac{g_4\left(b\right)-g_4\left(0\right)}{b}y\right),& \mathbf{x}\in {\mathrm{\Gamma }}_4.\end{array}\right.$$

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

$$h_1\left(0\right)=h_1\left(a\right)=h_2\left(0\right)=h_2\left(b\right)=h_3\left(0\right)=h_3\left(a\right)=h_4\left(0\right)=h_4\left(b\right)=0.$$

Finally, we define

$$w\left(x,y\right)=\left(c_0+c_{1}x+c_{2}y+c_3\mathit{xy}\right)+\left(h_1\left(x\right)+\frac{h_3\left(x\right)-h_1\left(x\right)}{b}y\right)+\left(h_4\left(y\right)+\frac{h_2\left(y\right)-h_4\left(y\right)}{a}x\right).$$

The reader should notice how the second term interpolates between the boundary values on ${\mathrm{\Gamma }}_1$ and ${\mathrm{\Gamma }}_3$, while the third term interpolates between the boundary values on ${\mathrm{\Gamma }}_2$ and ${\mathrm{\Gamma }}_4$. In order for these two terms not to interfere with each other, it is necessary that the boundary data be zero at the corners. It was for this reason that we transformed the $g_i$ ′s into the $h_i$ ′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

$$-\mathrm{\Delta }u=f\left(\mathbf{x}\right),\mathbf{x}\in \mathrm{\Omega }=\left(0,a\right)\times \left(0,b\right),\mathbf{n}\cdot \nabla u\left(\mathbf{x}\right)=g\left(\mathbf{x}\right)=\left\{\begin{array}{ll}{g}_1\left(\mathbf{x}\right),& \mathbf{x}\in {\mathrm{\Gamma }}_1,\\ g_2\left(\mathbf{x}\right),& \mathbf{x}\in {\mathrm{\Gamma }}_2,\\ g_3\left(\mathbf{x}\right),& \mathbf{x}\in {\mathrm{\Gamma }}_3,\\ g_4\left(\mathbf{x}\right),& \mathbf{x}\in {\mathrm{\Gamma }}_4,\end{array}\right.$$E40

where ${\mathrm{\Gamma }}_1$, ${\mathrm{\Gamma }}_2$, ${\mathrm{\Gamma }}_3$, and ${\mathrm{\Gamma }}_4$ are, respectively, the bottom, right, top, and left boundary edges of the rectangular domain $\mathrm{\Omega }=\left(0,a\right)\times \left(0,b\right)$. We first note that this is equivalent to

$$-u_y\left(\mathbf{x}\right)=g_1\left(x\right),\mathbf{x}\in {\mathrm{\Gamma }}_1,u_x\left(\mathbf{x}\right)=g_2\left(y\right),\mathbf{x}\in {\mathrm{\Gamma }}_2,u_y\left(\mathbf{x}\right)=g_3\left(x\right),\mathbf{x}\in {\mathrm{\Gamma }}_3,-u_x\left(\mathbf{x}\right)=g_4\left(y\right),\mathbf{x}\in {\mathrm{\Gamma }}_4.$$

We make the following observation: If there is a twice-continuously differentiable function $u$ satisfying the given Neumann conditions, then, since $u_{\mathit{xy}}=u_{\mathit{yx}}$, we have

$$-u_{\mathit{xy}}\left(x,0\right)=g_1^{\prime }\left(x\right),-u_{\mathit{yx}}\left(0,y\right)=g_{4^{\prime }}\left(y\right),$$

which together imply that $g_1^{\prime }\left(0\right)=g_{4^{\prime }}\left(0\right)$. By similar reasoning, we have all of the following conditions:

$$g_1^{\prime }\left(0\right)=g_{4^{\prime }}\left(0\right),g_1^{\prime }\left(0\right)=g_{4^{\prime }}\left(0\right),-g_1^{\prime }\left(a\right)=g_2^{\prime }\left(0\right),g_2^{\prime }\left(b\right)=g_3^{\prime }\left(a\right).$$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

$$\psi \left(x\right)=-\mathit{\alpha x}+\frac{\alpha +\beta }{2a}{x}^2\text{satisfies}{\psi }^{\prime }\left(0\right)=-\alpha ,{\psi }^{\prime }\left(a\right)=\beta .$$E42

The first step is to transform the boundary data $g_l\left(x\right)$ to a function $h_1\left(x\right)$ satisfying $h_{1^{\prime }}\left(0\right)=h_{1^{\prime }}\left(a\right)=0$, and similarly for $g_2$, $g_3$, $g_4$ and $h_2$, $h_3$, $h_4$. 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 $q\left(x,y\right)$ satisfying the conditions

$$u_{\mathit{xy}}\left(0,0\right)=-g_1^{\prime }\left(0\right),u_{\mathit{xy}}\left(a,0\right)=-g_1^{\prime }\left(a\right),u_{\mathit{xy}}\left(0,b\right)=-g_3^{\prime }\left(0\right),u_{\mathit{xy}}\left(a,b\right)=-g_2^{\prime }\left(b\right).$$

We can satisfy these conditions with a function of the form $q\left(x,y\right)=c_0\mathit{xy}+c_1{x}^{2}y+c_2{\mathit{xy}}^2+c_3{x}^2{y}^2$. The reader can verify that the necessary coefficients are

$$c_0=-g_1^{\prime }\left(0\right),c_1=\frac{g_1^{\prime }\left(0\right)-g_1^{\prime }\left(a\right)}{2a},c_2=\frac{g_3^{\prime }\left(0\right)+g_1^{\prime }\left(0\right)}{2b},c_3=\frac{g_2^{\prime }\left(b\right)+g_1^{\prime }\left(a\right)-g_3^{\prime }\left(0\right)-g_1^{\prime }\left(0\right)}{4\mathit{ab}}.$$

If $w$ is to satisfy the desired Neumann conditions, then $w-q=h_i$ on ${\mathrm{\Gamma }}_i$, $i=1-4$, where

$$h_1\left(x\right)=g_1\left(x\right)+c_{0}x+c_1{x}^2,h_2\left(y\right)=g_2\left(y\right)-\left(c_0+2{\mathit{ac}}_1\right)y-\left(c_2+2{\mathit{ac}}_3\right)y^2,$$

$$h_3\left(x\right)=g_3\left(x\right)-\left(c_0+2{\mathit{bc}}_2\right)x-\left(c_1+2{\mathit{bc}}_3\right)x^2,h_4\left(y\right)=g_4\left(y\right)+c_{0}y+c_2{y}^2.$$

We can now define $w-q$ by the interpolation described by (42):

$$w\left(x,y\right)=q\left(x,y\right)-h_1\left(x\right)y+\frac{h_3\left(x\right)+h_1\left(x\right)}{2b}{y}^2-h_4\left(y\right)x+\frac{h_2\left(y\right)+h_4\left(y\right)}{2a}{\mathit{yx}}^2.$$

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

3.11 Eigenvalue problem

Consider the following Eigenvalue Problem: Find $\lambda \in \mathbb{R}$ and $u$ such that

$$-\mathrm{\Delta }u=\mathit{\lambda u},\text{in}\mathrm{\Omega },\nabla u\cdot \mathbf{n}=0,\mathrm{on}\mathrm{\partial \Omega }.$$E43

In order to derive the weak formulation, we multiply (43) with $v\in V$, integrate over $\mathrm{\Omega }$ and use Green’s formula to obtain

$$\lambda {\int }_{\mathrm{\Omega }}\mathit{\text{uvdxdy}}=-{\int }_{\mathrm{\Omega }}v\mathrm{\Delta }\mathit{\text{udxdy}}={\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{\text{vdxdy}}-{\int }_{\mathrm{\partial \Omega }}v\nabla u\cdot \mathbf{n}\mathit{ds}={\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{\text{vdxdy}}.$$

We obtain the weak form: Find $u\in V$ such that

$${\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{\text{vdxdy}}=\lambda {\int }_{\mathrm{\Omega }}\mathit{\text{uvdxdy}},v\in V.$$E44

The FE method in 2D is defined as follows: Find ${\lambda }_h\in \mathbb{R}$ and $u_h\in V_h$ such that

$${\int }_{\mathrm{\Omega }}\nabla u_h\cdot \nabla v_h\mathit{\text{dxdy}}={\lambda }_h{\int }_{\mathrm{\Omega }}{u}_h{v}_h\mathit{\text{dxdy}},v_h\in V_h,$$E45

where $V_h=\left\{v\in V{\left|v\right|}_T\in P^1\left(T\right),\forall T\in {\mathcal{T}}_h\right\}$.

Implementation: Substituting $u_h={\sum }_{j=1}^N{c}_j{\varphi }_j$ into (45) and picking $v_h={\varphi }_i$, we obtain

$$\sum _{j=1}^N{c}_j\left({\int }_{\mathrm{\Omega }}\nabla {\varphi }_j\cdot \nabla {\varphi }_i\mathit{\text{dxdy}}-{\lambda }_h{\int }_{\mathrm{\Omega }}{\varphi }_i{\varphi }_j\mathit{\text{dxdy}}\right)=0,i=1,2,\dots ,N.$$

This leads to an algebraic system of the form $A\mathbf{c}={\lambda }_{h}M\mathbf{c}$, $i.e.$ an algebraic eigenvalue problem.

3.12 Error analysis

Consider the following model Problem: Find $u$ such that

$$-\mathrm{\Delta }u=f,\text{in}\mathrm{\Omega },u=0,\mathrm{on}\mathrm{\partial \Omega }.$$

The weak form: Find $u\in V_0$ such that

$${\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{\text{vdxdy}}={\int }_{\mathrm{\Omega }}\mathit{\text{fvdxdy}},v\in V_0.$$

The FE approximation is defined as follows: Find $u_h\in V_{h,0}$ such that

$${\int }_{\mathrm{\Omega }}\nabla u_h\cdot \nabla v_h\mathit{\text{dxdy}}={\int }_{\mathrm{\Omega }}{\mathit{fv}}_h\mathit{\text{dxdy}},v_h\in V_{h,0},$$

where $V_h=\left\{v\in V{\left|v\right|}_T\in P^1\left(T\right),\forall T\in {\mathcal{T}}_h\right\}$. Expressing $u_h={\sum }_{j=1}^N{c}_j{\varphi }_j$ and picking $v_h={\varphi }_i$, we obtain

$$\sum _{j=1}^N{c}_j\left({\int }_{\mathrm{\Omega }}\nabla {\varphi }_j\cdot \nabla {\varphi }_i\mathit{\text{dxdy}}\right)={\int }_{\mathrm{\Omega }}f{\varphi }_i\mathit{\text{dxdy}},i=1,2,\dots ,N.$$

This leads to system of the form $A\mathbf{c}=\mathbf{b},$ where the entries of $A\in {\mathbb{R}}^{N\times N}$ and $\mathbf{b}\in {\mathbb{R}}^N$ are

$$a_{\mathit{ij}}={\int }_{\mathrm{\Omega }}\nabla {\varphi }_j\cdot \nabla {\varphi }_i\mathit{\text{dxdy}},b_i={\int }_{\mathrm{\Omega }}f{\varphi }_i\mathit{\text{dxdy}},i,j=1,2,\dots ,N.$$

Theorem 3.6The stiffness matrix$A$is symmetric and positive definite.

Theorem 3.7 (Galerkin orthogonality)Let$u\in V_0$denote the weak solution and$u_h\in V_{h,0}$the corresponding FE method approximation. Then

$${\int }_{\mathrm{\Omega }}\nabla \left(u-u_h\right)\cdot \nabla v_h\mathit{\text{dxdy}}=0,v_h\in V_{h,0}.$$

Now, let ${\left|\left|\left|v\right|\right|\right|}^2={\int }_{\mathrm{\Omega }}\nabla v\cdot \nabla \mathit{\text{vdxdy}}={\int }_{\mathrm{\Omega }}{\left|\nabla v\right|}^2\mathit{\text{dxdy}}$ be the energy norm on $V_0$.

There are two different kinds of error estimates, a priori estimates, where the error is bounded in terms of the exact solution, and a posteriori error estimates, where the error is bounded in terms of the computed solution. Theorem 3.8 (A priori error bound)Let$u\in V_0$denote the weak solution and$u_h\in V_{h,0}$the corresponding FE method approximation. Then

$$\left|\left|\left|u-u_h\right|\right|\right|\le \left|\left|\left|u-v_h\right|\right|\right|,v_h\in V_{h,0}.$$

Theorem 3.9Let$u\in V_0$denote the weak solution and$u_h\in V_{h,0}$the corresponding FE method approximation. If$u\in C^2\left(\mathrm{\Omega }\right)$, then there exists$C$independent of$h_T$and$u$such that

$${\left|\left|\left|u-u_h\right|\right|\right|}_{L^2\left(\mathrm{\Omega }\right)}^2\le C\sum _{T\in {\mathcal{T}}_h}{h}_T^2{‖D^{2}u‖}_{L^2\left(T\right)}^2.$$

3.13 The FE method for elliptic problems with a convection term

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

$$-\nabla \cdot \left(a\nabla u\right)+\mathbf{b}\cdot \nabla u+\mathit{cu}=f,\text{in}\mathrm{\Omega },u=0,\mathrm{on}\mathrm{\partial \Omega }.$$E46

We seek a weak solution $u$ in $V_0=\left\{v\in V{\left|v\right|}_{\mathrm{\partial \Omega }}=0\right\}$. In order to derive the weak formulation, we multiply (46) with $v\in V_0$, integrate over $\mathrm{\Omega }$ and use Green’s formula to obtain

$$\begin{array}{cc}{\int }_{\mathrm{\Omega }}\mathit{\text{fvdxdy}}& =-{\int }_{\mathrm{\Omega }}v\nabla \cdot \left(a\nabla u\right)\mathit{\text{dxdy}}+{\int }_{\mathrm{\Omega }}v\mathbf{b}\cdot \nabla \mathit{\text{udxdy}}+{\int }_{\mathrm{\Omega }}\mathit{\text{cuvdxdy}}\\ & ={\int }_{\mathrm{\Omega }}a\nabla u\cdot \nabla \mathit{\text{vdxdy}}-{\int }_{\mathrm{\partial \Omega }}v\nabla u\cdot \mathbf{n}\mathit{ds}+{\int }_{\mathrm{\Omega }}v\mathbf{b}\cdot \nabla \mathit{\text{udxdy}}+{\int }_{\mathrm{\Omega }}\mathit{\text{cuvdxdy}}\\ & ={\int }_{\mathrm{\Omega }}a\nabla u\cdot \nabla \mathit{\text{vdxdy}}+{\int }_{\mathrm{\Omega }}v\mathbf{b}\cdot \nabla \mathit{\text{udxdy}}+{\int }_{\mathrm{\Omega }}\mathit{\text{cuvdxdy}}.\end{array}$$

Note that there is no need to apply Green’s formula to ${\int }_{\mathrm{\Omega }}v\mathbf{b}\cdot \nabla \mathit{\text{udxdy}}$. We obtain the weak form: Find $u\in V_0$ such that

$${\int }_{\mathrm{\Omega }}a\nabla u\cdot \nabla \mathit{\text{vdxdy}}+{\int }_{\mathrm{\Omega }}v\mathbf{b}\cdot \nabla \mathit{\text{udxdy}}+{\int }_{\mathrm{\Omega }}\mathit{\text{cuvdxdy}}={\int }_{\mathrm{\Omega }}\mathit{\text{fvdxdy}},v\in V_0.$$

The FE method in 2D is defined as follows: Find $u_h\in V_{h,0}=\left\{v\in V_h{\left|v\right|}_{\mathrm{\partial \Omega }}=0\right\}$ such that

$${\int }_{\mathrm{\Omega }}a\nabla u_h\cdot \nabla v_h\mathit{\text{dxdy}}+{\int }_{\mathrm{\Omega }}{v}_h\mathbf{b}\cdot \nabla u_h\mathit{\text{dxdy}}+{\int }_{\mathrm{\Omega }}{\mathit{cu}}_h{v}_h\mathit{\text{dxdy}}={\int }_{\mathrm{\Omega }}{\mathit{fv}}_h\mathit{\text{dxdy}},v_h\in V_{h,0},$$E47

where $V_h=\left\{v\in V{\left|v\right|}_T\in P^1\left(T\right),\forall T\in {\mathcal{T}}_h\right\}$.

Implementation: Substituting $u_h={\sum }_{j=1}^N{c}_j{\varphi }_j$ into (47) and picking $v_h={\varphi }_i$, we obtain

$$\sum _{j=1}^N{c}_j\left({\int }_{\mathrm{\Omega }}a\nabla {\varphi }_j\cdot \nabla {\varphi }_i\mathit{\text{dxdy}}+{\int }_{\mathrm{\Omega }}{\varphi }_i\mathbf{b}\cdot \nabla {\varphi }_j\mathit{\text{dxdy}}+{\int }_{\mathrm{\Omega }}c{\varphi }_i{\varphi }_j\mathit{\text{dxdy}}\right)={\int }_{\mathrm{\Omega }}f{\varphi }_i\mathit{\text{dxdy}},i=1,2,\dots ,N.$$

This gives us the system $\left(A+B+C\right)\mathbf{c}=\mathbf{d}$, where $\mathbf{c}={\left[c_1,\dots ,c_N\right]}^t\in {\mathbb{R}}^N$ is the unknown vector and the entries of $A,B,C\in {\mathbb{R}}^{N\times N}$ and $\mathbf{d}\in {\mathbb{R}}^N$ are given by

$$a_{\mathit{ij}}={\int }_{\mathrm{\Omega }}a\nabla {\varphi }_j\cdot \nabla {\varphi }_i\mathit{\text{dxdy}},b_{\mathit{ij}}={\int }_{\mathrm{\Omega }}{\varphi }_i\mathbf{b}\cdot \nabla {\varphi }_j\mathit{\text{dxdy}},c_{\mathit{ij}}={\int }_{\mathrm{\Omega }}c{\varphi }_i{\varphi }_j\mathit{\text{dxdy}},d_i={\int }_{\mathrm{\Omega }}f{\varphi }_i\mathit{\text{dxdy}},$$

for $i,j=1,2,\dots ,N$. Note that $B$ is not symmetric, $i.e.$$b_{\mathit{ij}}\ne b_{\mathit{ji}}$.

4. The FE method for the heat equation

Consider the following heat/diffusion problem: Find $u\left(\mathbf{x},t\right)$ such that

$$\dot{u}-\mathrm{\Delta }u=f,\text{in}\mathrm{\Omega }\subset {\mathbb{R}}^2,t\in \left[0,T\right],$$E48

$$u\left(\cdot ,t\right)=0,\mathrm{on}\mathrm{\partial \Omega }\text{and}t\in \left[0,T\right],$$E49

$$u\left(\mathbf{x},0\right)=u_0\left(\mathbf{x}\right),\text{for}\mathbf{x}\in \mathrm{\Omega }\text{and}t=0.$$E50

We seek a weak solution $u$ in $V_0=\left\{v{\left|‖v‖+‖\nabla v‖<\infty ,v\right|}_{\mathrm{\partial \Omega }}=0\right\}$. In order to derive the weak formulation, we multiply (48) with $v\in V_0$, integrate over $\mathrm{\Omega }$ and use Green’s formula to obtain, for $t\in \left[0,T\right]$,

$${\int }_{\mathrm{\Omega }}\mathit{fvd}\mathbf{x}={\int }_{\mathrm{\Omega }}\dot{u}\mathit{vd}\mathbf{x}+{\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{vd}\mathbf{x}-{\int }_{\mathrm{\partial \Omega }}v\nabla u\cdot \mathbf{n}\mathit{ds}={\int }_{\mathrm{\Omega }}\dot{u}\mathit{vd}\mathbf{x}+{\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{vd}\mathbf{x}.$$

The weak form therefore reads: Find $u\left(\cdot ,t\right)\in V_0$ such that for $t>0$

$${\int }_{\mathrm{\Omega }}\dot{u}\mathit{vd}\mathbf{x}+{\int }_{\mathrm{\Omega }}\nabla u\cdot \nabla \mathit{vd}\mathbf{x}={\int }_{\mathrm{\Omega }}\mathit{fvd}\mathbf{x},v\in V_0.$$E51

The semi-discrete FE method in 2D is defined as follows: Find $u_h\left(\cdot ,t\right)\in V_{h,0}=\left\{v\in V_h{\left|v\right|}_{\mathrm{\partial \Omega }}=0\right\}$ such that

$${\int }_{\mathrm{\Omega }}{\dot{u}}_h{v}_{h}d\mathbf{x}+{\int }_{\mathrm{\Omega }}\nabla u_h\cdot \nabla v_{h}d\mathbf{x}={\int }_{\mathrm{\Omega }}{\mathit{fv}}_{h}d\mathbf{x},v_h\in V_{h,0},$$E52

where $V_h=\left\{v\in V{\left|v\right|}_T\in P^1\left(T\right),\forall T\in {\mathcal{T}}_h\right\}$.

Implementation: Substituting $u_h\left(\mathbf{x},t\right)={\sum }_{j=1}^N{c}_j\left(t\right){\varphi }_j\left(\mathbf{x}\right)$ into (52) and choosing $v_h={\varphi }_i$, we obtain

$$\sum _{j=1}^N{\dot{c}}_j{\int }_{\mathrm{\Omega }}{\varphi }_j{\varphi }_{i}d\mathbf{x}+\sum _{j=1}^N{c}_j{\int }_{\mathrm{\Omega }}\nabla {\varphi }_j\cdot \nabla {\varphi }_{i}d\mathbf{x}={\int }_{\mathrm{\Omega }}f{\varphi }_{i}d\mathbf{x},i=1,2,\dots ,N.$$

This gives us the system of ODEs

$$M\dot{\mathbf{c}}\left(t\right)+A\left(t\right)\mathbf{c}\left(t\right)=\mathbf{b}\left(t\right),t\in \left(0,T\right],\mathbf{c}\left(0\right)={\mathbf{c}}_0,$$

where $\mathbf{c}={\left[c_1,c_2,\dots ,c_N\right]}^t={\left[u_h\left(N_1,t\right),\dots ,u_h\left(N_N,t\right)\right]}^t\in {\mathbb{R}}^N$ (here $N_i$ denotes the node that belongs to the basis function ${\varphi }_i$) is the unknown vector and the entries of $M,A\in {\mathbb{R}}^{N\times N}$ and $\mathbf{b}\in {\mathbb{R}}^N$ are given by

$$m_{\mathit{ij}}={\int }_{\mathrm{\Omega }}{\varphi }_i{\varphi }_{j}d\mathbf{x},a_{\mathit{ij}}={\int }_{\mathrm{\Omega }}\nabla {\varphi }_j\cdot \nabla {\varphi }_{i}d\mathbf{x},b_i={\int }_{\mathrm{\Omega }}f{\varphi }_{i}d\mathbf{x},i,j=1,2,\dots ,N.$$

Finally, the system of ODEs can be solved with $e.g.,$ the backward Euler method as follows: Let $0=t_0<t_1<\cdots <t_M=T$ be a discretization, let $k_m=t_m-t_{m-1}$ for $m=1,2,\dots ,M$ be the time step size and let ${\mathbf{c}}^m\approx \mathbf{c}\left(t_m\right)$ for $m=1,2,\dots ,M$ denote corresponding approximations. Then, we can compute ${\mathbf{c}}^m$ using

$$\left(M+k_m{A}_m\right){\mathbf{c}}^m=M{\mathbf{c}}^{m-1}+k_m{\mathbf{b}}_m,m=1,2,\dots ,M,$$

where ${\mathbf{c}}^0$ is obtained from $u_0\left(\mathbf{x}\right)$. We can either use ${\mathbf{c}}^0={\left[c_1\left(0\right),\dots ,c_N\left(0\right)\right]}^t={\left[u_0\left(N_1\right),\dots ,u_0\left(N_N\right)\right]}^t$, or we can let ${\mathbf{c}}^0$ to be the $L^2$-projection of $u_0$. We set $u_h^0={\sum }_{j=1}^N{c}_j^0{\varphi }_j\left(\mathbf{x}\right)$ and solve for $c_j^0$ using

$$\sum _{j=1}^N{c}_j^0{\int }_{\mathrm{\Omega }}{\varphi }_j{\varphi }_{i}d\mathbf{x}={\int }_{\mathrm{\Omega }}{u}_0{\varphi }_{i}d\mathbf{x},i=1,2,\dots ,N.$$

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

$$‖u\left(\cdot ,t\right)‖\le ‖u\left(\cdot ,0\right)‖+{\int }_0^t‖f\left(\cdot ,s\right)‖\mathit{ds},‖u_h^m‖\le ‖u_h^{m-1}‖+k_m‖f_m‖\le ‖u_h^0‖+\sum _{i=1}^m{k}_i‖f_i‖.$$

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 $u\left(\mathbf{x},t\right)$ such that

$$\ddot{u}-\nabla \cdot \left(\epsilon \nabla u\right)=f,\text{in}\mathrm{\Omega }\subset {\mathbb{R}}^2,t\in \left[0,T\right],$$E53

$$\mathbf{n}\cdot \nabla u\left(\cdot ,t\right)=0,\mathrm{on}\mathrm{\partial \Omega }\text{and}t\in \left[0,T\right],$$E54

$$u\left(\mathbf{x},0\right)=u_0\left(\mathbf{x}\right),\dot{u}\left(\mathbf{x},0\right)=v_0\left(\mathbf{x}\right),\text{for}\mathbf{x}\in \mathrm{\Omega }\text{and}t=0,$$E55

where $f$ is a given load, $\epsilon =\epsilon \left(\mathbf{x},t\right)$ is a positive parameter, $u_0$ and $v_0$ are a prescribed initial conditions, and $\mathrm{\Omega }$ is a bounded domain with boundary $\mathrm{\partial \Omega }$ and unit outward normal $\mathbf{n}$.

We seek a weak solution $u$ in $V=H^1\left(\mathrm{\Omega }\right)=\left\{v|‖v‖+‖\nabla v‖<\infty \right\}$. Multiplying the wave Eq. (53) with $v\in V$, integrating over $\mathrm{\Omega }$, and using Green’s formula, we obtain, for $t\in \left[0,T\right]$,

$$\begin{array}{cc}{\int }_{\mathrm{\Omega }}\mathit{fvd}\mathbf{x}={\int }_{\mathrm{\Omega }}\ddot{u}\mathit{vd}\mathbf{x}-{\int }_{\mathrm{\Omega }}v\nabla \cdot \left(\epsilon \nabla u\right)d\mathbf{x}& ={\int }_{\mathrm{\Omega }}\ddot{u}\mathit{vd}\mathbf{x}+{\int }_{\mathrm{\Omega }}\epsilon \nabla u\cdot \nabla \mathit{vd}\mathbf{x}-{\int }_{\mathrm{\partial \Omega }}\mathit{v\epsilon }\nabla u\cdot \mathbf{n}\mathit{ds}\\ & ={\int }_{\mathrm{\Omega }}\ddot{u}\mathit{vd}\mathbf{x}+{\int }_{\mathrm{\Omega }}\epsilon \nabla u\cdot \nabla \mathit{vd}\mathbf{x}.\end{array}$$

The weak form (variational formulation) therefore reads: Find $u\left(\cdot ,t\right)\in V=H^1\left(\mathrm{\Omega }\right)$ such that for all $t>0$

$${\int }_{\mathrm{\Omega }}\ddot{u}\mathit{vd}\mathbf{x}+{\int }_{\mathrm{\Omega }}\epsilon \nabla u\cdot \nabla \mathit{vd}\mathbf{x}={\int }_{\mathrm{\Omega }}\mathit{fvd}\mathbf{x},v\in V.$$E56

Let $V_h=\left\{v\in V{\left|v\right|}_T\in P^1\left(T\right),\forall T\in {\mathcal{T}}_h\right\}\subset V$ be the space of all continuous piecewise linear functions on a triangle mesh of $\mathrm{\Omega }$. The semi-discrete FE method in 2D is defined as follows: Find $u_h\left(\cdot ,t\right)\in V_h$ such that

$${\int }_{\mathrm{\Omega }}{\ddot{u}}_h{v}_{h}d\mathbf{x}+{\int }_{\mathrm{\Omega }}\epsilon \nabla u_h\cdot \nabla v_{h}d\mathbf{x}={\int }_{\mathrm{\Omega }}{\mathit{fv}}_{h}d\mathbf{x},v_h\in V_h.$$E57

Implementation: Substituting $u_h\left(\mathbf{x},t\right)={\sum }_{j=1}^N{c}_j\left(t\right){\varphi }_j\left(\mathbf{x}\right)$ into (57) and choosing $v_h={\varphi }_i$, we obtain

$$\sum _{j=1}^N{\ddot{c}}_j{\int }_{\mathrm{\Omega }}{\varphi }_j{\varphi }_{i}d\mathbf{x}+\sum _{j=1}^N{c}_j{\int }_{\mathrm{\Omega }}\epsilon \nabla {\varphi }_j\cdot \nabla {\varphi }_{i}d\mathbf{x}={\int }_{\mathrm{\Omega }}f{\varphi }_{i}d\mathbf{x},i=1,2,\dots ,N.$$

This gives us the system

$$M\ddot{\mathbf{c}}\left(t\right)+A\left(t\right)\mathbf{c}\left(t\right)=\mathbf{b}\left(t\right),t\in \left(0,T\right],$$E58

where $\mathbf{c}={\left[c_1,\dots ,c_N\right]}^t={\left[u_h\left(N_1,t\right),\dots ,u_h\left(N_N,t\right)\right]}^t\in {\mathbb{R}}^N$ (here $N_i$ denotes the node that belongs to the basis function ${\varphi }_i$) is the unknown vector and the entries of the mass and stiffness matrices $M,A\in {\mathbb{R}}^{N\times N}$ and the load vector $\mathbf{b}\in {\mathbb{R}}^N$ are given by

$$m_{\mathit{ij}}={\int }_{\mathrm{\Omega }}{\varphi }_i{\varphi }_{j}d\mathbf{x},a_{\mathit{ij}}={\int }_{\mathrm{\Omega }}\epsilon \nabla {\varphi }_j\cdot \nabla {\varphi }_{i}d\mathbf{x},b_i={\int }_{\mathrm{\Omega }}f{\varphi }_{i}d\mathbf{x},i,j=1,2,\dots ,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 $\mathbf{d}\left(t\right)=\dot{\mathbf{c}}\left(t\right)$, we get the new coupled system

$$M\dot{\mathbf{c}}\left(t\right)-M\mathbf{d}\left(t\right)=0,M\dot{\mathbf{d}}\left(t\right)+A\left(t\right)\mathbf{c}\left(t\right)=\mathbf{b}\left(t\right),t\in \left(0,T\right].$$

Let $\mathbf{w}={\left[\mathbf{c},\mathbf{d}\right]}^t$ then the system is equivalent to $\widehat{M}\dot{\mathbf{w}}\left(t\right)+\widehat{A}\left(t\right)\mathbf{w}\left(t\right)=\widehat{\mathbf{b}}\left(t\right)$, $t\in \left(0,T\right]$, where

$$\widehat{M}=\left[\begin{array}{cc}M& 0\\ 0& M\end{array}\right],\widehat{A}=\left[\begin{array}{cc}0& -M\\ A& 0\end{array}\right],\widehat{\mathbf{b}}=\left[\begin{array}{c}0\\ \mathbf{b}\end{array}\right].$$

Finally, the system of ODEs can be solved with $e.g.,$ the backward Euler method as follows: Let $0=t_0<t_1<\cdots <t_M=T$ be a discretization, let $k_m=t_m-t_{m-1}$ for $m=1,2,\dots ,M$ be the time step size and let ${\mathbf{w}}^m\approx \mathbf{w}\left(t_m\right)$ for $m=1,2,\dots ,M$ denote corresponding approximations. Then, we can compute ${\mathbf{w}}^m$ using

$$\left(\widehat{M}+k_m{\widehat{A}}_m\right){\mathbf{w}}^m=\widehat{M}{\mathbf{w}}^{m-1}+k_m{\widehat{\mathbf{b}}}_m,m=1,2,\dots ,M,$$

where ${\mathbf{w}}^0$ is obtained from $u_0\left(\mathbf{x}\right)$ and $v_0\left(\mathbf{x}\right)$.

There are several possible choices of initial data. We can either use ${\mathbf{w}}^0={\left[w_1\left(0\right),\dots ,c_{2N}\left(0\right)\right]}^t={\left[u_0\left(N_1\right),\dots ,u_0\left(N_N\right),v_0\left(N_1\right),\dots ,v_0\left(N_N\right)\right]}^t$, or we can let ${\mathbf{w}}^0={\left[{\mathbf{w}}_1^0,{\mathbf{w}}_2^0\right]}^t$, where ${\mathbf{w}}_1^0$ and ${\mathbf{w}}_2^0$ are the $L^2$-projection of $u_0$ and $v_0$, respectively. We set $w_{h,1}^0={\sum }_{j=1}^N{w}_{j,1}^0{\varphi }_j\left(\mathbf{x}\right)$ and $w_{h,2}^0={\sum }_{j=1}^N{w}_{j,2}^0{\varphi }_j\left(\mathbf{x}\right)$ and solve for $w_{j,1}^0,w_{j,2}^0$ using

$$\sum _{j=1}^N{w}_{j,1}^0{\int }_{\mathrm{\Omega }}{\varphi }_j{\varphi }_{i}d\mathbf{x}={\int }_{\mathrm{\Omega }}{u}_0{\varphi }_{i}d\mathbf{x},\sum _{j=1}^N{w}_{j,2}^0{\int }_{\mathrm{\Omega }}{\varphi }_j{\varphi }_{i}d\mathbf{x}={\int }_{\mathrm{\Omega }}{v}_0{\varphi }_{i}d\mathbf{x},i=1,2,\dots ,N.$$

We can also use Crank–Nicolson scheme

$$\left(\widehat{M}+\frac{k_m}{2}{\widehat{A}}_m\right){\mathbf{w}}^m=\left(\widehat{M}-\frac{k_m}{2}{\widehat{A}}_{m-1}\right){\mathbf{w}}^{m-1}+\frac{k_m}{2}\left({\widehat{\mathbf{b}}}_{m-1}+{\widehat{\mathbf{b}}}_m\right)\equiv {\mathbf{g}}_m.$$

Theorem 5.1 (Conservation of energy)If$f=0$, then

$${‖{\dot{u}}_h\left(\cdot ,t\right)‖}_{L^2\left(\mathrm{\Omega }\right)}^2+\epsilon {‖\nabla u_h\left(\cdot ,t\right)‖}_{L^2\left(\mathrm{\Omega }\right)}^2={‖\dot{u}\left(\cdot ,0\right)‖}_{L^2\left(\mathrm{\Omega }\right)}^2+\epsilon {‖\nabla u\left(\cdot ,0\right)‖}_{L^2\left(\mathrm{\Omega }\right)}^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.