**4. The FE method for the heat equation**

Consider the following heat/diffusion problem: Find *u*ð Þ **x**, *t* such that

$$
\dot{u} - \Delta u = f, \quad \text{in } \Omega \subset \mathbb{R}^2, \quad t \in [0, T], \tag{48}
$$

$$u(\cdot, t) = 0, \quad \text{on } \partial \Omega \text{ and } t \in [0, T], \tag{49}$$

$$u(\mathbf{x},0) = u\_0(\mathbf{x}), \quad \text{for } \mathbf{x} \in \Omega \text{ and } t = \mathbf{0}. \tag{50}$$

We seek a weak solution *<sup>u</sup>* in *<sup>V</sup>*<sup>0</sup> <sup>¼</sup> *v v* j j k k <sup>þ</sup> k k <sup>∇</sup>*<sup>v</sup>* <sup>&</sup>lt; <sup>∞</sup>, *<sup>v</sup>* <sup>∂</sup><sup>Ω</sup> <sup>¼</sup> <sup>0</sup> � �. In order to derive the weak formulation, we multiply (48) with *v*∈*V*0, integrate over Ω and use Green's formula to obtain, for *t*∈½ � 0, *T* ,

$$\int\_{\Omega} f v d\mathbf{x} = \int\_{\Omega} i v d\mathbf{x} + \int\_{\Omega} \nabla u \cdot \nabla v d\mathbf{x} - \int\_{\partial \Omega} v \nabla u \cdot \mathbf{n} ds = \int\_{\Omega} i v d\mathbf{x} + \int\_{\Omega} \nabla u \cdot \nabla v d\mathbf{x}.$$

The weak form therefore reads: Find *u*ð Þ �, *t* ∈*V*<sup>0</sup> such that for *t*>0

$$\int\_{\Omega} \dot{u}v d\mathbf{x} + \int\_{\Omega} \nabla u \cdot \nabla v d\mathbf{x} = \int\_{\Omega} fv d\mathbf{x}, \quad v \in V\_0. \tag{51}$$

The semi-discrete FE method in 2D is defined as follows: Find *uh*ð Þ �, *t* ∈*Vh*,0 ¼ *<sup>v</sup>*∈*Vh*j j *<sup>v</sup>* <sup>∂</sup><sup>Ω</sup> <sup>¼</sup> <sup>0</sup> � � such that

$$\int\_{\Omega} \dot{u}\_h v\_h d\mathbf{x} + \int\_{\Omega} \nabla u\_h \cdot \nabla v\_h d\mathbf{x} = \int\_{\Omega} f v\_h d\mathbf{x}, \quad v\_h \in V\_{h,0}, \tag{52}$$

where *Vh* <sup>¼</sup> *<sup>v</sup>*∈*V v* j j*<sup>T</sup>* <sup>∈</sup>*P*<sup>1</sup> ð Þ *T* , ∀ *T* ∈ T *<sup>h</sup>* � �.

**Implementation**: Substituting *uh*ð Þ¼ **<sup>x</sup>**, *<sup>t</sup>* <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>1</sup>*cj*ð Þ*<sup>t</sup> <sup>ϕ</sup><sup>j</sup>* ð Þ **x** into (52) and choosing *vh* ¼ *ϕi*, we obtain

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

This gives us the system of ODEs

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

where **<sup>c</sup>** <sup>¼</sup> ½ � *<sup>c</sup>*1,*c*2, … ,*cN <sup>t</sup>* <sup>¼</sup> ½ � *uh*ð Þ *<sup>N</sup>*1, *<sup>t</sup>* , … , *uh*ð Þ *NN*, *<sup>t</sup> <sup>t</sup>* ∈ *<sup>N</sup>* (here *Ni* denotes the node that belongs to the basis function *ϕi*) is the unknown vector and the entries of *M*, *A* ∈ *<sup>N</sup>*�*<sup>N</sup>* and **b**∈ *<sup>N</sup>* are given by

$$m\_{\vec{\eta}} = \int\_{\Omega} \phi\_i \phi\_j d\mathbf{x}, \quad a\_{\vec{\eta}} = \int\_{\Omega} \nabla \phi\_j \cdot \nabla \phi\_i d\mathbf{x}, \quad b\_i = \int\_{\Omega} f \phi\_i d\mathbf{x}, \quad i, j = 1, 2, \dots, N.$$

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

Finally, the system of ODEs can be solved with *e:g:*, the backward Euler method as follows: Let 0 ¼ *t*<sup>0</sup> < *t*<sup>1</sup> < ⋯ <*tM* ¼ *T* be a discretization, let *km* ¼ *tm* � *tm*�<sup>1</sup> for *<sup>m</sup>* <sup>¼</sup> 1, 2, … , *<sup>M</sup>* be the time step size and let **<sup>c</sup>***<sup>m</sup>* <sup>≈</sup>**c**ð Þ *tm* for *<sup>m</sup>* <sup>¼</sup> 1, 2, … , *<sup>M</sup>* denote corresponding approximations. Then, we can compute **c***<sup>m</sup>* using

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

where **<sup>c</sup>**<sup>0</sup> is obtained from *<sup>u</sup>*0ð Þ **<sup>x</sup>** . We can either use **<sup>c</sup>**<sup>0</sup> <sup>¼</sup> ½ � *<sup>c</sup>*1ð Þ <sup>0</sup> , … ,*cN*ð Þ <sup>0</sup> *<sup>t</sup>* <sup>¼</sup> ½ � *<sup>u</sup>*0ð Þ *<sup>N</sup>*<sup>1</sup> , … , *<sup>u</sup>*0ð Þ *NN <sup>t</sup>* , or we can let **c**<sup>0</sup> to be the *L*<sup>2</sup> -projection of *u*0. We set *u*<sup>0</sup> *h* ¼ P*<sup>N</sup> <sup>j</sup>*¼1*c*<sup>0</sup> *<sup>j</sup>ϕj*ð Þ **<sup>x</sup>** and solve for *<sup>c</sup>*<sup>0</sup> *<sup>j</sup>* using

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

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

$$\|\|u(\cdot,t)\|\| \le \|u(\cdot,0)\|\| + \int\_0^t \|\|f(\cdot,s)\|\|ds, \quad \left||u\_h^m\right|| \le \left||u\_h^{m-1}\right|| + k\_m \left||f\_m\right|| \le \left||u\_h^0\right|| + \sum\_{i=1}^m k\_i \left||f\_i\right||.$$

#### **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*ð Þ **x**, *t* such that

$$
\ddot{u} - \nabla \cdot (\varepsilon \nabla u) = f, \quad \text{in } \Omega \subset \mathbb{R}^2, \quad t \in [0, T], \tag{53}
$$

$$\mathbf{n} \cdot \nabla u(\cdot, t) = 0, \quad \text{on } \partial \Omega \text{ and } t \in [0, T], \tag{54}$$

$$u(\mathbf{x},0) = u\_0(\mathbf{x}), \quad \dot{u}(\mathbf{x},0) = v\_0(\mathbf{x}), \quad \text{for } \mathbf{x} \in \Omega \text{ and } t = 0,\tag{55}$$

where *f* is a given load, *ε* ¼ *ε*ð Þ **x**, *t* is a positive parameter, *u*<sup>0</sup> and *v*<sup>0</sup> are a prescribed initial conditions, and Ω is a bounded domain with boundary ∂Ω and unit outward normal **n**.

We seek a weak solution *<sup>u</sup>* in *<sup>V</sup>* <sup>¼</sup> *<sup>H</sup>*<sup>1</sup> ð Þ¼ Ω f g *v*j k k*v* þ k k ∇*v* < ∞ . Multiplying the wave Eq. (53) with *v*∈*V*, integrating over Ω, and using Green's formula, we obtain, for *t*∈ ½ � 0, *T* ,

$$\begin{aligned} \int\_{\Omega} f v d\mathbf{x} &= \int\_{\Omega} \ddot{u} v d\mathbf{x} - \int\_{\Omega} v \nabla \cdot (e \nabla u) d\mathbf{x} = \int\_{\Omega} \ddot{u} v d\mathbf{x} + \int\_{\Omega} e \nabla u \cdot \nabla v d\mathbf{x} - \int\_{\partial \Omega} v e \nabla u \cdot \mathbf{n} ds \\ &= \int\_{\Omega} \ddot{u} v d\mathbf{x} + \int\_{\Omega} e \nabla u \cdot \nabla v d\mathbf{x}. \end{aligned}$$

The weak form (variational formulation) therefore reads: Find *u*ð Þ �, *t* ∈*V* ¼ *H*1 ð Þ Ω such that for all *t*>0

$$
\int\_{\Omega} \ddot{u}v d\mathbf{x} + \int\_{\Omega} \varepsilon \nabla u \cdot \nabla v d\mathbf{x} = \int\_{\Omega} f v d\mathbf{x}, \quad v \in V. \tag{56}
$$

Let *Vh* <sup>¼</sup> *<sup>v</sup>*∈*V v* j j*<sup>T</sup>* <sup>∈</sup>*P*<sup>1</sup> ð Þ *T* , ∀ *T* ∈ T *<sup>h</sup>* � �⊂*V* be the space of all continuous piecewise linear functions on a triangle mesh of Ω. The semi-discrete FE method in 2D is defined as follows: Find *uh*ð Þ �, *t* ∈*Vh* such that

$$\int\_{\Omega} \ddot{u}\_h v\_h d\mathbf{x} + \int\_{\Omega} \varepsilon \nabla u\_h \cdot \nabla v\_h d\mathbf{x} = \int\_{\Omega} f v\_h d\mathbf{x}, \quad v\_h \in V\_h. \tag{57}$$

**Implementation**: Substituting *uh*ð Þ¼ **<sup>x</sup>**, *<sup>t</sup>* <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>1</sup>*cj*ð Þ*<sup>t</sup> <sup>ϕ</sup><sup>j</sup>* ð Þ **x** into (57) and choosing *vh* ¼ *ϕi*, we obtain

$$\sum\_{j=1}^{N} \ddot{\boldsymbol{c}}\_{j} \int\_{\Omega} \phi\_{j} \phi\_{i} d\mathbf{x} + \sum\_{j=1}^{N} \boldsymbol{c}\_{j} \int\_{\Omega} \boldsymbol{\varepsilon} \nabla \phi\_{j} \cdot \nabla \phi\_{i} d\mathbf{x} = \int\_{\Omega} f \phi\_{i} d\mathbf{x}, \quad i = 1, 2, \dots, N.$$

This gives us the system

$$M\ddot{\mathbf{c}}(t) + A(t)\mathbf{c}(t) = \mathbf{b}(t), \quad t \in (0, T], \tag{58}$$

where **<sup>c</sup>** <sup>¼</sup> ½ � *<sup>c</sup>*1, … ,*cN <sup>t</sup>* <sup>¼</sup> ½ � *uh*ð Þ *<sup>N</sup>*1, *<sup>t</sup>* , … , *uh*ð Þ *NN*, *<sup>t</sup> <sup>t</sup>* ∈ *<sup>N</sup>* (here *Ni* denotes the node that belongs to the basis function *ϕi*) is the unknown vector and the entries of the mass and stiffness matrices *M*, *A* ∈ *<sup>N</sup>*�*<sup>N</sup>* and the load vector **b**∈ *<sup>N</sup>* are given by

$$m\_{\vec{\eta}} = \int\_{\Omega} \phi\_i \phi\_j d\mathbf{x}, \quad a\_{\vec{\eta}} = \int\_{\Omega} \varepsilon \nabla \phi\_j \cdot \nabla \phi\_i d\mathbf{x}, \quad b\_i = \int\_{\tilde{\Omega}} f \phi\_i d\mathbf{x}, \quad 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 **d**ðÞ¼ *t* **c**\_ð Þ*t* , we get the new coupled system

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

Let **<sup>w</sup>** <sup>¼</sup> ½ � **<sup>c</sup>**, **<sup>d</sup>** *<sup>t</sup>* then the system is equivalent to *<sup>M</sup>*^ **<sup>w</sup>**\_ ðÞþ*<sup>t</sup> A t* ^ð Þ**w**ðÞ¼ *<sup>t</sup>* **<sup>b</sup>**^ð Þ*<sup>t</sup>* , *t*∈ð � 0, *T* , where

$$
\hat{M} = \begin{bmatrix} M & \mathbf{0} \\ \mathbf{0} & M \end{bmatrix}, \quad \hat{A} = \begin{bmatrix} \mathbf{0} & -M \\ A & \mathbf{0} \end{bmatrix}, \quad \hat{\mathbf{b}} = \begin{bmatrix} \mathbf{0} \\ \mathbf{b} \end{bmatrix}.
$$

Finally, the system of ODEs can be solved with *e:g:*, the backward Euler method as follows: Let 0 ¼ *t*<sup>0</sup> < *t*<sup>1</sup> < ⋯ <*tM* ¼ *T* be a discretization, let *km* ¼ *tm* � *tm*�<sup>1</sup> for *<sup>m</sup>* <sup>¼</sup> 1, 2, … , *<sup>M</sup>* be the time step size and let **<sup>w</sup>***<sup>m</sup>* <sup>≈</sup> **<sup>w</sup>**ð Þ *tm* for *<sup>m</sup>* <sup>¼</sup> 1, 2, … , *<sup>M</sup>* denote corresponding approximations. Then, we can compute **w***<sup>m</sup>* using

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

where **<sup>w</sup>**<sup>0</sup> is obtained from *<sup>u</sup>*0ð Þ **<sup>x</sup>** and *<sup>v</sup>*0ð Þ **<sup>x</sup>** .

There are several possible choices of initial data. We can either use **<sup>w</sup>**<sup>0</sup> <sup>¼</sup> ½ � *<sup>w</sup>*1ð Þ <sup>0</sup> , … ,*c*2*<sup>N</sup>*ð Þ <sup>0</sup> *<sup>t</sup>* <sup>¼</sup> ½ � *<sup>u</sup>*0ð Þ *<sup>N</sup>*<sup>1</sup> , … , *<sup>u</sup>*0ð Þ *NN* , *<sup>v</sup>*0ð Þ *<sup>N</sup>*<sup>1</sup> , … , *<sup>v</sup>*0ð Þ *NN <sup>t</sup>* , or we can let **<sup>w</sup>**<sup>0</sup> <sup>¼</sup> **w**<sup>0</sup> <sup>1</sup> , **w**<sup>0</sup> 2 � �*<sup>t</sup>* , where **w**<sup>0</sup> <sup>1</sup> and **w**<sup>0</sup> <sup>2</sup> are the *L*<sup>2</sup> -projection of *u*<sup>0</sup> and *v*0, respectively. We set *w*<sup>0</sup> *<sup>h</sup>*,1 <sup>¼</sup> <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>1</sup>*w*<sup>0</sup> *<sup>j</sup>*,1*ϕ<sup>j</sup>* ð Þ **<sup>x</sup>** and *<sup>w</sup>*<sup>0</sup> *<sup>h</sup>*,2 <sup>¼</sup> <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>1</sup>*w*<sup>0</sup> *<sup>j</sup>*,2*ϕ<sup>j</sup>* ð Þ **<sup>x</sup>** and solve for *<sup>w</sup>*<sup>0</sup> *<sup>j</sup>*,1, *w*<sup>0</sup> *<sup>j</sup>*,2 using

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

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

We can also use Crank–Nicolson scheme

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

**Theorem 5.1 (Conservation of energy)** *If f* ¼ 0, *then*

$$\left\|\dot{\boldsymbol{u}}\_{h}(\cdot,t)\right\|\_{L^{2}(\mathfrak{Q})}^{2} + \varepsilon \left\|\nabla\boldsymbol{u}\_{h}(\cdot,t)\right\|\_{L^{2}(\mathfrak{Q})}^{2} = \left\|\dot{\boldsymbol{u}}(\cdot,\mathbb{O})\right\|\_{L^{2}(\mathfrak{Q})}^{2} + \varepsilon \left\|\nabla\boldsymbol{u}(\cdot,\mathbb{O})\right\|\_{L^{2}(\mathfrak{Q})}^{2}.$$
