**3. The Scheme of the Finite Element Method**

We construct a scheme of the finite element method for finding an approximate generalized solution of the first boundary value Equation (1). We perform a triangulation of the domain Ω (see, for example, Figure 1).

We draw the curves Γ*j*, *j* = 0, ... , *n*, at distance *b j n r* , *j* = 0, ... , *n*, to the boundary *∂*Ω. Here *r* is the exponent of compression and *r* > 1; 0 < *b* < *<sup>δ</sup>*<sup>Ω</sup> <sup>2</sup> ≤ *d*, *δ*<sup>Ω</sup> is the diameter of the circle inscribed in Ω. In this case the line Γ*<sup>n</sup>* divides the domain Ω into two subdomains Ω<sup>1</sup> and Ω2. The subdomain Ω<sup>1</sup> is the outer domain on the boundary strip of width *b*, Ω<sup>2</sup> is the inner domain. On each curve Γ*j*, *j* = 0, ... , *n*, (Γ<sup>0</sup> = *∂*Ω, Γ*<sup>n</sup>* = *∂*Ω1) we fix *Mj* equidistant points, which we call the nodes. Here *Mj* = *lj*/ *b j n <sup>r</sup>* <sup>−</sup> *j*−1 *n r* <sup>+</sup> 1, *<sup>j</sup>* <sup>=</sup> 1, ... , *<sup>n</sup>*, *lj* is the length of the curve <sup>Γ</sup>*<sup>j</sup>* ([*x*] denotes the integer part of *x*) and *M*<sup>0</sup> = 2*M*1. All nodes on the curve Γ*j*, *j* = 0, ... , *n*, are connected by the broken line. Then, we connect each node on the curve Γ*j*−1, *j* = 1, ... , *n*, with closest nodes on the curve Γ*j*. As a

result, the subdomain Ω<sup>1</sup> is divided into triangles with the compression of nodes to the boundary *∂*Ω. The union of all triangles with vertices on <sup>Γ</sup>*j*−<sup>1</sup> and <sup>Γ</sup>*<sup>j</sup>* is a layer *<sup>Q</sup><sup>h</sup> <sup>j</sup>* . (In Figure 1 the subdomain Ω<sup>1</sup> is divided into the layers *Q<sup>h</sup>* <sup>1</sup>, ... , *<sup>Q</sup><sup>h</sup>* <sup>4</sup>, Γ<sup>4</sup> = *∂*Ω1). The parameter *h* denotes the greatest in length of the sides of the triangles in *Q<sup>h</sup> n*.

**Figure 1.** Triangulation of domain Ω.

The subdomain Ω<sup>2</sup> is divided quasi-uniformly into a finite number of the triangles. The sides of these triangles can not be greater than *h*. Moreover, the vertices of the triangles on the boundary *∂*Ω<sup>1</sup> belong to the set of vertices of the triangles in *∂*Ω2.

The algorithm and code description of this triangulation are given in [29].

Let <sup>Ω</sup>*<sup>h</sup>* be the union of closed triangles {*K*} = {*K*1, ... , *KNh* }, and *Kj*, *<sup>j</sup>* = 1, ... , *Nh*, is the finite element. The vertices *Gi*, *i* = 1, ... , *N*, of these triangles are the nodes of the triangulation. We denote by *N* the number of the internal nodes. To each node *Gi*, *i* = 1, ... , *N* , we assign the function *ϕi*(*x*) which is equal to 1 at the point *Gi* and zero at all other nodes, and *ϕi*(*x*) is linear on each triangle *K*. We denote by *<sup>V</sup><sup>h</sup>* the linear span {*ϕi*}*<sup>N</sup> <sup>i</sup>*=1. Next, we associate the following discrete problem with the constructed finite-dimensional space *<sup>V</sup><sup>h</sup>* <sup>⊂</sup> *<sup>W</sup>*˚ <sup>1</sup> 2,*α*(Ω): find the function *uh* ∈ *<sup>V</sup><sup>h</sup>* satisfying the equality

$$E(u\_h, w\_h) = (f, w\_h)$$

for any function *wh* ∈ *<sup>V</sup>h*.

An approximate (finite element) generalized solution will be found in the form

$$u\_h(\mathbf{x}) = \sum\_{i=1}^{N'} a\_i \varphi\_i(\mathbf{x}),$$

where *ai* = *uh*(*Gi*). We assume that *uh*(*x*) = 0, *<sup>x</sup>* ∈ <sup>Ω</sup> \ <sup>Ω</sup>*h*.

The coefficients *ai* are defined from system of equations

$$E(\mathfrak{u}\_{\hbar\prime}\mathfrak{p}\_{i}) = (f, \mathfrak{p}\_{i}), \quad i = 1, \ldots, N^{\prime}$$

or

$$Aa = F\_{\prime\prime}$$

where

$$a = (a\_1, \ldots, a\_{N'})^T, \quad F = (F\_1, \ldots, F\_{N'})^T, \quad \hat{A} = (A\_{\hat{i}\hat{j}})\_{\hat{i}}$$

*Symmetry* **2019**, *11*, 1455

$$A\_{i\bar{j}} = E(\varphi\_{i\prime}\varphi\_{\bar{j}}), \quad F\_{\bar{i}} = (f, \varphi\_{\bar{i}}), \quad i, j = 1, \dots, N'.$$

It is obvious that the approximate generalized solution of Problems (1)–(6) by the finite element method exists and is unique.

For the performed triangulation of the domain Ω with the exponent of the compression of the mesh *r* = <sup>1</sup> *<sup>λ</sup>* and for functions in the space *<sup>W</sup>*˚ <sup>2</sup> 2,*α*+*λ*−1(Ω) we have convergence estimates:

$$\begin{aligned} \left||\boldsymbol{u} - \boldsymbol{u}\_{\mathrm{h}}\right||\_{W^{1}\_{2,a}(\Omega)} &\leq \mathsf{C}\_{5} \boldsymbol{h} ||f\rho^{1+\alpha-\lambda}||\_{L\_{2}(\Omega)'} \\\\ \left||\boldsymbol{u} - \boldsymbol{u}\_{\mathrm{h}}\right||\_{L\_{2}(\Omega)} &\leq \mathsf{C}\_{6} \boldsymbol{h}^{2} ||f\rho^{1+\alpha-\lambda}||\_{L\_{2}(\Omega)}.\end{aligned}$$

Here, the positive constants *C*5, *C*<sup>6</sup> are independent of *u*, *uh*, *f* and *h*.

### **4. Numerical Experiments**

In this section we present the results of numerical experiments for two model problems.

Let Ω be a circle of unit radius with center at the point (2, 2). We consider the boundary value Equation (1) in the domain Ω. The right-hand side and coefficients of Equation (1) are given as

$$f(\mathbf{x}) = 4(1+\mu-t)\left( (\mu-2\alpha-t)\bar{\rho}^{\mu-2\alpha-t-1}(\mathbf{x})(1-\bar{\rho}(\mathbf{x})) - \bar{\rho}^{\mu-2\alpha-t}(\mathbf{x}) \right) + \bar{\rho}^{1+\mu-t}(\mathbf{x}),$$

$$a\_{11}(\mathbf{x}) = a\_{22}(\mathbf{x}) = \bar{\rho}^{-2\alpha}(\mathbf{x}), \quad a(\mathbf{x}) = 1, \quad \mathbf{a} \in \left( -\frac{1}{2}, \frac{1}{2} \right), \ \mu \in \left( \alpha, \frac{1}{2} \right), \ t < \frac{1}{2},$$

where *<sup>ρ</sup>*˜(*x*) = <sup>1</sup> − (*x*<sup>1</sup> − <sup>2</sup>)<sup>2</sup> − (*x*<sup>2</sup> − <sup>2</sup>)<sup>2</sup> and *<sup>ρ</sup>*˜(*x*) be a function that is infinitely differentiable and satisfies the following conditions:

*C*7*ρ*(*x*) ≤ *ρ*˜(*x*) ≤ *C*8*ρ*(*x*) for each point *x* of the domain Ω; 0 < *C*<sup>7</sup> ≤ *C*<sup>8</sup> < ∞. The exact solution of this problem is *u*(*x*) = *ρ*˜1−*μ*−*<sup>t</sup>* (*x*).

For finding an approximate solution of model problems we used mesh with the compression of nodes to the boundary (*Rc*), quasi-uniform mesh (*Rq*) and the finite element method scheme from paragraph three. For the mesh *Rc* we set the number of layers *n* and the exponent of compression of the mesh *r* = <sup>1</sup> *<sup>τ</sup>* , *τ* = *μ* − *α*.

We investigate the convergence rate of the approximate solution *uh* to the exact one in the norms of the spaces *L*2(Ω) and *W*<sup>1</sup> 2,*α*(Ω) on the mesh *Rc* and *Rq*. The absolute value of the error *e* = |*u* − *uh*| in the mesh nodes *Gi* on the mesh *Rc* and *Rq* is analyzed.

Model Problem 1. We set the parameters *α* = 0.01, *μ* = 0.49, *t* = 0.499, at which the solution, the coefficients and the right-hand side of the equation in Equation (1) have the form

$$u(\mathbf{x}) = \bar{\rho}^{0.991}(\mathbf{x}),$$

$$a\_{11}(\mathbf{x}) = a\_{22}(\mathbf{x}) = \bar{\rho}^{-0.02}(\mathbf{x}), \quad a(\mathbf{x}) = 1,$$

$$f(\mathbf{x}) = 0.114956 \cdot \bar{\rho}^{-1.029}(\mathbf{x}) + 3.849044 \cdot \bar{\rho}^{-0.029}(\mathbf{x}) + \bar{\rho}^{0.991}(\mathbf{x}),$$

the exponent of compression of the mesh *r* = 2.08(3).

In Table 1 we give the number of nodes and their percentage to the total number of mesh nodes *N*, in which the absolute value of the error *e* = |*u* − *uh*| is not less than the specified value of the limit error. In this Table the patterns of the absolute error distribution at the nodes of the *Rq* and *Rc* meshes are also showed. We present data on the *Rc* mesh for *N* nodes and on the *Rq* mesh for *N* and *<sup>N</sup>* <sup>2</sup> nodes.

In Table 2 we present the norms of the difference between an exact and an approximate solution *e <sup>L</sup>*2(Ω) = *u* − *uh <sup>L</sup>*2(Ω) and *e <sup>W</sup>*<sup>1</sup> 2,*α*(Ω) <sup>=</sup> *<sup>u</sup>* <sup>−</sup> *uh <sup>W</sup>*<sup>1</sup> 2,*α*(Ω) for *Rq* and *Rc* and find the ratios of the norms *β* when the mesh parameter *h* is reduced by a factor two. The value of the parameter *h* in the domain Ω<sup>2</sup> for *Rc* varies by changing number of layers *n*.


**Table 1.** Absolute value of the error *e* for Model Problem 1.

**Table 2.** The errors *e <sup>L</sup>*2(Ω) and *e <sup>W</sup>*<sup>1</sup> 2,*<sup>α</sup>*(Ω) for meshes *Rq* and *Rc* for Model Problem 1.


The distribution of the absolute values of the error *e* in the mesh nodes with a decrease in the *h* parameter by a factor of two on the meshes *Rq* and *Rc* is given in Table 3.

**Table 3.** The distribution of the error *e* on the grids *Rq* and *Rc* as *h* changes for Model Problem 1.

Model Problem 2. We set the parameters *α* = −0.49, *μ* = 0.01, *t* = 0.49, at which the solution, the coefficients and the right-hand side of Equation (1) have the form

$$u(\mathbf{x}) = \bar{\rho}^{0.52}(\mathbf{x}),$$

$$a\_{11}(\mathbf{x}) = a\_{22}(\mathbf{x}) = \bar{\rho}^{0.98}(\mathbf{x}), \quad a(\mathbf{x}) = 1,$$

$$f(\mathbf{x}) = -1.04 \cdot \bar{\rho}^{-0.5}(\mathbf{x}) + 3.12 \cdot \bar{\rho}^{0.5}(\mathbf{x}) + \bar{\rho}^{0.52}(\mathbf{x}),$$

the exponent of compression of the mesh *r* = 2.

A numerical analysis of this problem was carried out by analogy with Model Problem 1. The results of the research are presented in Tables 4–6.


**Table 4.** Absolute value of the error *e* for Model Problem 2.

**Table 5.** The error *e <sup>W</sup>*<sup>1</sup> 2,*<sup>α</sup>*(Ω) for meshes *Rq* and *Rc* for Model Problem 2.


**Table 6.** The distribution of the error *e* on the grids *Rq* and *Rc* as *h* changes for Model Problem 2.

In Figure 2a,b we present graphs of the error *e <sup>W</sup>*<sup>1</sup> 2,*α*(Ω) <sup>=</sup> *<sup>u</sup>* <sup>−</sup> *uh <sup>W</sup>*<sup>1</sup> 2,*α*(Ω) as a function of the parameter *h* on the grids *Rq* and *Rc* in a logarithmic scale. In the first case the parameter *h* decreases due to an increase in the number of layers *n* at a fixed value *b* = 1/64 (Figure 2a). In the second case *h* changes due to a decrease in the width of the border strip *b* at a fixed number *n* (Figure 2b).

**Figure 2.** The graph of the error *e <sup>W</sup>*<sup>1</sup> 2,*<sup>α</sup>*(Ω) on the grids *Rq* and *Rc* as *<sup>h</sup>* changes in a logarithmic scale for Model Problem 2. For *Rc* (**a**) *b* = *const* = 1/64, *n* is a variable; (**b**) *n* is a constant, *b* is a variable.
