**3.4** *L***2-projection**

Let Ω ⊂ <sup>2</sup> . We consider the space *L*<sup>2</sup> ð Þ¼ Ω *v*j Ð <sup>Ω</sup>*v*<sup>2</sup>ð Þ *<sup>x</sup>*, *<sup>y</sup> dxdy*<sup>&</sup>lt; <sup>∞</sup> � �. Let *u*∈ *L*<sup>2</sup> ð Þ <sup>Ω</sup> *:* We define the *<sup>L</sup>*<sup>2</sup> -projection *Ph* : *L*<sup>2</sup> ð Þ! Ω *Vh* ¼ *<sup>v</sup>*∈*C*<sup>0</sup>ð Þ <sup>Ω</sup> j j *<sup>v</sup> <sup>T</sup>* <sup>∈</sup>*P*<sup>1</sup> ð Þ *T* , ∀ *T* ∈ T *<sup>h</sup>* � � by *Phu* ∈*Vh* such that

$$\int\_{\Omega} (u - P\_h u) v\_h d\mathbf{x} d\mathbf{y} = \mathbf{0}, \quad \forall v\_h \in V\_h.$$

The problem of finding *Phu* ∈*Vh* is equivalent to solve the following linear system

$$\int\_{\Omega} (u - P\_h u) \phi\_i d\mathbf{x} dy = 0, \quad i = 1, 2 \dots, N,$$

where f g *ϕ<sup>i</sup> N <sup>i</sup>*¼<sup>1</sup> is a basis of *Vh*.

Since *Phu*∈*Vh* we can express it as *Phu* <sup>¼</sup> <sup>P</sup>*<sup>N</sup> <sup>i</sup>*¼<sup>1</sup>*ciϕi*ð Þ *<sup>x</sup>*, *<sup>y</sup>* , where *ci* <sup>∈</sup> . Therefore, to find *Phu* ∈*Vh* we need to find *c*1,*c*2, … ,*cN* ∈ such that

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

The problem can be expressed as a linear system of equations *M***c** ¼ **b**, where **<sup>c</sup>** <sup>¼</sup> ½ � *<sup>c</sup>*1,*c*2, … ,*cN <sup>t</sup>* and the entries of the matrix *<sup>M</sup>* <sup>∈</sup> *<sup>N</sup>*�*<sup>N</sup>* and the vector **<sup>b</sup>**<sup>∈</sup> *<sup>N</sup>* are given by

$$m\_{\vec{\eta}} = \int\_{\Omega} \phi\_i \phi\_j d\mathbf{x} dy, \quad b\_j = \int\_{\Omega} u \phi\_j d\mathbf{x} dy.$$

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

$$\int\_{T} f(\mathbf{x}, \mathbf{y}) d\mathbf{x} d\mathbf{y} \approx \sum\_{j=1}^{n} w\_j f\left(\overline{N}\_j\right),$$

where the *ω<sup>j</sup>* 0 s denote the weights and the *N <sup>j</sup>* � � <sup>0</sup> s the quadrature points.

**Lemma 3.1** *The mass matrix M with entries mij* <sup>¼</sup> <sup>Ð</sup> <sup>Ω</sup>*ϕiϕjdxdy is symmetric and positive definite.*

**Theorem 3.3** *For any u*∈ *L*<sup>2</sup> ð Þ <sup>Ω</sup> *the L*<sup>2</sup> -*projection Phu exists and is unique*.

#### **3.5 A priori error estimate**

**Theorem 3.4** *Let u*∈ *L*<sup>2</sup> ð Þ <sup>Ω</sup> *and let Phu be the L*<sup>2</sup> -*projection of u*, *then*

$$||\mathfrak{u} - P\_h \mathfrak{u}||\_{L^2(\mathfrak{Q})} \le ||\mathfrak{u} - \mathfrak{v}\_h||\_{L^2(\mathfrak{Q})}, \quad \forall \ \mathfrak{v}\_h \in V\_h.$$

**Theorem 3.5** *Suppose that u* ∈*C*<sup>2</sup> ð Þ <sup>Ω</sup> *with u*∈*C*<sup>2</sup> ð Þ *T for all T* ∈T *<sup>h</sup>*. *Then there exists a constant C such that*

$$\left\| \|u - P\_h u\|\|\_{L^2(\Omega)}^2 \le C \sum\_{T \in \mathcal{T}\_h} h\_T^4 \left\| D^2 u \right\|\|\_{L^2(T)}^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 (a\nabla u) + bu = f, \quad \text{in } \Omega, \quad a\nabla u \cdot \mathbf{n} = \kappa(\mathbf{g} - u), \quad \text{on } \partial\Omega,\tag{31}$$

where *a*> 0, *b*≥ 0, *κ* ≥ 0, *f* ∈*L*<sup>2</sup> ð Þ <sup>Ω</sup> and *<sup>g</sup>* <sup>∈</sup>*C*<sup>0</sup>ð Þ <sup>∂</sup><sup>Ω</sup> . We seek a weak solution *<sup>u</sup>* in

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

In order to derive the weak formulation, we multiply (31) with *v*∈*V*, integrate over Ω and use Green's formula to obtain

$$\begin{aligned} \int\_{\Omega} f v \mathrm{d}x dy &= -\int\_{\Omega} v \nabla \cdot (a \nabla u) dx dy + \int\_{\Omega} b u v d\mathrm{x} dy \\ &= \int\_{\Omega} a \nabla u \cdot \nabla v d\mathrm{x} dy - \int\_{\partial \Omega} v (a \nabla u) \cdot \mathbf{n} ds + \int\_{\Omega} b u v d\mathrm{x} dy \\ &= \int\_{\Omega} a \nabla u \cdot \nabla v d\mathrm{x} dy + \int\_{\Omega} b u v d\mathrm{x} dy + \int\_{\partial \Omega} \kappa (u - \mathrm{g}) v d\mathrm{s}. \end{aligned}$$

We obtain the weak form: Find *u*∈*V* such that

$$\int\_{\mathfrak{Q}} a \nabla u \cdot \nabla v \mathrm{d}x dy + \int\_{\mathfrak{Q}} buv \mathrm{d}x dy + \int\_{\partial \mathfrak{Q}} \kappa uv ds = \int\_{\tilde{\Omega}} f v \mathrm{d}x dy + \int\_{\partial \mathfrak{Q}} \kappa y ds, \quad v \in V. \tag{32}$$

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

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 *uh* ∈*Vh* such that

$$
\int\_{\Omega} a \nabla u\_h \cdot \nabla v\_h dx dy + \int\_{\Omega} b u\_h v\_h dx dy + \int\_{\partial \Omega} \kappa u\_h v\_h ds = \int\_{\Omega} f v\_h dx dy + \int\_{\partial \Omega} \kappa g v\_h ds, \quad v\_h \in V\_h,\tag{33}
$$

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**: Let *<sup>a</sup>* <sup>¼</sup> 1 and *<sup>b</sup>* <sup>¼</sup> *<sup>g</sup>* <sup>¼</sup> 0. Substituting *uh* <sup>¼</sup> <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>1</sup>*cjϕ<sup>j</sup>* into (33) and picking *vh* ¼ *ϕi*, we obtain

$$\sum\_{j=1}^{N} c\_j \left( \int\_{\Omega} \nabla \phi\_j \cdot \nabla \phi\_i dx dy + \int\_{\partial \Omega} \kappa \phi\_j \phi\_i ds \right) = \int\_{\Omega} f \phi\_i dx dy, \quad i = 1, 2, \dots, N.$$

This gives us the system ð Þ *<sup>A</sup>* <sup>þ</sup> *<sup>R</sup>* **<sup>c</sup>** <sup>¼</sup> **<sup>b</sup>**, where **<sup>c</sup>** <sup>¼</sup> ½ � *<sup>c</sup>*1,*c*2, … ,*cN <sup>t</sup>* ∈ *<sup>N</sup>* is the unknown vector and the entries of *A* ∈ *<sup>N</sup>*�*<sup>N</sup>*, *R*∈ *<sup>N</sup>*�*<sup>N</sup>*, and **b**∈ *<sup>N</sup>* are given by

$$\mathcal{A}\_{\vec{\eta}} = \int\_{\Omega} \nabla \phi\_j \cdot \nabla \phi\_i d\mathbf{x} dy, \quad r\_{\vec{\eta}} = \int\_{\partial \Omega} \kappa \phi\_j \phi\_i d\mathbf{s}, \quad b\_i = \int\_{\Omega} f \phi\_i d\mathbf{x} dy, \quad 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\_{ij}^T = \int\_{\Omega} \nabla \phi\_j \cdot \nabla \phi\_i d\mathbf{x} dy, \quad i, j = 1, 2, 3.$$

where *T* is an arbitrary triangle with vertices *Ni* ¼ *xi*, *yi* � � and *ϕ<sup>i</sup>* are the hat functions *i:e:*, *ϕ<sup>j</sup>* ð Þ¼ *Ni δij*. Let *ϕi*ð Þ¼ *x*, *y α<sup>i</sup>* þ *βix* þ *γiy*, for *i* ¼ 1, 2, 3. Then, we compute *αi*, *βi*, *γ<sup>i</sup>* by

$$
\begin{bmatrix} \mathbf{1} & \mathbf{x}\_1 & \mathbf{y}\_1 \\ \mathbf{1} & \mathbf{x}\_2 & \mathbf{y}\_2 \\ \mathbf{1} & \mathbf{x}\_3 & \mathbf{y}\_3 \end{bmatrix} \begin{bmatrix} a\_1 \\ \boldsymbol{\beta}\_1 \\ \boldsymbol{\gamma}\_1 \end{bmatrix} = \begin{bmatrix} \mathbf{1} \\ \mathbf{0} \\ \mathbf{0} \end{bmatrix}, \quad \begin{bmatrix} \mathbf{1} & \mathbf{x}\_1 & \mathbf{y}\_1 \\ \mathbf{1} & \mathbf{x}\_2 & \mathbf{y}\_2 \\ \mathbf{1} & \mathbf{x}\_3 & \mathbf{y}\_3 \end{bmatrix} \begin{bmatrix} a\_2 \\ \boldsymbol{\beta}\_2 \\ \boldsymbol{\gamma}\_2 \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \mathbf{1} \\ \mathbf{0} \end{bmatrix}, \quad \begin{bmatrix} \mathbf{1} & \mathbf{x}\_1 & \mathbf{y}\_1 \\ \mathbf{1} & \mathbf{x}\_2 & \mathbf{y}\_2 \\ \mathbf{1} & \mathbf{x}\_3 & \mathbf{y}\_3 \end{bmatrix} \begin{bmatrix} \boldsymbol{\alpha}\_3 \\ \boldsymbol{\beta}\_3 \\ \boldsymbol{\gamma}\_3 \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \\ \mathbf{1} \end{bmatrix}.
$$

In general we have *Bα<sup>i</sup>* ¼ **e***<sup>i</sup>* for *i* ¼ 1, 2, 3, where

$$B = \begin{bmatrix} \mathbf{1} & \mathbf{x}\_1 & \mathbf{y}\_1 \\ \mathbf{1} & \mathbf{x}\_2 & \mathbf{y}\_2 \\ \mathbf{1} & \mathbf{x}\_3 & \mathbf{y}\_3 \end{bmatrix}, \quad a\_i = \begin{bmatrix} a\_i \\ \boldsymbol{\beta}\_i \\ \boldsymbol{\gamma}\_i \end{bmatrix}, \quad \mathbf{e}\_1 = \begin{bmatrix} \mathbf{1} \\ \mathbf{0} \\ \mathbf{0} \end{bmatrix}, \quad \mathbf{e}\_2 = \begin{bmatrix} \mathbf{0} \\ \mathbf{1} \\ \mathbf{0} \end{bmatrix}, \quad \mathbf{e}\_3 = \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \\ \mathbf{1} \end{bmatrix}.$$

Furthermore, we obviously have <sup>∇</sup>*ϕ<sup>i</sup>* <sup>¼</sup> *<sup>β</sup>i*, *<sup>γ</sup><sup>i</sup>* ½ �*<sup>t</sup>* , which gives

$$a\_{ij}^T = \int\_{\Omega} \left(\beta\_i \beta\_j + \gamma\_i \gamma\_j\right) d\mathbf{x} = \left(\beta\_i \beta\_j + \gamma\_i \gamma\_j\right)|T|, \quad i, j = 1, 2, 3.$$

**Assembly of boundary matrix** *R*: Let Γ*out <sup>h</sup>* denote the set of boundary edges of the triangulation, *i:e:* Γ*out <sup>h</sup>* ¼ f g *E*j *E* ¼ *T* ∩ ∂Ω, *for T* ∈T *<sup>h</sup>* . Assume that *κ* is constant on *E*. For an edge *E*∈Γ*out <sup>h</sup>* , we define *RE* ∈ <sup>2</sup>�<sup>2</sup> by the entries

*Finite Element Methods and Their Applications*

$$r\_{\vec{\eta}}^E = \int\_E \kappa \phi\_j \phi\_i ds = \frac{\kappa}{6} \left( 1 + \delta\_{\vec{\eta}} \right) |E|, \quad i, j = 1, 2, 3$$

where ∣*E*∣ is the length of *E* and *δij* is 1 for *i* ¼ *j* and 0 else.

**Assembly of load vector**: We use a corner quadrature rule for approximating the integral. We obtain for *T* ∈T *<sup>h</sup>*

$$b\_i^T = \int\_T f \phi\_i d\mathbf{x} dy \approx \frac{|T|}{3} f(N\_i), \quad i = 1, 2, \dots, N.$$

Given *<sup>A</sup>*, *<sup>R</sup>* and **<sup>b</sup>**, we can solve ð Þ *<sup>A</sup>* <sup>þ</sup> *<sup>R</sup>* **<sup>c</sup>** <sup>¼</sup> **<sup>b</sup>** and write *uh* <sup>¼</sup> <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>1</sup>*cjϕ<sup>j</sup>* .

#### **3.7 The Dirichlet problem**

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

$$- \Delta u = f, \quad \text{in } \Omega, \quad u = \emptyset, \quad \text{on } \partial \Omega,\tag{34}$$

where *f* ∈*L*<sup>2</sup> ð Þ <sup>Ω</sup> and *<sup>g</sup>* <sup>∈</sup>*C*<sup>0</sup>ð Þ <sup>∂</sup><sup>Ω</sup> . We seek a weak solution *<sup>u</sup>* in *Vg* <sup>¼</sup> *<sup>v</sup>*∈*V v* j j<sup>∂</sup><sup>Ω</sup> <sup>¼</sup> *<sup>g</sup>* � �*:* Multiplying (34) by a test function *<sup>v</sup>*∈*V*<sup>0</sup> and integrating over <sup>Ω</sup>, we get

$$
\int\_{\Omega} f v \mathrm{d}x \mathrm{d}y = -\int\_{\Omega} v \Delta u \mathrm{d}x \mathrm{d}y = \int\_{\Omega} \nabla u \cdot \nabla v \mathrm{d}x \mathrm{d}y - \int\_{\partial \Omega} v \nabla u \cdot \mathbf{n} \mathrm{d}s = \int\_{\Omega} \nabla u \cdot \nabla v \mathrm{d}x \mathrm{d}y.
$$

So the weak problem reads: Find *u*∈*Vg* such that

$$\int\_{\Omega} \nabla u \cdot \nabla v \mathrm{d}x dy = \int\_{\Omega} f v \mathrm{d}x dy, \quad v \in V\_0.$$

Assume that *g* is piecewise linear on ∂Ω with respect to the triangulation. Then the FE method in 2D is defined as follows: Find *uh* <sup>∈</sup>*Vh*,*<sup>g</sup>* <sup>¼</sup> *<sup>v</sup>*∈*Vh*j j *<sup>v</sup>* <sup>∂</sup><sup>Ω</sup> <sup>¼</sup> *<sup>g</sup>* � � such that

$$\int\_{\Omega} \nabla u\_h \cdot \nabla v\_h d\mathbf{x} dy = \int\_{\Omega} f v\_h d\mathbf{x} dy, \quad v\_h \in V\_{h,0} \dots$$

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

$$
\begin{bmatrix} A\_{0,0} & A\_{0, \emptyset} \\ A\_{\emptyset,0} & A\_{0, \emptyset} \end{bmatrix} \begin{bmatrix} \mathbf{c}\_{0} \\ \mathbf{c}\_{1} \end{bmatrix} = \begin{bmatrix} \mathbf{b}\_{0} \\ \mathbf{b}\_{1} \end{bmatrix},
$$

where *A*0,0 ∈ ð Þ� *<sup>N</sup>*�*<sup>J</sup>* ð Þ *<sup>N</sup>*�*<sup>J</sup>* , *Ag*,*<sup>g</sup>* ∈ *<sup>J</sup>*�*<sup>J</sup>* , *A*0,*<sup>g</sup>* ∈ ð Þ� *<sup>N</sup>*�*<sup>J</sup> <sup>J</sup>* , *Ag*,0 ∈ *<sup>J</sup>*�ð Þ *<sup>N</sup>*�*<sup>J</sup>* . Note that **c**<sup>1</sup> ∈ *<sup>J</sup>* is known (it contains the values of *g* in the boundary nodes). We can therefore solve the simplified problem reading: find **<sup>c</sup>**<sup>0</sup> <sup>∈</sup> *<sup>N</sup>*�*<sup>J</sup>* with *<sup>A</sup>*0,0**c**<sup>0</sup> <sup>¼</sup> **<sup>b</sup>**<sup>0</sup> � *<sup>A</sup>*0,*<sup>g</sup>***c**1.

#### **3.8 The Neumann problem**

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

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

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

where *f* ∈*L*<sup>2</sup> ð Þ <sup>Ω</sup> and *<sup>g</sup>* <sup>∈</sup>*C*0ð Þ <sup>∂</sup><sup>Ω</sup> . Let us try to seek a solution to this problem in the space *<sup>V</sup>* <sup>¼</sup> *<sup>v</sup>*<sup>j</sup> k k*<sup>v</sup> <sup>L</sup>*2ð Þ <sup>Ω</sup> <sup>þ</sup> k k <sup>∇</sup>*<sup>v</sup> <sup>L</sup>*2ð Þ <sup>Ω</sup> <sup>&</sup>lt; <sup>∞</sup> n o. Multiplying (35) by a test function *v*∈*V*, integrating over Ω, and using Green's formula, we get

$$\begin{aligned} \int\_{\mathfrak{U}} f v \mathrm{d}x dy &= -\int\_{\mathfrak{U}} v \Delta u \mathrm{d}x dy = \int\_{\mathfrak{U}} \nabla u \cdot \nabla v \mathrm{d}x dy - \int\_{\partial \mathfrak{U}} v \nabla u \cdot \mathbf{n} ds \\ &= \int\_{\mathfrak{U}} \nabla u \cdot \nabla v \mathrm{d}x dy - \int\_{\partial \mathfrak{U}} v g ds. \end{aligned}$$

Thus, the variational formulation reads: find *u*∈*V* such that

$$\int\_{\Omega} \nabla u \cdot \nabla v \mathrm{d}x dy - \int\_{\partial \Omega} v g ds = \int\_{\Omega} f v d\mathbf{x} dy, \quad \forall \ v \in V.$$

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

$$0 = \int\_{\Omega} \nabla u \cdot \nabla \mathbf{1} dx dy = \int\_{\Omega} f d\mathbf{x} dy + \int\_{\partial \Omega} g ds \, \mathbf{J}$$

Therefore we need to assume the following compatibility condition

$$\int\_{\Omega} f d\mathbf{x} dy + \int\_{\partial\Omega} g ds = 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*∈ . To fix this constant and obtain a unique solution a common trick is to impose the additional constraint Ð <sup>Ω</sup>*udxdy* ¼ 0. We therefore define the weak solution space

$$
\hat{V} = \left\{ v \in V \mid \int\_{\Omega} v dx dy = 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*∈*V*^ such that

$$\int\_{\Omega} \nabla u \cdot \nabla v \mathrm{d}x \mathrm{d}y - \int\_{\partial \Omega} v g \mathrm{d}s = \int\_{\Omega} f v \mathrm{d}x \mathrm{d}y, \quad \forall \ v \in V.$$

Now, the FE method takes the form: find *uh* ∈*V*^ *<sup>h</sup>* ⊂*V*^ such that

$$\int\_{\Omega} \nabla u\_h \cdot \nabla v\_h d\mathbf{x} dy - \int\_{\partial \Omega} v\_h \mathbf{g} ds = \int\_{\Omega} f v\_h d\mathbf{x} dy, \quad \forall \ v\_h \in \hat{V}\_h, \mathbf{x}$$

where *V*^ *<sup>h</sup>* 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 twodimensional finite elements. Suppose Ω is a domain in either <sup>2</sup> or <sup>3</sup> and assume that ∂Ω has been partitioned into two disjoint sets: ∂Ω ¼ Γ<sup>1</sup> ∪ Γ2. We consider the following BVP:

$$-\nabla \cdot (\kappa(\mathbf{x}) \nabla u) = f(\mathbf{x}), \quad \mathbf{x} \in \ \Omega, \quad u = 0, \quad \mathbf{x} \in \ \Gamma\_1, \quad \nabla u \cdot \mathbf{n} = 0, \quad \mathbf{x} \in \ \Gamma\_2,\tag{36}$$

where *f* ∈*L*<sup>2</sup> ð Þ Ω . 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

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

Multiplying (36) by a test function *v*∈*V*^ and integrating over Ω, we get

$$\int\_{\Omega} f v \mathrm{d}x dy = -\int\_{\Omega} v \nabla \cdot (\kappa(\mathbf{x}) \nabla u) d\mathbf{x} dy = \int\_{\partial \Omega} \kappa(\mathbf{x}) \nabla u \cdot \nabla v d\mathbf{x} dy - \int\_{\partial \Omega} \kappa(\mathbf{x}) v \nabla u \cdot \mathbf{n} ds$$

$$= \int\_{\Omega} \kappa(\mathbf{x}) \nabla u \cdot \nabla v d\mathbf{x} dy - \int\_{\Gamma\_1} \kappa(\mathbf{x}) v \nabla u \cdot \mathbf{n} ds - \int\_{\Gamma\_2} \kappa(\mathbf{x}) v \nabla u \cdot \mathbf{n} ds$$

$$= \int\_{\Omega} \kappa(\mathbf{x}) \nabla u \cdot \nabla v d\mathbf{x} dy,$$

since *<sup>v</sup>* <sup>¼</sup> 0 on <sup>Γ</sup><sup>1</sup> and <sup>∇</sup>*<sup>u</sup>* � **<sup>n</sup>** on <sup>Γ</sup>1. Thus the weak form of (36) is: Find *<sup>u</sup>*∈*V*^ such that

$$\int\_{\Omega} \kappa(\mathbf{x}) \nabla u \cdot \nabla v \mathrm{d}x \mathrm{d}y = \int\_{\Omega} f v \mathrm{d}x \mathrm{d}y, \quad v \in \hat{V}. \tag{37}$$

We now restrict our discussion once more to two-dimensional polygonal domains. To apply the FE method, we must choose an approximating subspace of *V*^ . Since the boundary conditions are mixed, there are at least two points where the boundary conditions change from Dirichlet to Neumann. We will make the assumption that the mesh is chosen so that all such points are nodes (and that all such nodes belong to Γ1, that is, that Γ<sup>1</sup> includes its "endpoints"). We can then choose the approximating subspace of *V*^ as follows:

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

A basis for *Vh* is formed by including all basis functions corresponding to interior boundary nodes that do not belong to Γ1. If the BVP includes only Neumann conditions, then the stiffness matrix will be singular, reflecting the fact that BVP either does not have a solution or has infinitely many solutions. Special care must be taken to compute a meaningful solution to the resulting linear system.

#### **3.10 The method of shifting the data**

#### *3.10.1 Inhomogeneous Dirichlet conditions on a rectangle*

In a two-dimensional problem, inhomogeneous boundary conditions are handled just as in one dimension. Inhomogeneous Dirichlet conditions are addressed via the method of shifting the data (with a specially chosen piecewise linear function), while inhomogeneous Neumann conditions are taken into account directly when deriving the weak form. Both types of boundary conditions lead to a change in the load vector.

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

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

$$-\Delta u = f(\mathbf{x}), \quad \mathbf{x} \in \ \Omega = (0, a) \times (0, b), \quad u(\mathbf{x}) = g(\mathbf{x}) = \begin{cases} g\_1(\mathbf{x}), & \mathbf{x} \in \Gamma\_1, \\ g\_2(\mathbf{x}), & \mathbf{x} \in \Gamma\_2, \\ g\_3(\mathbf{x}), & \mathbf{x} \in \Gamma\_3, \\ g\_4(\mathbf{x}), & \mathbf{x} \in \Gamma\_4, \end{cases} \tag{38}$$

where Γ1, Γ2, Γ3, and Γ<sup>4</sup> are, respectively, the bottom, right, top, and left boundary edges of the rectangular domain Ω ¼ ð Þ� 0, *a* ð Þ 0, *b* . We will assume that the boundary data are continuous, so

$$\mathbf{g\_1(0)} = \mathbf{g\_4(0)}, \quad \mathbf{g\_1(a)} = \mathbf{g\_2(0)}, \quad \mathbf{g\_2(b)} = \mathbf{g\_3(a)}, \quad \mathbf{g\_3(0)} = \mathbf{g\_4(b)}.$$

Suppose we find a function *w* defined on Ω and satisfying *w*ð Þ¼ **x** *g*ð Þ **x** for all **x**∈ ∂Ω. We then define *v* ¼ *u* � *w* and note that

$$-\Delta v = -\Delta u + \Delta w = f(\mathbf{x}) + \Delta w = \hat{f}(\mathbf{x}),$$

and *v*ð Þ¼ **x** *u*ð Þ� **x** *w*ð Þ¼ **x** 0 for all **x**∈ ∂Ω. We can then solve

$$- \Delta v = \hat{f}(\mathbf{x}), \quad \mathbf{x} \in \ \Omega, \quad v(\mathbf{x}) = 0, \quad \mathbf{x} \in \partial \Omega. \tag{39}$$

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 x*ð Þ¼ , *y c*<sup>0</sup> þ *c*1*x* þ *c*2*y* þ *c*3*xy*, which assumes the desired boundary values at the corners:

$$\begin{array}{lll} q(\mathbf{0},\mathbf{0}) = \mathbf{g}\_1(\mathbf{0}) = \mathbf{g}\_4(\mathbf{0}), & q(a,\mathbf{0}) = \mathbf{g}\_1(a) = \mathbf{g}\_2(\mathbf{0}), \\ \mathbf{-g}\_3(a), & q(\mathbf{0},b) = \mathbf{g}\_3(\mathbf{0}) = \mathbf{g}\_4(b). \end{array}$$

A direct calculation shows that

$$\begin{aligned} \mathbf{c}\_{0} &= \mathbf{g}\_{1}(\mathbf{0}), \quad \mathbf{c}\_{1} = \frac{\mathbf{g}\_{1}(a) - \mathbf{g}\_{1}(\mathbf{0})}{a}, \quad \mathbf{c}\_{2} = \frac{\mathbf{g}\_{4}(b) - \mathbf{g}\_{4}(\mathbf{0})}{b}, \\\mathbf{c}\_{3} &= \frac{\mathbf{g}\_{2}(b) + \mathbf{g}\_{1}(\mathbf{0}) - \mathbf{g}\_{1}(a) - \mathbf{g}\_{4}(b)}{ab}. \end{aligned}$$

We then define

$$h(\mathbf{x}) = \begin{cases} h\_1(\mathbf{x}) = \mathbf{g}\_1(\mathbf{x}) - \left(\mathbf{g}\_1(\mathbf{0}) + \frac{\mathbf{g}\_1(a) - \mathbf{g}\_1(\mathbf{0})}{a}\mathbf{x}\right), & \mathbf{x} \in \Gamma\_1, \\ h\_2(\mathbf{y}) = \mathbf{g}\_2(\mathbf{y}) - \left(\mathbf{g}\_2(\mathbf{0}) + \frac{\mathbf{g}\_2(b) - \mathbf{g}\_2(\mathbf{0})}{b}\mathbf{y}\right), & \mathbf{x} \in \Gamma\_2, \\ h\_3(\mathbf{x}) = \mathbf{g}\_3(\mathbf{x}) - \left(\mathbf{g}\_3(\mathbf{0}) + \frac{\mathbf{g}\_3(a) - \mathbf{g}\_3(\mathbf{0})}{a}\mathbf{x}\right), & \mathbf{x} \in \Gamma\_3, \end{cases}$$

$$\begin{cases} \mathbf{u}^{(2)} \\ \begin{aligned} h\_3(\mathbf{x}) &= \mathbf{g}\_3(\mathbf{x}) - \left( \mathbf{g}\_3(\mathbf{0}) + \frac{\mathbf{g}\_3(a) - \mathbf{g}\_3(\mathbf{0})}{a} \mathbf{x} \right), \quad \mathbf{x} \in \Gamma\_3, \\\ h\_4(\mathbf{y}) &= \mathbf{g}\_4(\mathbf{y}) - \left( \mathbf{g}\_4(\mathbf{0}) + \frac{\mathbf{g}\_4(b) - \mathbf{g}\_4(\mathbf{0})}{b} \mathbf{y} \right), \quad \mathbf{x} \in \Gamma\_4. \end{aligned} \end{cases}$$

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

$$h\_1(\mathbf{0}) = h\_1(a) = h\_2(\mathbf{0}) = h\_2(b) = h\_3(\mathbf{0}) = h\_3(a) = h\_4(\mathbf{0}) = h\_4(b) = \mathbf{0}.$$

Finally, we define

$$\begin{aligned} w(\mathbf{x}, \boldsymbol{y}) &= (c\_0 + c\_1 \mathbf{x} + c\_2 \mathbf{y} + c\_3 \mathbf{x} \mathbf{y}) + \left( h\_1(\mathbf{x}) + \frac{h\_3(\mathbf{x}) - h\_1(\mathbf{x})}{b} \mathbf{y} \right), \\ &+ \left( h\_4(\mathbf{y}) + \frac{h\_2(\mathbf{y}) - h\_4(\mathbf{y})}{a} \mathbf{x} \right). \end{aligned}$$

The reader should notice how the second term interpolates between the boundary values on Γ<sup>1</sup> and Γ3, while the third term interpolates between the boundary values on Γ<sup>2</sup> and Γ4. In order for these two terms not to interfere with each other, it is necessary that the boundary data be zero at the corners. It was for this reason that we transformed the *gi* 0 s into the *hi* 0 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

$$-\Delta u = f(\mathbf{x}), \quad \mathbf{x} \in \ \Omega = (0, a) \times (0, b), \quad \mathbf{n} \cdot \nabla u(\mathbf{x}) = \mathbf{g}(\mathbf{x}) = \begin{cases} g\_1(\mathbf{x}), & \mathbf{x} \in \Gamma\_1, \\ g\_2(\mathbf{x}), & \mathbf{x} \in \Gamma\_2, \\ g\_3(\mathbf{x}), & \mathbf{x} \in \Gamma\_3, \\ g\_4(\mathbf{x}), & \mathbf{x} \in \Gamma\_4, \end{cases} \tag{40}$$

where Γ1, Γ2, Γ3, and Γ<sup>4</sup> are, respectively, the bottom, right, top, and left boundary edges of the rectangular domain Ω ¼ ð Þ� 0, *a* ð Þ 0, *b* . We first note that this is equivalent to

$$\begin{array}{cccc} -u\_{\mathcal{I}}(\mathbf{x}) = \mathbf{g}\_{1}(\mathbf{x}), & \mathbf{x} \in \Gamma\_{1}, \quad u\_{\mathbf{x}}(\mathbf{x}) = \mathbf{g}\_{2}(\mathbf{y}), & \mathbf{x} \in \Gamma\_{2}, \quad u\_{\mathcal{I}}(\mathbf{x}) = \\ -\mathbf{g}\_{3}(\mathbf{x}), & \mathbf{x} \in \Gamma\_{3}, & -u\_{\mathbf{x}}(\mathbf{x}) = \mathbf{g}\_{4}(\mathbf{y}), & \mathbf{x} \in \Gamma\_{4}. \end{array}$$

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

$$-\mathfrak{u}\_{\mathfrak{x}\mathfrak{y}}(\mathfrak{x},\mathbf{0}) = \mathfrak{g}\_{\mathbf{1}}'(\mathfrak{x}), \quad -\mathfrak{u}\_{\mathfrak{x}\mathfrak{x}}(\mathbf{0},\mathfrak{y}) = \mathfrak{g}\_{\mathcal{A}'}(\mathfrak{y}),$$

which together imply that *g*<sup>0</sup> <sup>1</sup>ð Þ¼ 0 *g*40ð Þ 0 . By similar reasoning, we have all of the following conditions:

$$\mathbf{g}'\_1(\mathbf{0}) = \mathbf{g}\_{4'}(\mathbf{0}), \quad \mathbf{g}'\_1(\mathbf{0}) = \mathbf{g}\_{4'}(\mathbf{0}), \quad -\mathbf{g}'\_1(\mathbf{a}) = \mathbf{g}'\_2(\mathbf{0}), \quad \mathbf{g}'\_2(\mathbf{b}) = \mathbf{g}'\_3(\mathbf{a}).\tag{41}$$

We will assume that (41) holds.

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

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(\mathbf{x}) = -a\mathbf{x} + \frac{a+\beta}{2a}\mathbf{x}^2 \quad \text{satisfies} \quad \psi'(\mathbf{0}) = -a, \quad \psi'(a) = \beta. \tag{42}
$$

The first step is to transform the boundary data *gl* ð Þ *x* to a function *h*1ð Þ *x* satisfying *h*10ð Þ¼ 0 *h*10ð Þ¼ *a* 0, and similarly for *g*2, *g*3, *g*<sup>4</sup> 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 x*ð Þ , *y* satisfying the conditions

$$
\mu\_{\mathbf{xy}}(\mathbf{0}, \mathbf{0}) = -\mathsf{g}\_1'(\mathbf{0}), \quad \mu\_{\mathbf{xy}}(\mathbf{a}, \mathbf{0}) = -\mathsf{g}\_1'(a), \quad \mu\_{\mathbf{xy}}(\mathbf{0}, b) = -\mathsf{g}\_3'(\mathbf{0}), \ \mathfrak{u}\_{\mathbf{xy}}(\mathbf{a}, b) = -\mathsf{g}\_2'(b).
$$

We can satisfy these conditions with a function of the form *q x*ð Þ¼ , *y c*0*xy* þ *<sup>c</sup>*1*x*<sup>2</sup>*<sup>y</sup>* <sup>þ</sup> *<sup>c</sup>*2*xy*<sup>2</sup> <sup>þ</sup> *<sup>c</sup>*3*x*<sup>2</sup>*y*2. The reader can verify that the necessary coefficients are

$$\begin{aligned} c\_0 &= -\mathbf{g}\_1'(\mathbf{0}), \quad c\_1 = \frac{\mathbf{g}\_1'(\mathbf{0}) - \mathbf{g}\_1'(a)}{2a}, \quad c\_2 = \frac{\mathbf{g}\_3'(\mathbf{0}) + \mathbf{g}\_1'(\mathbf{0})}{2b}, \\\ c\_3 &= \frac{\mathbf{g}\_2'(b) + \mathbf{g}\_1'(a) - \mathbf{g}\_3'(\mathbf{0}) - \mathbf{g}\_1'(\mathbf{0})}{4ab}. \end{aligned}$$

If *w* is to satisfy the desired Neumann conditions, then *w* � *q* ¼ *hi* on Γ*i*, *i* ¼ 1 � 4, where

$$h\_1(\mathbf{x}) = \mathbf{g}\_1(\mathbf{x}) + c\_0 \mathbf{x} + c\_1 \mathbf{x}^2, \quad h\_2(\mathbf{y}) = \mathbf{g}\_2(\mathbf{y}) - (c\_0 + 2ac\_1)\mathbf{y} - (c\_2 + 2ac\_3)\mathbf{y}^2,$$

$$h\_3(\mathbf{x}) = \mathbf{g}\_3(\mathbf{x}) - (c\_0 + 2bc\_2)\mathbf{x} - (c\_1 + 2bc\_3)\mathbf{x}^2, \quad h\_4(\mathbf{y}) = \mathbf{g}\_4(\mathbf{y}) + c\_0 \mathbf{y} + c\_2 \mathbf{y}^2.$$

We can now define *w* � *q* by the interpolation described by (42):

$$w(\mathbf{x}, \mathbf{y}) = q(\mathbf{x}, \mathbf{y}) - h\_1(\mathbf{x})\mathbf{y} + \frac{h\_3(\mathbf{x}) + h\_1(\mathbf{x})}{2\mathbf{b}}\mathbf{y}^2 - h\_4(\mathbf{y})\mathbf{x} + \frac{h\_2(\mathbf{y}) + h\_4(\mathbf{y})}{2a}\mathbf{y}\mathbf{x}^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 *λ*∈ and *u* such that

$$-\Delta u = \lambda u,\quad\text{in }\Omega,\quad\nabla u \cdot \mathbf{n} = 0,\quad\text{on }\partial\Omega.\tag{43}$$

In order to derive the weak formulation, we multiply (43) with *v*∈*V*, integrate over Ω and use Green's formula to obtain

$$
\lambda \int\_{\Omega} u v \mathrm{d}x \mathrm{d}y = -\int\_{\Omega} v \Delta u \mathrm{d}x \mathrm{d}y = \int\_{\Omega} \nabla u \cdot \nabla v \mathrm{d}x \mathrm{d}y - \int\_{\partial \Omega} v \nabla u \cdot \mathbf{n} \mathrm{d}s = \int\_{\Omega} \nabla u \cdot \nabla v \mathrm{d}x \mathrm{d}y.
$$

We obtain the weak form: Find *u*∈*V* such that

$$\int\_{\Omega} \nabla u \cdot \nabla v \mathrm{d}x dy = \lambda \int\_{\Omega} u v \mathrm{d}x dy, \quad v \in V. \tag{44}$$

The FE method in 2D is defined as follows: Find *λ<sup>h</sup>* ∈ and *uh* ∈*Vh* such that

$$
\int\_{\Omega} \nabla u\_h \cdot \nabla v\_h d\mathbf{x} dy = \lambda\_h \int\_{\Omega} u\_h v\_h d\mathbf{x} dy, \quad v\_h \in V\_h,\tag{45}
$$

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>¼</sup> <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>1</sup>*cjϕ<sup>j</sup>* into (45) and picking *vh* <sup>¼</sup> *<sup>ϕ</sup>i*, we obtain

$$\sum\_{j=1}^{N} c\_j \left( \int\_{\Omega} \nabla \phi\_j \cdot \nabla \phi\_i d\mathbf{x} dy - \lambda\_h \int\_{\Omega} \phi\_i \phi\_j d\mathbf{x} dy \right) = 0, \quad i = 1, 2, \dots, N.$$

This leads to an algebraic system of the form *A***c** ¼ *λhM***c**, *i:e:* an algebraic eigenvalue problem.

#### **3.12 Error analysis**

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

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

The weak form: Find *u*∈*V*<sup>0</sup> such that

$$\int\_{\Omega} \nabla u \cdot \nabla v \mathrm{d}x dy = \int\_{\Omega} f v \mathrm{d}x dy, \quad v \in V\_0.$$

The FE approximation is defined as follows: Find *uh* ∈*Vh*,0 such that

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

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>* � �. Expressing *uh* <sup>¼</sup> <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>1</sup>*cjϕ<sup>j</sup>* and picking *vh* ¼ *ϕi*, we obtain

$$\sum\_{j=1}^{N} c\_j \left( \int\_{\Omega} \nabla \phi\_j \cdot \nabla \phi\_i d\mathbf{x} dy \right) = \int\_{\tilde{\Omega}} f \phi\_i d\mathbf{x} dy, \quad i = 1, 2, \dots, N.$$

This leads to system of the form *<sup>A</sup>***<sup>c</sup>** <sup>¼</sup> **<sup>b</sup>**, where the entries of *<sup>A</sup>* <sup>∈</sup> *<sup>N</sup>*�*<sup>N</sup>* and **b**∈ *<sup>N</sup>* are

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

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

**Theorem 3.7 (Galerkin orthogonality)** *Let u*∈*V*<sup>0</sup> *denote the weak solution and uh* ∈*Vh*,0 *the corresponding FE method approximation. Then*

$$\int\_{\Omega} \nabla (u - u\_h) \cdot \nabla v\_h d\mathbf{x} dy = \mathbf{0}, \quad v\_h \in V\_{h,0} \dots$$

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

Now, let j j j j j j *v* <sup>2</sup> <sup>¼</sup> <sup>Ð</sup> <sup>Ω</sup>∇*<sup>v</sup>* � <sup>∇</sup>*vdxdy* <sup>¼</sup> <sup>Ð</sup> <sup>Ω</sup>j j ∇*v* 2 *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*∈*V*<sup>0</sup> *denote the weak solution and uh* ∈*Vh*,0 *the corresponding FE method approximation. Then*

$$|||u - u\_h||| \le |||u - v\_h|||, \quad v\_h \in V\_{h,0} \dots$$

**Theorem 3.9** *Let u* ∈*V*<sup>0</sup> *denote the weak solution and uh* ∈*Vh*,0 *the corresponding FE method approximation. If u* ∈*C*<sup>2</sup> ð Þ Ω , *then there exists C independent of hT and u such that*

$$|||\boldsymbol{u} - \boldsymbol{u}\_h|||\_{L^2(\Omega)}^2 \leq \mathcal{C} \sum\_{T \in \mathcal{T}\_h} h\_T^2 \left\| D^2 \boldsymbol{u} \right\|\_{L^2(T)}^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 (a\nabla u) + \mathbf{b} \cdot \nabla u + cu = f, \quad \text{in } \Omega, \quad u = 0, \quad \text{on } \partial\Omega. \tag{46}$$

We seek a weak solution *<sup>u</sup>* in *<sup>V</sup>*<sup>0</sup> <sup>¼</sup> *<sup>v</sup>*∈*V v* j j<sup>∂</sup><sup>Ω</sup> <sup>¼</sup> <sup>0</sup> � �. In order to derive the weak formulation, we multiply (46) with *v*∈*V*0, integrate over Ω and use Green's formula to obtain

$$\int\_{\Omega} f v \mathbf{d} x dy = -\int\_{\Omega} v \nabla \cdot (a \nabla u) dx dy + \int\_{\Omega} v \mathbf{b} \cdot \nabla u dx dy + \int\_{\Omega} cuv dx dy$$

$$= \int\_{\Omega} a \nabla u \cdot \nabla v dx dy - \int\_{\partial \Omega} v \nabla u \cdot \mathbf{n} ds + \int\_{\Omega} v \mathbf{b} \cdot \nabla u dx dy + \int\_{\Omega} cuv dx dy$$

$$= \int\_{\Omega} a \nabla u \cdot \nabla v dx dy + \int\_{\Omega} v \mathbf{b} \cdot \nabla u dx dy + \int\_{\Omega} cuv dx dy.$$

Note that there is no need to apply Green's formula to Ð <sup>Ω</sup>*v***b** � ∇*udxdy*. We obtain the weak form: Find *u*∈*V*<sup>0</sup> such that

$$\int\_{\Omega} a \nabla u \cdot \nabla v \mathrm{d}x dy + \int\_{\Omega} v \mathbf{b} \cdot \nabla u \mathrm{d}x dy + \int\_{\Omega} cuv d\mathbf{x} dy = \int\_{\Omega} f v \mathrm{d}x dy, \quad v \in V\_0.$$

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

$$\int\_{\Omega} a \nabla u\_h \cdot \nabla v\_h d\mathbf{x} dy + \int\_{\Omega} v\_h \mathbf{b} \cdot \nabla u\_h d\mathbf{x} dy + \int\_{\Omega} c u\_h v\_h d\mathbf{x} dy = \int\_{\Omega} f v\_h d\mathbf{x} dy, \quad v\_h \in V\_{h,0}, \tag{47}$$

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>¼</sup> <sup>P</sup>*<sup>N</sup> <sup>j</sup>*¼<sup>1</sup>*cjϕ<sup>j</sup>* into (47) and picking *vh* <sup>¼</sup> *<sup>ϕ</sup>i*, we obtain

$$\sum\_{j=1}^{N} \mathbf{c}\_{j} \left( \int\_{\Omega} \mathbf{a} \nabla \phi\_{j} \cdot \nabla \phi\_{i} \mathbf{d} \mathbf{x} dy + \int\_{\Omega} \phi\_{i} \mathbf{b} \cdot \nabla \phi\_{j} \mathbf{d} \mathbf{x} dy + \int\_{\Omega} \mathbf{c} \phi\_{i} \phi\_{j} \mathbf{d} \mathbf{x} dy \right) = \int\_{\Omega} f \phi\_{i} \mathbf{d} \mathbf{x} dy, \quad i = 1, 2, \dots, N.$$

This gives us the system ð Þ *<sup>A</sup>* <sup>þ</sup> *<sup>B</sup>* <sup>þ</sup> *<sup>C</sup>* **<sup>c</sup>** <sup>¼</sup> **<sup>d</sup>**, where **<sup>c</sup>** <sup>¼</sup> ½ � *<sup>c</sup>*1, … ,*cN <sup>t</sup>* ∈ *<sup>N</sup>* is the unknown vector and the entries of *A*, *B*,*C*∈ *<sup>N</sup>*�*<sup>N</sup>* and **d**∈ *<sup>N</sup>* are given by

$$a\_{\overrightarrow{\eta}} = \int\_{\Omega} a \nabla \phi\_j \cdot \nabla \phi\_i \mathbf{dx} dy, \quad b\_{\overrightarrow{\eta}} = \int\_{\Omega} \phi\_i \mathbf{b} \cdot \nabla \phi\_j \mathbf{dx} dy, \quad c\_{\overrightarrow{\eta}} = \int\_{\Omega} c \phi\_i \phi\_j \mathbf{dx} dy, \quad d\_i = \int\_{\Omega} f \phi\_i \mathbf{dx} dy,\tag{10.15}$$

$$c\_{\overrightarrow{\eta}} \cdot \nabla \phi\_i \cdot \nabla \phi\_j \cdot \nabla \phi\_i \cdot \nabla \phi\_j \cdot \nabla \phi\_i \cdot \nabla \phi\_j \cdot \nabla \phi\_i \cdot \nabla \phi\_j \cdot \nabla \phi\_i \cdot \nabla \phi\_j \cdot \nabla \phi\_i \cdot \nabla \phi\_j \cdot \nabla \phi\_i \cdot \nabla \phi\_j$$

for *i*, *j* ¼ 1, 2, … , *N*. Note that *B* is not symmetric, *i:e: bij* 6¼ *bji*.
