**Abstract**

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

**Keywords:** the finite element method, initial-value problems, boundary-value problems, Laplace equation, heat equation, wave equation

### **1. Introduction**

#### **1.1 An overview of the finite element method**

Differential equations arise in many disciplines such as engineering, mathematics, sciences, economics, and many other fields. Unfortunately solutions to differential equations can rarely be expressed by closed formulas and numerical methods are needed to approximate their solutions. There are many numerical methods for approximating the solution to differential equations including the finite difference (FD), finite element (FE), finite volume (FV), spectral, and discontinuous Galerkin (DG) methods. These methods are used when the mathematical equations are too complicated to be solved analytically.

The FE method has become the standard numerical scheme for approximating the solution to many mathematical problems; see [1–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:


COMSOL Multiphysics (known as FEMLAB before 2005) is a commercial FE software package designed to address a wide range of physical phenomena. It is widely used in science and industry for research and development. It excels at modeling almost any multi-physics problem by solving the governing set of PDEs via the FE method. This software package is able to solve one, two and threedimensional problems. It comes with a modern graphical user interface to set up simulation models and can be scripted from Matlab or via its native Java API.

In this chapter, we introduce the FE method for several one-dimensional and two-dimensional model problems. Although the FE method has been extensively used in the field of structural mechanics, it has been successfully applied to solve several other types of engineering problems, such as heat conduction, fluid dynamics, seepage flow, and electric and magnetic fields. These applications prompted mathematicians to use this technique for the solution of complicated problems. For illustration, we will use simple one-dimensional and two-dimensional model problems to introduce the FE method.

#### **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' = f(\boldsymbol{\kappa}), \quad \boldsymbol{\kappa} \in [a, b], \quad u(a) = u\_0. \tag{1}$$

*A Brief Summary of the Finite Element Method for Differential Equations DOI: http://dx.doi.org/10.5772/intechopen.95423*

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

	- by multiplying the ODE in (1) by a test function *v x*ð Þ∈*V*<sup>0</sup> ¼
	- *v*∈*L*<sup>2</sup> ½ � *<sup>a</sup>*, *<sup>b</sup>* : <sup>∥</sup>*v*∥<sup>2</sup> <sup>þ</sup> <sup>∥</sup>*v*<sup>0</sup> <sup>∥</sup><sup>2</sup> <sup>&</sup>lt; <sup>∞</sup>, *v a*ð Þ¼ <sup>0</sup> n o, where <sup>∥</sup>*v*∥<sup>2</sup> <sup>¼</sup> <sup>Ð</sup> *<sup>b</sup> <sup>a</sup> <sup>v</sup>*2ð Þ *<sup>x</sup> dx*,

integrating from *a* to *b*, using integration by parts, and applying *v a*ð Þ¼ 0, to get

$$\int\_{a}^{b} \! f v d\mathbf{x} = \int\_{a}^{b} u' v d\mathbf{x} = -\int\_{a}^{b} u v' d\mathbf{x} + u(b)v(b) - u(a)v(a) = -\int\_{a}^{b} u v' d\mathbf{x} + u(b)v(b).$$


Define the FE space as the set of all continuous piecewise linear polynomials *Vh* ¼ f*<sup>v</sup>* : *<sup>v</sup>*j*Ii* <sup>∈</sup>*P*<sup>1</sup> ð Þ *Ii* , *<sup>i</sup>* <sup>¼</sup> 1, 2, … , *<sup>N</sup>*, *v a*ð Þ¼ <sup>0</sup>g, where *<sup>P</sup>*<sup>1</sup> ð Þ *Ii* is the space of polynomials of degree ≤1 on *Ii*. Functions in *Vh* are linear on each *Ii*, and continuous on the whole interval ½ � *a*, *b* . An example of such a function is shown in **Figure 1**.

We remark that any function *v*∈*Vh* is uniquely determined by its nodal values *v x*ð Þ*<sup>i</sup>* .

4.Construct a set of basis functions based on the triangulation. Since *Vh* has finite dimension, we can find one set of basis functions. A basis for *Vh*

is *ϕ<sup>j</sup>* n o*<sup>N</sup> j*¼0 , where *ϕ<sup>j</sup>* ∈*Vh* are linearly independent. Then

**Figure 1.** *A continuous piecewise linear function v.*

*Vh* <sup>¼</sup> *vh*ð Þ *<sup>x</sup>* <sup>∈</sup>*V*, *vh*ð Þ¼ *<sup>x</sup>* <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>0</sup>*cjϕ<sup>j</sup>* ð Þ *x* n o is the space spanned by the basis functions f g *ϕ<sup>i</sup> N <sup>i</sup>*¼0. The simplest finite dimensional space is the piecewise continuous linear function space defined over the triangulation.

*Vh* ¼ f*vh*ð Þ *x* ∈*V*, *vh*ð Þ *x* is piecewise continuous linear over a, b ½ � with *vh*ð Þ¼ *a* 0g*:*

There are infinite number of sets of basis functions. We should choose a set of basis functions that are simple, have compact (minimum) support (that is, zero almost everywhere except for a small region), and meet the regularity requirement, that is, they have to be continuous, and differentiable except at nodal points. The simplest ones are the so-called hat functions satisfying *ϕi*ð Þ¼ *xi* 1 and *ϕ<sup>i</sup> x <sup>j</sup>* � � <sup>¼</sup> 0 for *<sup>i</sup>* 6¼ *<sup>j</sup>*. The analytic form is (see **Figure 2**)

$$\begin{aligned} \phi\_0(\mathbf{x}) &= \begin{cases} \frac{\mathbf{x}\_1 - \mathbf{x}}{h}, & \mathbf{x} \in I\_1, \\ 0, & \text{else}, \end{cases}, \quad \phi\_N(\mathbf{x}) = \begin{cases} \frac{\mathbf{x} - \mathbf{x}\_{N-1}}{h}, & \mathbf{x} \in I\_N, \\ 0, & \text{else}, \end{cases}, \\\ \phi\_i(\mathbf{x}) &= \begin{cases} \frac{\mathbf{x} - \mathbf{x}\_{i-1}}{h}, & \mathbf{x} \in I\_i, \\ \frac{\mathbf{x}\_{i+1} - \mathbf{x}}{h}, & \mathbf{x} \in I\_{i+1}, \\ 0, & \text{else}. \end{cases} \end{aligned}$$

5.Approximate the exact solution *u* by a continuous piecewise linear function *uh*ð Þ *x* . The FE method consists of finding *uh* ∈*Vh* such that

$$-\int\_{a}^{b} u\_{h}v' d\mathfrak{x} + u\_{h}(b)v(b) = \int\_{a}^{b} fv d\mathfrak{x}, \quad \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 *uh*ð Þ¼ *<sup>x</sup>* <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>0</sup>*cjϕ<sup>j</sup>* ð Þ *x* , where *c*0,*c*1, … ,*cN* are unknowns, and choosing *v* ¼ *ϕi*, *i* ¼ 1, 2, … , *N* to get

**Figure 2.** *A typical hat function ϕ<sup>i</sup> on a mesh. Also shown is the half hat functions ϕ*<sup>0</sup> *and ϕN.*

*A Brief Summary of the Finite Element Method for Differential Equations DOI: http://dx.doi.org/10.5772/intechopen.95423*

since *uh*ð Þ¼ *b cN*. Note that using the hat functions, we have *uh*ð Þ¼ *x*<sup>0</sup> 0 and *uh*ð Þ¼ *xi* P*<sup>N</sup> <sup>j</sup>*¼<sup>0</sup>*cjϕ<sup>j</sup>* ð Þ¼ *xi ciϕi*ð Þ¼ *xi ci* for *i* ¼ 1, 2, … , *N*. Thus, we get the following linear system

$$-\sum\_{j=1}^{N} c\_j \int\_a^b \phi\_j \phi\_i' d\mathbf{x} + c\_N \phi\_i(b) = \int\_a^b \phi\_0 \phi\_i' d\mathbf{x}, \quad i = 1, 2, \dots, N.$$

Finally, we solve the linear system for *c*1, … ,*cN*. We note that for *i* ¼ 1, 2, … , *N* � 1, we have

$$\int\_{a}^{b} \phi\_{i} \phi\_{i}^{\prime} d\mathbf{x} = \int\_{\mathcal{X}\_{i-1}}^{\mathcal{X}\_{i+1}} \phi\_{i} \phi\_{i}^{\prime} d\mathbf{x} = \frac{1}{h\_{i}} \int\_{\mathcal{X}\_{i-1}}^{\mathcal{X}\_{i}} \left(\frac{\boldsymbol{\varkappa} - \boldsymbol{\varkappa}\_{i-1}}{h\_{i}}\right) d\mathbf{x} - \frac{1}{h\_{i}} \int\_{\mathcal{X}\_{i}}^{\mathcal{X}\_{i+1}} \left(\frac{\boldsymbol{\varkappa}\_{i+1} - \boldsymbol{\varkappa}}{h\_{i}}\right) d\mathbf{x} = \mathbf{0} \cdot \boldsymbol{\Lambda}$$

However, for *i* ¼ *N*, we have

$$\begin{split} \int\_{a}^{b} \phi\_{N} \phi\_{N} ' d\mathbf{x} &= \int\_{\begin{subarray}{c} \mathbf{x}\_{N-1} \\ \mathbf{x} \end{subarray}}^{\mathbf{x}\_{N}} \phi\_{N} \phi\_{N} ' d\mathbf{x} = \int\_{\begin{subarray}{c} \mathbf{x}\_{N-1} \\ \mathbf{x}\_{N} \end{subarray}}^{\mathbf{x}\_{N}} \left( \frac{\mathbf{x} - \mathbf{x}\_{N-1}}{h\_{N}} \right) \left( \frac{\mathbf{x} - \mathbf{x}\_{N-1}}{h\_{N}} \right) d\mathbf{x} \\ &= \frac{1}{h\_{N}} \int\_{\begin{subarray}{c} \mathbf{x}\_{N-1} \end{subarray}}^{\mathbf{x}\_{N}} \frac{\mathbf{x} - \mathbf{x}\_{N-1}}{h\_{N}} d\mathbf{x} = \frac{1}{2} . \end{split}$$

Similarly, for *i* ¼ 1, 2, … , *N*, we have

$$\begin{split} \int\_{a}^{b} \phi\_{i-1} \phi\_{i}^{\prime} dx &= \int\_{\mathbf{x}\_{i-1}}^{\mathbf{x}\_{i}} \phi\_{i-1} \phi\_{i}^{\prime} dx = \int\_{\mathbf{x}\_{i-1}}^{\mathbf{x}\_{i}} \left( \frac{\mathbf{x}\_{i} - \mathbf{x}}{h\_{i}} \right) \left( \frac{\mathbf{x} - \mathbf{x}\_{i-1}}{h\_{i}} \right) dx = \frac{1}{h\_{i}} \int\_{\mathbf{x}\_{i-1}}^{\mathbf{x}\_{i}} \frac{\mathbf{x}\_{i} - \mathbf{x}}{h\_{i}} dx = \frac{1}{2}, \\ \int\_{a}^{b} \phi\_{i+1} \phi\_{i}^{\prime} dx &= \int\_{\mathbf{x}\_{i}}^{\mathbf{x}\_{i+1}} \phi\_{i+1} \phi\_{i}^{\prime} dx = \int\_{\mathbf{x}\_{i}}^{\mathbf{x}\_{i+1}} \left( \frac{\mathbf{x} - \mathbf{x}\_{i}}{h\_{i+1}} \right) \left( \frac{\mathbf{x}\_{i+1} - \mathbf{x}}{h\_{i+1}} \right) dx = -\frac{1}{h\_{i+1}} \int\_{\mathbf{x}\_{i}}^{\mathbf{x}\_{i+1}} \frac{\mathbf{x} - \mathbf{x}\_{i}}{h\_{i+1}} dx \\ &= -\frac{1}{2}. \end{split}$$

We next calculate Ð *<sup>b</sup> <sup>a</sup> fϕidx*. Since it depends on *f*, we cannot generally expect to calculate it exactly. However, we can approximate it using a quadrature rule. Using the Trapezoidal rule Ð *<sup>b</sup> <sup>a</sup> f x*ð Þ*dx*<sup>≈</sup> *<sup>b</sup>*�*<sup>a</sup>* <sup>2</sup> ð Þ *f a*ð Þþ *f b*ð Þ and using *ϕi*ð Þ¼ *xi*�<sup>1</sup> *ϕi*ð Þ¼ *xi*þ<sup>1</sup> 0 and *ϕi*ð Þ¼ *xi* 1, we get

$$\int\_{a}^{b} f \phi\_{i} d\mathbf{x} = \int\_{x\_{i-1}}^{x\_{i}} f \phi\_{i} d\mathbf{x} + \int\_{x\_{i}}^{x\_{i+1}} f \phi\_{i} d\mathbf{x} \approx \frac{h\_{i} + h\_{i+1}}{2} f(\mathbf{x}\_{i}), \quad i = 1, 2, \dots, N - 1,$$

$$\int\_{a}^{b} f(\mathbf{x}) \phi\_{N} d\mathbf{x} = \int\_{x\_{N-1}}^{x\_{N}} f(\mathbf{x}) \phi\_{N} d\mathbf{x} \approx \frac{h\_{N}}{2} \left( f(\mathbf{x}\_{N-1}) \phi\_{N}(\mathbf{x}\_{N-1}) + f(\mathbf{x}\_{N}) \phi\_{N}(\mathbf{x}\_{N}) \right) = \frac{h\_{N}}{2} f(\mathbf{x}\_{N}).$$

Thus, we obtain the following linear system of equations

$$
\begin{bmatrix}
0 & \frac{1}{2} & 0 & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots & -\frac{1}{2} & 0 & \frac{1}{2} \\
0 & \cdots & 0 & -\frac{1}{2} & \frac{1}{2}
\end{bmatrix}
\begin{bmatrix}
c\_1 \\ c\_2 \\ \vdots \\ c\_{N-1} \\ c\_N
\end{bmatrix} = 
\begin{bmatrix}
\frac{h\_1 + h\_2}{2} f(\mathbf{x}\_1) \\
\frac{h\_2 + h\_3}{2} f(\mathbf{x}\_2) \\
\vdots \\
\frac{h\_{N-1} + h\_N}{2} f(\mathbf{x}\_{N-1}) \\
\frac{h\_N}{2} f(\mathbf{x}\_N)
\end{bmatrix}.
$$

The determinant of the above matrix is <sup>1</sup> <sup>2</sup>*<sup>N</sup>*. Thus, the system has a unique solution *c*1,*c*2, … ,*cN*.

**Remark 2.1** *Suppose that u a*ð Þ¼ *<sup>u</sup>*0*, then we let uh*ð Þ¼ *<sup>x</sup>* <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>0</sup>*cjϕ<sup>j</sup>* ð Þ *x . Since u*<sup>0</sup> ¼ *uh*ð Þ¼ *x*<sup>0</sup> P*<sup>N</sup> <sup>j</sup>*¼<sup>1</sup>*cjϕ<sup>j</sup>* ð Þ¼ *x*<sup>0</sup> *c*0*ϕ*0ð Þ¼ *x*<sup>0</sup> *c*0*, we only need to find c*1,*c*2, … ,*cN. Choosing v* ¼ *ϕi*, *i* ¼ 1, 2, … , *N, we get the following linear system*

$$-\sum\_{j=1}^{N} c\_j \int\_a^b \phi\_j \phi\_i' d\mathbf{x} + c\_N \phi\_i(b) = \int\_a^b f \phi\_i d\mathbf{x} + u\_0 \int\_a^b \phi\_0 \phi\_i' d\mathbf{x}, \quad i = 1, 2, \dots, N.$$

Finally, we solve the linear system for *c*1, … ,*cN*. We note that Ð *<sup>b</sup> <sup>a</sup> ϕ*0*ϕ*<sup>0</sup> *i dx* ¼ 0 for *i* ¼ 2, … , *N* and

$$\int\_{a}^{b} \phi\_{0} \phi\_{1} ^{\prime} d\boldsymbol{\infty} = \int\_{\boldsymbol{\infty}\_{0}}^{\boldsymbol{\infty}\_{1}} \left( \frac{\boldsymbol{\infty}\_{1} - \boldsymbol{\infty}}{h\_{1}} \right) \left( \frac{\boldsymbol{\infty} - \boldsymbol{\infty}\_{0}}{h\_{1}} \right)^{\prime} d\boldsymbol{\infty} = \frac{1}{h\_{1}} \int\_{\boldsymbol{\infty}\_{0}}^{\boldsymbol{\infty}\_{1}} \frac{\boldsymbol{\infty}\_{1} - \boldsymbol{\infty}}{h\_{1}} d\boldsymbol{\infty} = \frac{1}{2} \dots$$

Following the same steps used for the case *u a*ð Þ¼ 0, we obtain the following linear system of equations

$$
\begin{bmatrix} 0 & \frac{1}{2} & 0 & \cdots & 0 \\ -\frac{1}{2} & 0 & \frac{1}{2} & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & -\frac{1}{2} & 0 & \frac{1}{2} \\ 0 & \cdots & 0 & -\frac{1}{2} & \frac{1}{2} \end{bmatrix} \begin{bmatrix} c\_1 \\ c\_2 \\ \vdots \\ c\_{N-1} \\ c\_N \end{bmatrix} = \begin{bmatrix} \frac{h\_1 + h\_2}{2} f(\mathbf{x}\_1) + \frac{u\_0}{2} \\ \frac{h\_2 + h\_3}{2} f(\mathbf{x}\_2) \\ \vdots \\ \frac{h\_{N-1} + h\_N}{2} f(\mathbf{x}\_{N-1}) \\ \frac{h\_N}{2} f(\mathbf{x}\_N) \end{bmatrix}.$$

#### **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' = f(\kappa, u), \quad \kappa \in [a, b], \quad u(a) = u\_0. \tag{2}$$

The FE method consists of finding *uh* <sup>∈</sup>*Vh* ¼ f*<sup>v</sup>* : *<sup>v</sup>*j*Ii* <sup>∈</sup> *<sup>P</sup>*<sup>1</sup> ð Þ *Ii* , *i* ¼ 1, 2, … , *N*, *v a*ð Þ¼ 0g, such that

$$u\_h(b)v(b) - \int\_a^b u\_h v' d\mathfrak{x} = \int\_a^b f(\mathfrak{x}, u\_h) v d\mathfrak{x}, \quad \forall \ v \in V\_h.$$

Writing *uh*ð Þ¼ *<sup>x</sup>* <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>0</sup>*cjϕ<sup>j</sup>* ð Þ *x* and choosing *v* ¼ *ϕi*, *i* ¼ 1, 2, … , *N*, we get

$$\begin{aligned} c\_N \left( \phi\_i - \int\_a^b \phi\_N \phi\_i' dx \right) - \sum\_{j=0}^{N-1} c\_j \int\_a^b \phi\_j \phi\_i' dx - \int\_a^b f \left( \mathbf{x}, \sum\_{j=0}^N c\_j \phi\_j \right) \phi\_i d\mathbf{x} &= \mathbf{0}, \\ i &= 1, 2, \dots, N, \end{aligned}$$

where *uh*ð Þ¼ *x*<sup>0</sup> *c*<sup>0</sup> ¼ *u*0. Finally, we solve the nonlinear system for *c*1,*c*2, … ,*cN* using *e:g:*, Newton's method for systems of nonlinear equations. The system can be written as *Fi*ð Þ¼ *c*1,*c*2, … ,*cN* 0, *i* ¼ 1, 2, … , *N*, where

*A Brief Summary of the Finite Element Method for Differential Equations DOI: http://dx.doi.org/10.5772/intechopen.95423*

$$\begin{aligned} F\_i &= c\_N \left( \phi\_i - \int\_a^b \phi\_N \phi\_i' d\mathbf{x} \right) - \sum\_{j=0}^{N-1} c\_j \int\_a^b \phi\_j \phi\_i' d\mathbf{x} - \int\_a^b f \left( \mathbf{x}, \sum\_{j=0}^N c\_j \phi\_j \right) \phi\_i d\mathbf{x}, \\\ i &= 1, 2, \dots, N. \end{aligned}$$

Let *<sup>α</sup><sup>i</sup>* <sup>¼</sup> <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>0</sup>*cj* Ð *b <sup>a</sup> ϕjϕ*<sup>0</sup> *i dx* and *<sup>β</sup><sup>i</sup>* <sup>¼</sup> <sup>Ð</sup> *<sup>b</sup> <sup>a</sup> f x*, P*<sup>N</sup> <sup>j</sup>*¼<sup>0</sup>*cjϕ<sup>j</sup>* � �*ϕidx*. Then, for *<sup>i</sup>* <sup>¼</sup> 1, 2, … , *N* � 1,

*α<sup>i</sup>* ¼ *ci*�<sup>1</sup> ð*xi xi*�<sup>1</sup> *ϕi*�1*ϕ*<sup>0</sup> *i dx* þ *ci* ð*xi xi*�<sup>1</sup> *ϕiϕ*<sup>0</sup> *i dx* þ ð*xi*þ<sup>1</sup> *xi ϕiϕ*<sup>0</sup> *i dx* � � <sup>þ</sup> *ci*þ<sup>1</sup> ð*xi*þ<sup>1</sup> *xi ϕi*þ1*ϕ*<sup>0</sup> *i dx* ¼ *ci*�<sup>1</sup> ð*xi xi*�<sup>1</sup> *xi* � *x h*2 *i dx* þ *ci* ð*xi xi*�<sup>1</sup> *x* � *xi*�<sup>1</sup> *h*2 *i dx* � ð*xi*þ<sup>1</sup> *xi xi*þ<sup>1</sup> � *x h*2 *i*þ1 *dx* ! � *ci*þ<sup>1</sup> ð*xi*þ<sup>1</sup> *xi x* � *xi h*2 *i*þ1 *dx* ¼ 1 2 *ci*�<sup>1</sup> þ *ci* 1 2 � 1 2 � � � 1 2 *ci*þ<sup>1</sup> <sup>¼</sup> <sup>1</sup> 2 *ci*�<sup>1</sup> � <sup>1</sup> 2 *ci*þ1, *α<sup>N</sup>* ¼ *cN*�<sup>1</sup> ð*xN xN*�<sup>1</sup> *ϕ<sup>N</sup>*�<sup>1</sup>*ϕ<sup>N</sup>* 0 *dx* þ *cN* ð*xN xN*�<sup>1</sup> *ϕNϕ<sup>N</sup>* 0 *dx* ¼ *cN*�<sup>1</sup> ð*xN xN*�<sup>1</sup> *xN* � *x h*2 *N dx* þ *cN* ð*xN xN*�<sup>1</sup> *x* � *xN*�<sup>1</sup> *h*2 *N dx* <sup>¼</sup> <sup>1</sup> 2 *cN*�<sup>1</sup> þ 1 2 *cN:*

Similarly,

$$\beta\_i = \int\_{\mathbf{x}\_{i-1}}^{\mathbf{x}\_{i+1}} f\left(\mathbf{x}, \sum\_{j=0}^{N} c\_j \phi\_j\right) \phi\_i d\mathbf{x} = \int\_{\mathbf{x}\_{i-1}}^{\mathbf{x}\_i} f\left(\mathbf{x}, \sum\_{j=0}^{N} c\_j \phi\_j\right) \phi\_i d\mathbf{x} + \int\_{\mathbf{x}\_i}^{\mathbf{x}\_{i+1}} f\left(\mathbf{x}, \sum\_{j=0}^{N} c\_j \phi\_j\right) \phi\_i d\mathbf{x} \dots d\mathbf{x}$$

Using Simpson's Rule Ð *<sup>b</sup> <sup>a</sup> f x*ð Þ*dx*<sup>≈</sup> *<sup>b</sup>*�*<sup>a</sup>* <sup>6</sup> *f a*ð Þþ <sup>4</sup>*<sup>f</sup> <sup>a</sup>*þ*<sup>b</sup>* 2 � � <sup>þ</sup> *f b*ð Þ � �, and using *<sup>ϕ</sup>i*ð Þ¼ *xi*�<sup>1</sup> *<sup>ϕ</sup>i*ð Þ¼ *xi*þ<sup>1</sup> 0, *<sup>ϕ</sup>i*ð Þ¼ *xi* 1, <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>0</sup>*cjϕ<sup>j</sup> xi*�<sup>1</sup> <sup>þ</sup> *hi* 2 � � <sup>¼</sup> *ci*�1þ*ci* <sup>2</sup> , *<sup>ϕ</sup><sup>i</sup> xi*�<sup>1</sup> <sup>þ</sup> *hi* 2 � � <sup>¼</sup> <sup>1</sup> 2 , P*<sup>N</sup> <sup>j</sup>*¼<sup>0</sup>*cjϕj*ð Þ¼ *xi ci*, we have, for *<sup>i</sup>* <sup>¼</sup> 1, 2, … , *<sup>N</sup>* � 1,

$$
\beta\_i \approx \frac{h\_i}{3} f\left(\mathbf{x}\_{i-1} + \frac{h\_i}{2}, \frac{c\_{i-1} + c\_i}{2}\right) + \frac{h\_i + h\_{i+1}}{6} f(\mathbf{x}\_i, c\_i) + \frac{h\_{i+1}}{3} f\left(\mathbf{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( \mathfrak{F} \left( \mathfrak{x}\_{N-1} + \frac{h}{2}, \frac{c\_{N-1} + c\_N}{2} \right) + f(\mathfrak{x}\_N, c\_N) \right).
$$

Next, we compute the Jacobian matrix with entries

$$J\_{i,j} = \frac{\partial F\_i}{\partial c\_j} = \int\_a^b \phi\_j \phi\_i' d\mathbf{x} - \int\_a^b f\_u \left(\mathbf{x}, \sum\_{j=0}^N c\_j \phi\_j\right) \phi\_j \phi\_i \, d\mathbf{x} = a\_{i,j} - b\_{i,j}, \ \ i = 1, 2, \dots, N.$$

We already computed the entries *ai*,*<sup>j</sup>* as

$$a\_{i,i-1} = \int\_{a}^{b} \phi\_{i-1} \phi\_i' d\mathbf{x} = \frac{1}{2}, \quad a\_{i,i} = \int\_{a}^{b} \phi\_i \phi\_i' d\mathbf{x} = 0, \quad i = 1, 2, \dots, N - 1,$$

$$a\_{N,N} = \int\_{a}^{b} \phi\_N \phi\_N' d\mathbf{x} = \frac{1}{2}, \quad a\_{i,i+1} = \int\_{a}^{b} \phi\_{i+1} \phi\_i' d\mathbf{x} = -\frac{1}{2}.$$

Using Simpson's Rule, we get

*bi*,*i*�<sup>1</sup> ¼ ð*xi xi*�<sup>1</sup> *ϕ<sup>i</sup>*�<sup>1</sup>*ϕ<sup>i</sup> f <sup>u</sup> x*, X *N j*¼0 *cjϕ<sup>j</sup>* !*dx*<sup>≈</sup> *hi* <sup>6</sup> *<sup>f</sup> <sup>u</sup> xi*�<sup>1</sup> <sup>þ</sup> *hi* 2 , *ci*�<sup>1</sup> þ *ci* 2 � �, *bi*,*i*þ<sup>1</sup> ¼ ð*xi*þ<sup>1</sup> *xi ϕ<sup>i</sup>*þ<sup>1</sup>*ϕ<sup>i</sup> f <sup>u</sup> x*, X *N j*¼0 *cjϕ<sup>j</sup>* !*dx*<sup>≈</sup> *hi*þ<sup>1</sup> <sup>6</sup> *<sup>f</sup> <sup>u</sup> xi* <sup>þ</sup> *hi*þ<sup>1</sup> 2 , *ci* þ *ci*þ<sup>1</sup> 2 � �, *bi*,*<sup>i</sup>* ¼ ð*xi xi*�<sup>1</sup> *ϕ*2 *<sup>i</sup> f <sup>u</sup> x*, X *N j*¼0 *cjϕ<sup>j</sup>* !*dx* þ ð*xi*þ<sup>1</sup> *xi ϕ*2 *<sup>i</sup> f <sup>u</sup> x*, X *N j*¼0 *cjϕ<sup>j</sup>* !*dx* <sup>≈</sup> *hi* <sup>6</sup> *<sup>f</sup> <sup>u</sup> xi*�<sup>1</sup> <sup>þ</sup> *hi* 2 , *ci*�<sup>1</sup> þ *ci* 2 � � þ *hi* þ *hi*þ<sup>1</sup> <sup>6</sup> *<sup>f</sup> <sup>u</sup> xi* ð Þþ ,*ci hi*þ<sup>1</sup> <sup>6</sup> *<sup>f</sup> <sup>u</sup> xi* <sup>þ</sup> *hi*þ<sup>1</sup> 2 , *ci* þ *ci*þ<sup>1</sup> 2 � �, *bN*,*<sup>N</sup>* ¼ ð*xN xN*�<sup>1</sup> *ϕ*2 *<sup>N</sup> f <sup>u</sup> x*, X *N j*¼0 *cjϕ<sup>j</sup>* !*dx*<sup>≈</sup> *hN* <sup>6</sup> *<sup>f</sup> <sup>u</sup> xN*�<sup>1</sup> <sup>þ</sup> *h* 2 , *cN*�<sup>1</sup> þ *cN* 2 � � þ *f <sup>u</sup>*ð Þ *xN*,*cN* � �*:*

#### **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* ∈*C*<sup>2</sup> ½ � *a*, *b* such that

$$-u'' + q(x)u = f(x), \quad x \in \Omega = (a, b), \quad u(a) = u(b) = 0,\tag{3}$$

where *u* : Ω ¼ ½ �! *a*, *b* is the sought solution, *q x*ð Þ≥ 0 is a continuous function on ½ � *<sup>a</sup>*, *<sup>b</sup>* , and *<sup>f</sup>* <sup>∈</sup>*L*<sup>2</sup> ½ � *a*, *b* . Under these assumptions, (3) has a unique solution *u* ∈*C*<sup>2</sup> ½ � *a*, *b* . For general *q x*ð Þ, it is impossible to find an explicit form of the solution. Therefore, our goal is to obtain a numerical solution via the FE method.

#### *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: Ð *<sup>b</sup> <sup>a</sup> u*<sup>0</sup> *v*0 *dx* <sup>þ</sup> <sup>Ð</sup> *<sup>b</sup> <sup>a</sup> quvdx* <sup>¼</sup> <sup>Ð</sup> *<sup>b</sup> <sup>a</sup> fvdx*, for any test function *v* in the Sobolev space *H*<sup>1</sup> <sup>0</sup>½ �¼ *<sup>a</sup>*, *<sup>b</sup> <sup>v</sup>*<sup>∈</sup> *<sup>L</sup>*<sup>2</sup> ½ � *<sup>a</sup>*, *<sup>b</sup>* : �

k k*v* <sup>2</sup> <sup>þ</sup> *<sup>v</sup>*<sup>0</sup> k k<sup>2</sup> <sup>&</sup>lt; <sup>∞</sup>, *v a*ð Þ¼ *v b*ð Þ¼ <sup>0</sup> � , where k k*v* <sup>2</sup> <sup>¼</sup> <sup>Ð</sup> *<sup>b</sup> <sup>a</sup> <sup>v</sup>*<sup>2</sup>ð Þ *<sup>x</sup> 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 x*ð Þ<sup>∈</sup> *<sup>H</sup>*<sup>1</sup> <sup>0</sup>½ � *a*,*b* Ð *b a* 1 <sup>2</sup> *<sup>v</sup>*<sup>0</sup> ð Þ<sup>2</sup> <sup>þ</sup> <sup>1</sup> <sup>2</sup> *qv*<sup>2</sup> � *fv* � �*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*<sup>00</sup> *exists and continuous on a*½ � , *b* . *Then we have the following mathematical equivalences*.

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

*A Brief Summary of the Finite Element Method for Differential Equations DOI: http://dx.doi.org/10.5772/intechopen.95423*

#### *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*<sup>1</sup> <sup>0</sup> <sup>¼</sup> *<sup>v</sup>*<sup>∈</sup> *<sup>L</sup>*<sup>2</sup> ð Þ *a*, *b* : k k*v* <sup>2</sup> <sup>þ</sup> *<sup>v</sup>*<sup>0</sup> k k<sup>2</sup> <sup>&</sup>lt; <sup>∞</sup>, *v a*ð Þ¼ *v b*ð Þ¼ <sup>0</sup> n o. Multiplying (3) by a test function *v*∈ *H*<sup>1</sup> 0, integrating from *a* to *b*, and using integration by parts, we get

$$\int\_{a}^{b} fv d\mathbf{x} = \int\_{a}^{b} -u''v d\mathbf{x} + \int\_{a}^{b} quv d\mathbf{x} = \int\_{a}^{b} u'v' d\mathbf{x} + \int\_{a}^{b} quv d\mathbf{x},$$

since *v a*ð Þ¼ *v b*ð Þ¼ 0. Hence, the weak (or variational) form of (3) reads: Find *u*∈ *H*<sup>1</sup> 0, such that

$$\int\_{a}^{b} u'v' d\mathbf{x} + \int\_{a}^{b} quv d\mathbf{x} = \int\_{a}^{b} fv d\mathbf{x}, \quad \forall \ v \in H\_0^1. \tag{4}$$

We want to find *u*∈ *H*<sup>1</sup> <sup>0</sup> that satisfies (4). We note that a solution *u* to (4) is less regular than the solution *u* (3). Indeed, (4) has only *u*<sup>0</sup> whereas (3) contains *u*00. Furthermore, we can easily verify the following:


From now on, we use the notation k k*v* ¼ k k*v* <sup>Ω</sup>, where Ω ¼ ½ � *a*, *b* .

**The FE formulation**: The FE method is based on the variational form (4). We note that the space *H*<sup>1</sup> <sup>0</sup> contains many functions and it is therefore just as hard to find a function *u*∈ *H*<sup>1</sup> <sup>0</sup> 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*<sup>0</sup> <*x*<sup>1</sup> < ⋯ <*xN* ¼ *b* be a regular partition of ½ � *a*, *b* . Suppose that the length of *Ii* ¼ *xi*�1, *xi* ½ � is *hi* ¼ *xi* � *xi*�1. We define *h* ¼ max *i*¼1, 2, … , *N hi* to be the mesh size. We wish to construct a subspace *Vh* <sup>⊂</sup>*<sup>V</sup>* <sup>¼</sup> *<sup>H</sup>*<sup>1</sup> 0. Since *Vh* has finite dimension, we can find one set of basis functions *ϕ<sup>j</sup>* n o*<sup>N</sup>*�<sup>1</sup> *j*¼1 for *Vh*, where *ϕ<sup>j</sup>* ∈*Vh*, *j* ¼ 1, 2, … , *N* � 1 are linearly independent. We remark that *Vh* is the space spanned by the basis functions *<sup>i</sup>:e:*, *Vh* <sup>¼</sup> *vh*ð Þ *<sup>x</sup>* , *vh*ð Þ¼ *<sup>x</sup>* <sup>P</sup>*<sup>N</sup>*�<sup>1</sup> *<sup>j</sup>*¼<sup>1</sup> *cjϕj*ð Þ *<sup>x</sup>* n o. The FE method consists of choosing a basis for the subspace *Vh* that satisfies the following properties


It is natural to obtain an approximation *uh* to *u* as follows: Find *uh* ∈*Vh* such that

$$\int\_{a}^{b} u\_{h}' v' d\mathfrak{x} + \int\_{a}^{b} qu\_{h} v d\mathfrak{x} = \int\_{a}^{b} f v d\mathfrak{x}, \quad \forall \ v \in V\_{h}. \tag{5}$$

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

**FE approximation using Lagrange** <sup>1</sup> **elements**: The simplest finite dimensional space is the piecewise continuous linear function space defined over the triangulation

*V*1 *<sup>h</sup>*,0 ¼ f*vh* ∈*V*, *vh* is piecewise continuous linear over a, b ½ � with *vh*ð Þ¼ *a vh*ð Þ¼ *b* 0g*:*

It is easy to show that *V*<sup>1</sup> *<sup>h</sup>*,0 has a finite dimension even although there are infinite number of elements in *V*<sup>1</sup> *<sup>h</sup>*,0. The approximation of the FE method is therefore to look for an approximation *uh* within a small (finite dimensional) subspace *V*<sup>1</sup> *<sup>h</sup>*,0 ¼ *v*∈*V*<sup>1</sup> *<sup>h</sup>*<sup>j</sup> *v a*ð Þ¼ *v b*ð Þ¼ <sup>0</sup> � � of *<sup>H</sup>*<sup>1</sup> 0, consisting of piecewise linear polynomials, where *V*1 *<sup>h</sup>* <sup>¼</sup> *<sup>v</sup>*∈*C*<sup>0</sup>½ � *<sup>a</sup>*, *<sup>b</sup>* j j *<sup>v</sup> Ii* <sup>∈</sup>*P*<sup>1</sup> ð Þ *Ii* n o.

Let *V*<sup>1</sup> *<sup>h</sup>*,0 be the space of all continuous piecewise linear functions, which vanish at the end points *a* and *b*. There are many types of basis functions f g *ϕ<sup>i</sup> N*�1 *<sup>i</sup>*¼<sup>1</sup> . The simplest ones are the so-called hat functions satisfying *ϕ<sup>i</sup> x <sup>j</sup>* � � <sup>¼</sup> *<sup>δ</sup>ij*, where *<sup>δ</sup>ij* is the Kronecker symbol. Note especially that there is no need to construct hat functions *ϕ*<sup>0</sup> and *ϕ<sup>N</sup>* since any function of *V*<sup>1</sup> *<sup>h</sup>*,0 must vanish at the end points *x*<sup>0</sup> ¼ *a* and *xN* ¼ *b*.

The explicit expressions for the hat function *ϕi*ð Þ *x* and its derivative *ϕ*<sup>0</sup> *i* ð Þ *x* are given by

$$\phi\_i(\mathbf{x}) = \begin{cases} 0, & a \le \mathbf{x} \le \mathbf{x}\_{i-1}, \\ \frac{\mathbf{x} - \mathbf{x}\_{i-1}}{h\_i}, & \mathbf{x}\_{i-1} \le \mathbf{x} \le \mathbf{x}\_i, \\ \frac{\mathbf{x}\_{i+1} - \mathbf{x}}{h\_{i+1}}, & \mathbf{x}\_i \le \mathbf{x} \le \mathbf{x}\_{i+1}, \\ 0, & \mathbf{x}\_{i+1} \le \mathbf{x} \le b, \end{cases}, \quad \phi'\_i(\mathbf{x}) = \begin{cases} 0, & a < \mathbf{x} < \mathbf{x}\_{i-1}, \\ \frac{1}{h\_i}, & \mathbf{x}\_{i-1} < \mathbf{x} < \mathbf{x}\_i, \\ -\frac{1}{h\_{i+1}}, & \mathbf{x}\_i < \mathbf{x} < \mathbf{x}\_{i+1}, \\ 0, & \mathbf{x}\_{i+1} < \mathbf{x} < b, \end{cases},$$

for *<sup>i</sup>* <sup>¼</sup> 1, 2, … , *<sup>N</sup>* � 1. The FE approximation of (4) thus reads: Find *<sup>u</sup>* <sup>∈</sup>*V*<sup>1</sup> *<sup>h</sup>*,0, such that

$$\int\_{a}^{b} u'\_{h} v' d\mathfrak{x} + \int\_{a}^{b} qu\_{h} v d\mathfrak{x} = \int\_{a}^{b} f v d\mathfrak{x}, \quad \forall \ v \in V^{1}\_{h,0}. \tag{6}$$

We call *uh* the FE approximation of *u*. We say that (6) is the Galerkin approximation of (4) and the method used to find *uh* ∈*V*<sup>1</sup> *<sup>h</sup>*,0 is called Galerkin method.

It can be shown that (6) is equivalent to the *N* � 1 equations

*A Brief Summary of the Finite Element Method for Differential Equations DOI: http://dx.doi.org/10.5772/intechopen.95423*

$$\int\_{a}^{b} u\_{h}^{\prime} \phi\_{i}^{\prime} d\mathfrak{x} + \int\_{a}^{b} qu\_{h} \phi\_{i} d\mathfrak{x} = \int\_{a}^{b} f \phi\_{i} d\mathfrak{x}, \quad i = 1, 2, \dots, N - 1. \tag{7}$$

**Derivation of the discrete system**: Since *uh* ∈*V*<sup>1</sup> *<sup>h</sup>*,0, we can express it as a linear combination of hat functions *i:e:*,

$$u\_h = \sum\_{j=1}^{N-1} c\_j \phi\_j(\mathbf{x}),\tag{8}$$

where *cj* are real numbers to be determined. We note that the coefficients *cj*, *j* ¼ 1, 2, … , *N* � 1 are the *N* � 1 nodal values of *uh* to be determined. Note that the index is only from 1 to *N* � 1, because of the zero boundary conditions. We remark that *uh*ð Þ¼ *a uh*ð Þ¼ *b* 0 and *uh*ð Þ¼ *xi ci*. So *ci* is an approximate solution to the exact solution at *x* ¼ *xi*.

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

Substituting (8) into (7) yields

$$\sum\_{j=1}^{N-1} c\_j \left( \int\_a^b \phi\_i' \phi\_j' d\mathbf{x} + \int\_a^b q \phi\_i \phi\_j d\mathbf{x} \right) = \int\_a^b f \phi\_i d\mathbf{x}, \quad i = 1, 2, \dots, N - 1. \tag{9}$$

The problem (7) is now equivalent to the following: Find the real numbers *c*1,*c*2, … ,*cN*�<sup>1</sup> 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},\tag{10}$$

where **<sup>c</sup>** <sup>¼</sup> ½ � *<sup>c</sup>*1,*c*2, … ,*cN*�<sup>1</sup> *<sup>t</sup>* <sup>∈</sup> *<sup>N</sup>*�<sup>1</sup> is the unknown vector, *<sup>A</sup>* is an ð Þ� *<sup>N</sup>* � <sup>1</sup> ð Þ *N* � 1 matrix, the so-called stiffness matrix when *q* ¼ 0, with entries

$$a\_{i\bar{j}} = \int\_{a}^{b} (\phi\_i' \phi\_j' + q \phi\_i \phi\_j) dx, \quad i, j = 1, 2, \dots, N - 1,\tag{11}$$

and **b**∈ *<sup>N</sup>*�<sup>1</sup> , the so-called load vector, has entries

$$b\_i = \int\_a^b f\phi\_i d\mathbf{x}, \quad i = \mathbf{1}, 2, \dots, N - \mathbf{1}. \tag{12}$$

To obtain the approximate solution we need to solve the linear system for the unknown vector **c**. We note that *aij* ¼ *a ϕi*, *ϕ<sup>j</sup>* � � and *bi* <sup>¼</sup> *<sup>f</sup>*, *<sup>ϕ</sup><sup>i</sup>* ð Þ, where *a u*ð Þ¼ , *<sup>v</sup>* Ð *b <sup>a</sup> u*<sup>0</sup> *<sup>v</sup>* ð Þ <sup>0</sup> <sup>þ</sup> *quv dx* is a bi-linear and ð Þ¼ *<sup>f</sup>*, *<sup>v</sup>* <sup>Ð</sup> *<sup>b</sup> <sup>a</sup> 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

$$\min\_{\nu(\mathbf{x}) \in H\_0^1[a,b]} F(\nu), \quad \text{where} \ F(\nu) = \int\_a^b \left(\frac{1}{2}(\nu')^2 + \frac{1}{2}q\nu^2 - f\nu\right) d\mathbf{x}.$$

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

$$F(u\_h) = \int\_{a}^{b} \left(\frac{1}{2} \left(\sum\_{j=1}^{N-1} c\_j \phi\_j'(\mathbf{x})\right)^2 + \frac{1}{2} q \left(\sum\_{j=1}^{N-1} c\_j \phi\_j(\mathbf{x})\right)^2 - f\left(\sum\_{j=1}^{N-1} c\_j \phi\_j(\mathbf{x})\right)\right) d\mathbf{x}, \quad u\_h$$

which is a multi-variable function of *c*1,*c*2, … ,*cN*�<sup>1</sup> and can be written as *F u*ð Þ¼ *<sup>h</sup> F c*ð Þ 1,*c*2, … ,*cN*�<sup>1</sup> . The necessary conditions for a global minimum are *<sup>∂</sup><sup>F</sup> <sup>∂</sup>ci* ¼ 0, *i* ¼ 1, 2, … , *N* � 1. Taking the partial derivatives directly with respect to *ci*, we get

$$\int\_{a}^{b} \left(\phi\_i'(\mathbf{x}) \sum\_{j=1}^{N-1} c\_j \phi\_j'(\mathbf{x}) + q \phi\_i(\mathbf{x}) \sum\_{j=1}^{N-1} c\_j \phi\_j(\mathbf{x}) - f \phi\_i(\mathbf{x})\right) d\mathbf{x} = \mathbf{0}, \quad 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(\phi\_i'(\mathbf{x})\phi\_j'(\mathbf{x}) + q\phi\_i(\mathbf{x})\phi\_j(\mathbf{x})\right) d\mathbf{x} = \int\_a^b f\phi\_i(\mathbf{x})d\mathbf{x} = \mathbf{0}, \quad i = 1, 2, \dots, N-1, \quad j = 1, 2, \dots, N$$

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 *<sup>a</sup>*^*<sup>i</sup>*,*<sup>j</sup>* <sup>¼</sup> <sup>Ð</sup> *<sup>b</sup> <sup>a</sup> ϕ*<sup>0</sup> *i ϕ*0 *j dx*. For ∣*i* � *j*∣>1, we have *a*^*<sup>i</sup>*,*<sup>j</sup>* ¼ 0, since *ϕ<sup>i</sup>* and *ϕ<sup>j</sup>* lack overlapping support. However, if *i* ¼ *j*, then

$$\hat{a}\_{i,i} = \int\_{a}^{b} \left(\phi\_{i}^{\prime}\right)^{2} d\mathbf{x} = \int\_{x\_{i-1}}^{x\_{i}} \left(\frac{\mathbf{1}}{h\_{i}}\right)^{2} d\mathbf{x} + \int\_{x\_{i}}^{x\_{i+1}} \left(-\frac{1}{h\_{i+1}}\right)^{2} d\mathbf{x} = \frac{1}{h\_{i}} + \frac{1}{h\_{i+1}}, \quad i, j = 1, 2, \dots, N - 1.$$

Furthermore, if *j* ¼ *i* þ 1, then

$$\hat{a}\_{i,i+1} = \int\_{a}^{b} \phi\_{i}^{\prime} \phi\_{i+1}^{\prime} dx = \int\_{x\_{i}}^{x\_{i+1}} \left( -\frac{1}{h\_{i+1}} \right) \left( \frac{1}{h\_{i+1}} \right) dx = -\frac{1}{h\_{i+1}}, \quad i, j = 1, 2, \dots, N-2. \tag{13}$$

By symmetry, we also have

$$\hat{a}\_{i+1,i} = \int\_a^b \phi'\_{i+1} \phi'\_i d\infty = -\frac{1}{h\_{i+1}}, \quad i, j = 1, 2, \dots, N - 2.$$

To obtain *<sup>a</sup>*~*<sup>i</sup>*,*<sup>j</sup>* <sup>¼</sup> <sup>Ð</sup> *<sup>b</sup> <sup>a</sup> <sup>q</sup>ϕiϕjdx* and *bi* <sup>¼</sup> <sup>Ð</sup> *<sup>b</sup> <sup>a</sup> fϕidx*, we use the composite trapezoidal rule

$$\int\_{a}^{b} f(\mathbf{x})d\mathbf{x} = \sum\_{i=1}^{N} \int\_{\mathbf{x}\_{i-1}}^{\mathbf{x}\_{i}} f(\mathbf{x})d\mathbf{x} \approx \frac{1}{2} \left[ h\_{\mathbf{f}}f(\mathbf{x}\_{0}) + \sum\_{i=1}^{N-1} (h\_{i} + h\_{i+1})f(\mathbf{x}\_{i}) + h\_{\mathbf{N}}f(\mathbf{x}\_{N}) \right].$$

So, we can easily verify that

$$\tilde{a}\_{i\bar{j}} = \int\_{a}^{b} q \phi\_i \phi\_j d\mathbf{x} \approx \begin{cases} \frac{q\_i}{2} (h\_i + h\_{i+1}), & i = j \\ 0, & i \neq j \end{cases}, \quad b\_i = \int\_{a}^{b} f \phi\_i d\mathbf{x} \approx \frac{1}{2} (h\_i + h\_{i+1}) f\_i,$$

*A Brief Summary of the Finite Element Method for Differential Equations DOI: http://dx.doi.org/10.5772/intechopen.95423*

where *qi* ¼ *q x*ð Þ*<sup>i</sup>* and *fi* ¼ *f x*ð Þ*<sup>i</sup>* . Thus, the matrix *A* ¼ *a*^*i*,*<sup>j</sup>* þ *a*~*i*,*<sup>j</sup>* � � is tridiagonal and has the form

$$A = \begin{bmatrix} \frac{1}{h\_1} + \frac{1}{h\_2} + \frac{q\_1}{2}(h\_1 + h\_2) & -\frac{1}{h\_2} & 0 & \cdots & 0\\ -\frac{1}{h\_2} & \frac{1}{h\_2} + \frac{1}{h\_3} + \frac{q\_2}{2}(h\_2 + h\_3) & -\frac{1}{h\_3} & \ddots & & 0\\ 0 & -\frac{1}{h\_3} & \ddots & \cdots & 0\\ \vdots & \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}(h\_{N-1} + h\_N) \end{bmatrix}.$$

Finally, we obtain the following system: *c*<sup>0</sup> ¼ *cN* ¼ 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(h\_i + h\_{i+1})}{2}c\_i = \frac{1}{2}(h\_i + h\_{i+1})f\_i, \quad i = 1, 2, \dots, N - 1, 2$$

**Remark 2.2** *Suppose that the partition is uniform i:e:*, *hi* <sup>¼</sup> *<sup>h</sup>* <sup>¼</sup> *<sup>b</sup>*�*<sup>a</sup> <sup>N</sup> for all i* ¼ 1, 2, … , *N. Then the stiffness matrix A and the load vector* **b** *have the form:*

$$A = \begin{bmatrix} \frac{2}{h} + hq\_1 & -\frac{1}{h} & 0 & \cdots & 0\\ -\frac{1}{h} & \frac{2}{h} + hq\_2 & -\frac{1}{h} & \ddots & & 0\\ 0 & -\frac{1}{h} & \ddots & \cdots & 0\\ \vdots & \ddots & \ddots & \ddots & -\frac{1}{h\_{N-1}}\\ 0 & \cdots & 0 & -\frac{1}{h} & \frac{2}{h} + hq\_{N-1} \end{bmatrix}, \quad \mathbf{b} = h \begin{bmatrix} f\_1\\ f\_2\\ f\_3\\ \vdots\\ f\_{N-1} \end{bmatrix}.$$

Finally, we obtain the following system: *c*<sup>0</sup> ¼ *cN* ¼ 0 and

$$\frac{-c\_{j-1} + 2c\_j - c\_{j+1}}{h} + hq\_i c\_j = hf\_i \Rightarrow -\frac{c\_{j-1} - 2c\_j + c\_{j+1}}{h^2} + q\_i c\_j = f\_i, \quad i = 1, 2, \dots, N - 1, \dots$$

which is the same system obtained using the finite difference method, where *u*<sup>00</sup> is approximated using the second-order midpoint formula *u*<sup>00</sup> *x <sup>j</sup>* � �≈ *u x*ð Þ *<sup>j</sup>*�<sup>1</sup> �2*u x*ð Þ*<sup>j</sup>* <sup>þ</sup>*u x*ð Þ *<sup>j</sup>*þ<sup>1</sup> *<sup>h</sup>*<sup>2</sup> . 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.1** *The matrix A with entries ai*,*<sup>j</sup>* <sup>¼</sup> <sup>Ð</sup> *<sup>b</sup> <sup>a</sup> ϕ*<sup>0</sup> *i ϕ*0 *j dx is symmetric positive definite i:e:*, *ai*,*<sup>j</sup>* ¼ *aj*,*<sup>i</sup> and*

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

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

Next, we state a general convergence result for the Galerkin method. We first define the following norm and semi-norm: For *v*∈ *H*<sup>1</sup> 0, we define

$$||v|| = \left(\int\_{a}^{b} v^2(\mathfrak{x})d\mathfrak{x}\right)^{1/2}, \quad |v|\_1 = ||v'|| = \left(\int\_{a}^{b} (v'(\mathfrak{x}))^2 d\mathfrak{x}\right)^{1/2}.$$

**Theorem 2.3** *Suppose that q x*ð Þ≥ 0, ∀ *x*∈ ½ � *a*, *b* . *Let u be the solution to (4) and uh be the solution to (6). Then there exists a constant C such that*

$$\left| \left| (u - u\_h)' \right| \right| \le C \left| \left| (u - v\_h)' \right| \right|, \quad \forall \ v\_h \in V^1\_{h,0}, \tag{14}$$

where *C* is given by *C* ¼ 1 þ max *<sup>x</sup>*∈½ � *<sup>a</sup>*,*<sup>b</sup>* ∣*q x*ð Þ∣, which is independent of the choice of *V*<sup>1</sup> *<sup>h</sup>*,0.

**Remark 2.3** *From (14), taking the minimum over all vh* ∈*V*<sup>1</sup> *<sup>h</sup>*,0*, we get* ð Þ *<sup>u</sup>* � *uh* <sup>0</sup> � � � � ≤*C* min *vh* ∈*Vh*,0 ð Þ *<sup>u</sup>* � *vh* <sup>0</sup> � � � �*. Thus, u*j j � *uh* <sup>1</sup> ≤*C* min *vh* <sup>∈</sup>*Vh*,0 j j *u* � *vh* <sup>1</sup>*, where C* ¼ 1 þ max *<sup>x</sup>*∈½ � *<sup>a</sup>*,*<sup>b</sup>* ∣*q x*ð Þ∣*.*

Next, we study the convergence of *uh* to *u*. Let *u* ∈ *H*<sup>1</sup> 0. Define the piecewise linear interpolant by

$$\pi u = \sum\_{j=1}^{N} u\left(\mathbf{x}\_{j}\right) \phi\_{j}(\mathbf{x}) \in V^{1}\_{h,0}, \quad \mathbf{x} \in [a, b].$$

Since *πu*∈*V*<sup>1</sup> *<sup>h</sup>*,0, the estimate (14) gives

$$\left| \left| (\boldsymbol{u} - \boldsymbol{u}\_h)' \right| \right| \le C \left| \left| (\boldsymbol{u} - \boldsymbol{\pi}\boldsymbol{u})' \right| \right|.$$

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

**Theorem 2.4** (**A priori error estimate**) *Suppose that q x*ð Þ≥0 ∀ *x*∈½ � *a*, *b* . *Let u be the solution to (4) and uh be the solution to (6). Then there exists a constant C such that*

$$\left\| \left( (\boldsymbol{\mu} - \boldsymbol{u}\_h)' \right\| ^2 \leq \mathcal{C} \sum\_{i=1}^N h\_i^2 \| \boldsymbol{u}'' \| \_{I\_i}^2,$$

where *C* is a constant independent of *h*. Consequently, if *h* ¼ max *ihi*, then

$$\left\| \left( (\boldsymbol{\mu} - \boldsymbol{\mu}\_h)' \right) \right\|^2 \le Ch^2 \left\| \left\| \boldsymbol{\mu}'' \right\| ^2.$$

#### **Remark 2.4**


### *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):


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

**Nonhomogeneous Dirichlet boundary conditions**: Let us consider the following two-point BVP: find *u*∈*C*<sup>2</sup> ð Þ *a*, *b* such that

$$-u'' = f(x), \quad x \in (a, b), \quad u(a) = a, \quad u(b) = \beta,\tag{15}$$

where *α* and *β* are given constants and *f* ∈*C a*ð Þ , *b* is a given function. In this case, the admissible function space *H*<sup>1</sup> <sup>0</sup> ¼ *v* : k k*v* <sup>2</sup> <sup>þ</sup> *<sup>v</sup>*<sup>0</sup> k k<sup>2</sup> <sup>&</sup>lt; <sup>∞</sup>, *v a*ð Þ¼ *v b*ð Þ¼ <sup>0</sup> n o and the FE space *V*<sup>1</sup> *<sup>h</sup>*,0 defined earlier remain the same. Multiplying (15) by a test function *v*∈ *H*<sup>1</sup> <sup>0</sup> and integrating by parts gives

$$\int\_{a}^{b} fv d\mathbf{x} = \int\_{a}^{b} -u''v d\mathbf{x} = -u'(b)v(b) + u'(a)v(a) + \int\_{a}^{b} u'v' d\mathbf{x} = \int\_{a}^{b} u'v' d\mathbf{x},$$

since *v a*ð Þ¼ *v b*ð Þ¼ 0. Hence, the weak or variational form of (15) reads: Given *u a*ð Þ¼ *<sup>α</sup>*, *u b*ð Þ¼ *<sup>β</sup>*, find *<sup>u</sup>*<sup>∈</sup> *<sup>H</sup>*<sup>1</sup> <sup>¼</sup> *<sup>v</sup>* : k k*<sup>v</sup>* <sup>2</sup> <sup>þ</sup> *<sup>v</sup>*<sup>0</sup> k k<sup>2</sup> <sup>&</sup>lt; <sup>∞</sup> n o, such that

$$\int\_{a}^{b} u'v' d\mathbf{x} = \int\_{a}^{b} fv d\mathbf{x}, \quad \forall \ v \in H\_0^1. \tag{16}$$

Let *V*<sup>1</sup> *<sup>h</sup>* and *V*<sup>1</sup> *<sup>h</sup>*,0, respectively, be the space of all continuous piecewise linear functions and the space of all continuous piecewise linear functions which vanish at the endpoints *a* and *b*. We also let *a* ¼ *x*<sup>0</sup> <*x*<sup>1</sup> < ⋯ < *xN* ¼ *b* be a uniform partition of the interval ½ � *a*, *b* . Moreover let f g *ϕ<sup>i</sup>* be the set of hat basis functions of *Vh* associated with the *N* þ 1 nodes *x <sup>j</sup>*, *j* ¼ 0, 1, … , *N*, such that *ϕ<sup>i</sup> x <sup>j</sup>* � � <sup>¼</sup> *<sup>δ</sup>ij*. The FE approximation of (16) thus reads: Find *uh* <sup>∈</sup>*V*<sup>1</sup> *<sup>h</sup>* such that *uh*ð Þ¼ *a α*, *uh*ð Þ¼ *b β*, and

$$\int\_{a}^{b} u'\_{h} v' d\mathfrak{x} = \int\_{a}^{b} f v d\mathfrak{x}, \quad \forall \ v \in V\_{h,0}^{1}. \tag{17}$$

It can be shown that (17) is equivalent to the *N* � 1 equations

$$\int\_{a}^{b} u'\_{h} \phi'\_{i} d\mathbf{x} = \int\_{a}^{b} f \phi\_{i} d\mathbf{x}, \quad i = 1, 2, \dots, N - 1. \tag{18}$$

Expanding *uh* as a linear combination of hat functions

$$u\_h = \sum\_{j=0}^{N} c\_j \phi\_j(\mathbf{x}) = a\phi\_0(\mathbf{x}) + \sum\_{j=1}^{N-1} c\_j \phi\_j(\mathbf{x}) + \beta \phi\_N(\mathbf{x}),\tag{19}$$

where the coefficients *cj*, *j* ¼ 1, 2, … , *N* � 1 are the *N* � 1 nodal values of *uh* to be determined.

Substituting (19) into (18) yields

$$\sum\_{j=1}^{N-1} c\_j \left( \int\_a^b \phi\_i' \phi\_j' d\mathbf{x} \right) = \int\_a^b \left( f \phi\_i - a \phi\_0' \phi\_i' - \beta \phi\_N' \phi\_i' \right) d\mathbf{x}, \quad i = 1, 2, \dots, N - 1, 2$$

which is a ð Þ� *N* � 1 ð Þ *N* � 1 system of equations for *cj*. In matrix form we write

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

where *A* is a ð Þ� *N* � 1 ð Þ *N* � 1 matrix, the so-called stiffness matrix, with entries

$$a\_{i,j} = \int\_{a}^{b} \phi\_i' \phi\_j' d\mathbf{x}, \quad i, j = 1, 2, \dots, N - 1,\tag{21}$$

**<sup>c</sup>** <sup>¼</sup> ½ � *<sup>c</sup>*1,*c*2, … ,*cN*�<sup>1</sup> *<sup>t</sup>* is a ð Þ *<sup>N</sup>* � <sup>1</sup> vector containing the unknown coefficients *cj*, *j* ¼ 1, 2, … , *N* � 1, and **b** is a ð Þ *N* � 1 vector, the so-called load vector, with entries

$$b\_i = \int\_a^b (f\phi\_i - a\phi\_0'\phi\_i' - \beta\phi\_N'\phi\_i')d\mathbf{x}, \quad i = 1, 2, \dots, N - 1. \tag{22}$$

**Computer Implementation**: The explicit expression for a hat function *ϕi*ð Þ *x* is given by

$$\begin{aligned} \phi\_i(\mathbf{x}) &= \begin{cases} 0, & a \le \mathbf{x} \le \mathbf{x}\_{i-1}, \\ \frac{\mathbf{x} - \mathbf{x}\_{i-1}}{h\_i}, & \mathbf{x}\_{i-1} < \mathbf{x} \le \mathbf{x}\_i, \\ \frac{\mathbf{x}\_{i+1} - \mathbf{x}}{h\_{i+1}}, & \mathbf{x}\_i < \mathbf{x} \le \mathbf{x}\_{i+1}, \end{cases}, \quad i = 1, 2, \dots, N - 1, \\\ 0, & \quad \mathbf{x}\_{i+1} < \mathbf{x} \le b, \\\ \phi\_0(\mathbf{x}) &= \begin{cases} \frac{\mathbf{x}\_1 - \mathbf{x}}{h\_1}, & \mathbf{x}\_0 < \mathbf{x} \le \mathbf{x}\_1, \\ 0, & \mathbf{x}\_1 < \mathbf{x} \le b, \end{cases}, \quad \phi\_N(\mathbf{x}) = \begin{cases} 0, & \mathbf{x}\_0 < \mathbf{x} \le \mathbf{x}\_{N-1}, \\ \frac{\mathbf{x} - \mathbf{x}\_{N-1}}{h\_N}, & \mathbf{x}\_{N-1} < \mathbf{x} \le b. \end{cases} \end{aligned}$$

For simplicity we assume the partition is uniform so that *hi* ¼ *h* for *i* ¼ 1, 2, … , *N*. Hence the derivative *ϕ*<sup>0</sup> *i* ð Þ *<sup>x</sup>* is either � <sup>1</sup> *h*, 1 *<sup>h</sup>*, or 0 depending on the interval.

It is straightforward to calculate the entries of the stiffness matrix. For ∣*i* � *j*∣> 1, we have *ai*,*<sup>j</sup>* ¼ 0, since *ϕ<sup>i</sup>* and *ϕ<sup>j</sup>* lack overlapping support. However, if *i* ¼ *j*, then

*A Brief Summary of the Finite Element Method for Differential Equations DOI: http://dx.doi.org/10.5772/intechopen.95423*

$$\mathfrak{a}\_{i,j} = \int\_{a}^{b} \left(\phi\_{i}^{\prime}\right)^{2} d\mathfrak{x} = \int\_{\chi\_{i-1}}^{\chi\_{i}} \left(\frac{1}{h}\right)^{2} d\mathfrak{x} + \int\_{\chi\_{i}}^{\chi\_{i+1}} \left(-\frac{1}{h}\right)^{2} d\mathfrak{x} = \frac{2}{h}, \quad i, j = 1, 2, \dots, N - 1, 2$$

where we have used that *xi* � *xi*�<sup>1</sup> ¼ *xi*þ<sup>1</sup> � *xi* ¼ *h*. Furthermore, if *j* ¼ *i* þ 1, then

$$a\_{i,i+1} = \int\_{a}^{b} \phi\_i' \phi\_{i+1}' d\mathbf{x} = \int\_{x\_i}^{x\_{i+1}} \left(-\frac{1}{h}\right) \left(\frac{1}{h}\right) d\mathbf{x} = -\frac{1}{h}, \quad i, j = 1, 2, \dots, N - 2.$$

Changing *i* to *i* � 1 we also have

$$\mathfrak{a}\_{i-1,i} = \int\_{a}^{b} \phi\_{i-1}^{\prime} \phi\_{i} d\mathbf{x} = \int\_{\mathbf{x}\_{i-1}}^{\mathbf{x}\_{i}} \left(\frac{1}{h}\right) \left(-\frac{1}{h}\right) d\mathbf{x} = -\frac{1}{h}, \quad i, j = 2, 3, \dots, N - 1.$$

Thus the stiffness matrix is

$$A = \frac{1}{h} \begin{bmatrix} 2 & -1 & 0 & \cdots & 0 \\ -1 & 2 & -1 & \ddots & 0 \\ 0 & -1 & 2 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & -1 \\ 0 & \cdots & 0 & -1 & 2 \end{bmatrix}.$$

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

$$\begin{split} b\_{1} &= \int\_{a}^{b} \left( f\phi\_{1} - a\phi\_{0}^{t}\phi\_{1}^{t} - \beta\phi\_{N}^{t}\phi\_{1}^{t} \right) dx = \int\_{x\_{0}}^{x\_{1}} \left( f\phi\_{1} - a\left( -\frac{1}{h} \right) \left( \frac{1}{h} \right) \right) dx + \int\_{x\_{1}}^{x\_{1}} f\phi\_{1} \\ &= \frac{a}{h} + \int\_{x\_{0}}^{x\_{1}} f\phi\_{1} \approx \frac{a}{h} + hf(\mathbf{x}\_{1}), \end{split}$$

$$\begin{split} b\_{i} &= \int\_{a}^{b} \left( f\phi\_{i} - a\phi\_{0}^{t}\phi\_{i}^{t} - \beta\phi\_{N}^{t}\phi\_{i}^{t} \right) dx = \int\_{x\_{i-1}}^{x\_{i+1}} f\phi\_{i} dx \approx hf(\mathbf{x}\_{i}), \quad i = 2, \dots, N-2, \\ b\_{N-1} &= \int\_{a}^{b} \left( f\phi\_{N-1}^{t} - a\phi\_{0}^{t}\phi\_{N-1} - \beta\phi\_{N}^{t}\phi\_{N-1}^{t} \right) dx \\ &= \int\_{x\_{N-1}}^{x\_{N-1}} f\phi\_{N-1} dx + \int\_{x\_{N-1}}^{x\_{N}} \left( f\phi\_{N-1} - \beta\left( \frac{1}{h} \right) \left( -\frac{1}{h} \right) \right) dx = \int\_{x\_{N-1}}^{x\_{N}} f\phi\_{N-1} dx + \frac{\rho}{h} \approx hf(\mathbf{x}\_{N-1}) + \frac{\rho}{h}, \end{split}$$

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

$$\frac{1}{h} \begin{bmatrix} 2 & -1 & 0 & \cdots & 0 \\ -1 & 2 & -1 & \ddots & 0 \\ 0 & -1 & 2 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & -1 \\ 0 & \cdots & 0 & -1 & 2 \end{bmatrix} \begin{bmatrix} c\_1 \\ c\_2 \\ c\_3 \\ \vdots \\ c\_{N-1} \end{bmatrix} = \begin{bmatrix} hf(\mathbf{x}\_1) + \frac{a}{h} \\ hf(\mathbf{x}\_2) \\ \vdots \\ hf(\mathbf{x}\_N) \\ \vdots \\ hf(\mathbf{x}\_{N-1}) + \frac{\theta}{h} \end{bmatrix}$$

We note that *uh*ð Þ¼ *a α* ¼ *u a*ð Þ and *uh*ð Þ¼ *b β* ¼ *u b*ð Þ. Therefore, we see that the system matrix *A* remains the same, and only the first and last entries of the load vector **b** need to be modified because of the definition of the basis functions *ϕ*<sup>0</sup> f g , … , *ϕ<sup>N</sup>* . An alternative approach is to use all the basis functions *ϕ*<sup>0</sup> f g , … , *ϕ<sup>N</sup>* to form a larger system of equation, *i:e:*, and ð Þ� *N* þ 1 ð Þ *N* þ 1 system. The procedure for inserting the boundary conditions into the system equation is: enter zeros in the first and ð Þ *N* þ 1 -th rows of the system matrix *A* except for unity in the main diagonal positions of these two rows, and enter *α* and *β* in the first and ð Þ *N* þ 1 -th rows of the vector **b**, respectively.

**General boundary conditions**: Let us consider the following two-point BVP: find *u*∈*C*<sup>2</sup> ð Þ *a*, *b* such that

$$-u'' = f(\mathbf{x}), \quad \mathbf{x} \in [a, b], \quad u(a) = a, \quad \gamma u(b) + u'(b) = \beta,\tag{23}$$

where *α*, *β* and *γ* are given numbers and *f* ∈*C a*ð Þ , *b* is a given function. The boundary condition at *x* ¼ *b* is called a Robin boundary condition (combination and *u* and *u*<sup>0</sup> is prescribed at *x* ¼ *b*). In this case, the admissible function space is modified to

$$H\_0^1 = \left\{ v : \left\| v \right\|^2 + \left\| v' \right\|^2 < \infty, \ v(a) = \mathbf{0} \right\}.$$

Multiplying (23) by a function *v*∈ *H*<sup>1</sup> <sup>0</sup> and integrating by parts gives

$$\int\_{a}^{b} fv d\mathbf{x} = \int\_{a}^{b} -u''v d\mathbf{x} = -u'(b)v(b) + u'(a)v(a) + \int\_{a}^{b} u'v' d\mathbf{x}$$

$$= -(\beta - \gamma u(b))v(b) + u'(a)v(a) + \int\_{a}^{b} u'v' d\mathbf{x}.$$

Since *v a*ð Þ¼ 0, we are left with

$$\int\_{a}^{b} u'v' d\mathfrak{x} + \mathfrak{y}u(b)v(b) = \int\_{a}^{b} f\nu d\mathfrak{x} + \beta\nu(b).$$

Hence, the weak or variational form of (23) reads: Given *u a*ð Þ¼ *α*, find the approximate solution *u*∈ *H*<sup>1</sup> 0, such that

$$\int\_{a}^{b} u'v' d\mathbf{x} + \chi u(b)v(b) = \int\_{a}^{b} fv d\mathbf{x} + \beta v(b), \quad \forall \ v \in H\_0^1. \tag{24}$$

The FE space *V*<sup>1</sup> *<sup>h</sup>* is now the set of all continuous piecewise linear functions which vanish at the end point *a*. The FE approximation of (24) thus reads: Find the piecewise linear approximation *uh* to the solution *u* satisfies

$$\int\_{a}^{b} u'\_h v' d\mathbf{x} + \gamma u\_h(b) v(b) = \int\_{a}^{b} f v d\mathbf{x} + \beta v(b), \quad \forall \ v \in V\_h^1,\tag{25}$$

with *uh*ð Þ¼ *a α*. 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*∈*C*<sup>2</sup> ð Þ *a*, *b* such that *A Brief Summary of the Finite Element Method for Differential Equations DOI: http://dx.doi.org/10.5772/intechopen.95423*

$$\begin{aligned} -(p(\mathbf{x})u')' &= f(\mathbf{x}), & \mathbf{x} \in I &= [a, b], \quad p(a)u'(a) = \kappa\_0(u(a) - a), \\ p(b)u'(b) &= \kappa\_1(u(b) - \beta), \end{aligned} \tag{26}$$

where *<sup>p</sup>* <sup>¼</sup> *p x*ð Þ with *p x*ð Þ≥*p*<sup>0</sup> <sup>&</sup>gt;0, *<sup>f</sup>* <sup>∈</sup>*L*<sup>2</sup> ð Þ*I* , *κ*0, *κ*<sup>1</sup> ≥ 0, and *α*, *β* are given numbers. Let

$$V = \left\{ \boldsymbol{\nu} \in \mathcal{C}^0(I) : \left\| \boldsymbol{\nu} \right\|^2 + \left\| \boldsymbol{\nu}' \right\|^2 < \infty \right\}.$$

Multiplying (26) by a function *v*∈*V* and integrating by parts gives

$$\int\_{a}^{b} f v d\mathbf{x} = \int\_{a}^{b} -(pu')' v d\mathbf{x} = \int\_{a}^{b} pu' v' d\mathbf{x} - p(b)u'(b)v(b) + p(a)u'(a)v(a)$$

$$= \int\_{a}^{b} pu' v' d\mathbf{x} - \kappa\_1 (u(b) - \beta) v(b) + \kappa\_0 (u(a) - a)v(a).$$

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

$$\int\_{a}^{b} p u' v' d\mathbf{x} - \kappa\_1 u(b) v(b) + \kappa\_0 u(a) v(a) = \int\_{a}^{b} f v d\mathbf{x} - \kappa\_1 \theta v(b) + \kappa\_0 a v(a), \quad \forall \ v \in V.$$

The FE method consists of finding *uh* <sup>∈</sup>*Vh* <sup>¼</sup> *<sup>v</sup>*∈*C*<sup>0</sup>ð Þ *<sup>a</sup>*, *<sup>b</sup>* j j *<sup>v</sup> Ii* <sup>∈</sup>*P*<sup>1</sup> ð Þ *Ii* n o such that

$$\int\_{a}^{b} p u\_h' v' dx - \kappa \mathbf{u} u\_h(b) v(b) + \kappa u u\_h(a) v(a) = \int\_{a}^{b} f v dx - \kappa \boldsymbol{\chi} \boldsymbol{\beta} v(b) + \kappa \boldsymbol{\alpha} w(a), \quad \forall v \in V\_h. \tag{27}$$

**Implementation**: We need to assemble a stiffness matrix *A* and a load vector *b*. Substituting *uh* <sup>¼</sup> <sup>P</sup>*<sup>N</sup> <sup>i</sup>*¼<sup>0</sup>*ciϕ<sup>i</sup>* into (27) and taking *<sup>v</sup>* <sup>¼</sup> *<sup>ϕ</sup><sup>j</sup>* for *<sup>j</sup>* <sup>¼</sup> 0, 1, … , *<sup>N</sup>* yields

$$\sum\_{i=0}^{N} \int\_{a}^{b} p\phi\_{i}^{\prime} \phi\_{j}^{\prime} d\mathbf{x} - \kappa\_{1}\phi\_{i}(b)\phi\_{j}(b) + \kappa\_{0}\phi\_{i}(a)\phi\_{j}(a) = \int\_{a}^{b} f\phi\_{j} d\mathbf{x} - \kappa\_{1}\beta\phi\_{j}(b) + \kappa\_{0}a\phi\_{j}(a),$$
  $\forall j = 0, 1, \dots, N.$ 

which is a ð Þ� *N* þ 1 ð Þ *N* þ 1 system of equations for *ci*. In matrix form we write *<sup>A</sup>***<sup>c</sup>** <sup>¼</sup> **<sup>b</sup>**, where **<sup>c</sup>** <sup>¼</sup> ½ � *<sup>c</sup>*0, … ,*cN <sup>t</sup>* is a ð Þ *<sup>N</sup>* <sup>þ</sup> <sup>1</sup> vector containing the unknown coefficients *ci*, *i* ¼ 0, 1, … , *N*, *A* is a ð Þ� *N* þ 1 ð Þ *N* þ 1 matrix with entries

$$a\_{i,j} = \int\_a^b p\phi\_i'\phi\_j' d\mathbf{x} - \kappa\_1\phi\_i(b)\phi\_j(b) + \kappa\_0\phi\_i(a)\phi\_j(a), \quad i, j = 0, 1, \dots, N, \varepsilon$$

and **b** is a ð Þ *N* þ 1 vector with entries

$$b\_j = \int\_a^b f\phi\_j d\mathbf{x} - \kappa\_1 \beta \phi\_j(b) + \kappa\_0 a \phi\_j(a), \quad j = 0, 1, \dots, N.$$

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

$$A = \begin{bmatrix} \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\\ \vdots & \ddots & \ddots & \ddots & -\frac{1}{h\_N}\\ & & & &\\ 0 & \cdots & 0 & -\frac{1}{h\_N} & \frac{1}{h\_N} - \kappa\_1 \end{bmatrix}, \quad \mathbf{b} = \begin{bmatrix} h\_1\\ \frac{h\_1}{2}f\_0 + \kappa\_0 a\\ h\_1 + h\_2\\ \vdots\\ h\_{N-1} + h\_N\\ \frac{h\_{N-1} + h\_N}{2}f\_{N-1}\\ \frac{h\_N}{2}f\_N - \kappa\_1 \boldsymbol{\beta} \end{bmatrix}.$$
