INVERTED FINITE ELEMENTS APPROXIMATION
OF THE NEUMANN PROBLEM FOR SECOND
ORDER ELLIPTIC EQUATIONS IN EXTERIOR
TWO-DIMENSIONAL DOMAINS

R Belbaki, S K Bhowmik, T Z Boulmezaoud, N Kerdid, S Mziou

<figure><img src="image_1.png" /></figure>

To cite this version:

R Belbaki, S K Bhowmik, T Z Boulmezaoud, N Kerdid, S Mziou. INVERTED FINITE ELEMENTS APPROXIMATION OF THE NEUMANN PROBLEM FOR SECOND ORDER ELLIPTIC EQUATIONS IN EXTERIOR TWO-DIMENSIONAL DOMAINS. 2025. hal-04524673

HAL Id: hal-04524673

https://hal.science/hal-04524673v1

Preprint submitted on 22 Jan 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
INVERTED FINITE ELEMENTS APPROXIMATION OF THE NEUMANN
PROBLEM FOR SECOND ORDER ELLIPTIC EQUATIONS IN EXTERIOR
TWO-DIMENSIONAL DOMAINS

R. BELBAKI<fnref n="4" />, S. K. BHOWMIK<fnref n="3" />, T. Z. BOULMEZAOUD<fnref n="2" />, N. KERDID<fnref n="1" />, AND S. MZIOU<fnref n="5" />

ABSTRACT. We use inverted finite elements method for approximating solutions of second order elliptic equations with non-constant coefficients varying to infinity in the exterior of a 2D bounded obstacle, when a Neumann boundary condition is considered. After proposing an appropriate functional framework for the deployment of the method, we analyse its convergence and detail its implementation. Numerical tests performed after implementation confirm convergence and high efficiency of the method.

# 1. INTRODUCTION

In this paper we deal with the approximation of the second order elliptic equation

$$ (1) \qquad -\operatorname{div} (\sigma \nabla u) = f \quad \text{in } \mathbb{R}^2 \setminus \bar{\omega}, $$

when a Neumann boundary condition is imposed on its boundary

$$ (2) \qquad \sigma \nabla u \cdot n = g \quad \text{on } \partial \omega. $$

Here $\omega$ designates a bounded open set of $\mathbb{R}^2$ having a lispchitzian boundary $\partial\omega$, $\mathbb{R}^2\setminus\bar{\omega}$ is its exterior and $n$ is the unit normal vector pointing towards $\omega$ and outwards $\mathbb{R}^2\setminus\bar{\omega}$. Formally, we also complete equation (1) with a condition at infinity which we momentarily write as

$$ (3) \qquad u(x) = \circ(\log(|x|)^{1/2}) \quad \text{when } |x| \to +\infty. $$

The exact meaning of this condition will be clarified later. Actually, since the domain $\mathbb{R}^2\setminus\bar{\omega}$ is unbounded, appropriate weighted spaces will be used for describing the decay or the growth of functions at infinity.

Following [5], [8], [7], and [1], our approach consists to use Inverted Finite Elements Method (IFEM) introduced by one of the authors [5]. This method was used with a significant success in solving 2D Dirichlet's problem even when the coefficient $\sigma$ is indefinitely varying (see [1]). It has also been used to solve one-dimensional problems, with or without singularities (see [8] and [7]), and fully three-dimensional problems (see, e. g., [5] and [6]). This short paper is in the continuation of paper [1] in which the authors use IFEM for solving the same equation with a Dirichlet's boundary condition. Although the approach is somewhat similar, there are substantial differences between the two systems of equations as well as important technical details about the method, which fully justify this work.

The outline of this paper is as follows. In Section 2 we introduce an appropriate family of

<footnote n="4">⁴ Ecole Normale Supérieure de Kouba, Algiers, Algeria.</footnote>

<footnote n="3">³ Department of Mathematics, University of Dhaka, Bangladesh</footnote>

<footnote n="2">² Université Paris-Saclay, UVSQ, LMV, Versailles, France.</footnote>

<footnote n="1">¹ IMSIU, College of Sciences, Department of mathematics and Statistics, P. O. Box 90950, Riyadh, 11623, Saudi Arabia.</footnote>

<footnote n="5">⁵ Laboratoire Lammda, ESSTHS, Université de Sousse, Tunisia</footnote>

<footnote>Date: Wednesday 22nd January, 2025.</footnote>

<footnote>Key words and phrases. Inverted finite elements, exterior domain, Poisson equation, unbounded domains.</footnote>
weighted spaces and we prove the well posedness of the problem. A weak formulation of equation (1)-(2) is also given. In Section 3, inverted finite element method is adapted to the problem. The last section is devoted to some computational results which prove the efficiency of the method.

## 2. FUNCTION SPACES AND VARIATIONAL FORMULATION.

Our intention here is to set up and establish an adequate formulation of the system (1)-(2) that is suitable for its discretization by the inverted finite elements method. In the sequel, we set

$$ \omega_{\text{ext}} = \mathbb{R}^2 \setminus \bar{\omega}. $$

Following [15] and [4, 2], we base our approach on particular weighted Sobolev spaces. Since the 2-dimensional case is often considered as critical (see, e. g., [12]), the weight we use contains a logarithmic factor. Set

$$ \varrho(x) = (|x|^2 + 1)^{-1} (\log(2 + |x|^2))^{-2}, x \in \mathbb{R}^2, $$

and define the space $L_{\varrho}^2(\omega_{\text{ext}})$ of all the measurable functions on $\omega_{\text{ext}}$ satisfying

$$ \int_{\omega_{\text{ext}}} \varrho(x) |v|^2 dx = \int_{\omega_{\text{ext}}} \frac{|v|^2}{(|x|^2 + 1)(\log(2 + |x|^2))^2} dx < +\infty. $$

This weighted space is equipped with the natural norm

$$ \|v\|_{L_{\varrho}^2(\omega_{\text{ext}})} = \left( \int_{\omega_{\text{ext}}} \varrho(x) |v|^2 dx \right)^{1/2}. $$

Similarly, let $L_{1/\varrho}^2(\omega_{\text{ext}})$ be the space of (generalized) functions satisfying

$$ \int_{\omega_{\text{ext}}} \frac{|v|^2}{\varrho(x)} dx = \int_{\omega_{\text{ext}}} (|x|^2 + 1)(\log(2 + |x|^2))^2 |v|^2 dx < +\infty. $$

Obviously,

$$ L_{1/\varrho}^2(\omega_{\text{ext}}) \hookrightarrow L^2(\omega_{\text{ext}}) \hookrightarrow L_{\varrho}^2(\omega_{\text{ext}}). $$

For each function $\lambda \in L_{1/\varrho}^2(\omega_{\text{ext}})$, the map $\varphi_{\lambda}: v \in L_{\varrho}^2(\omega_{\text{ext}}) \mapsto \int_{\omega_{\text{ext}}} \lambda v dx$ belongs to $L_{\varrho}^2(\omega_{\text{ext}})'$, the dual of $L_{\varrho}^2(\omega_{\text{ext}})$ (thanks to Cauchy-Schwarz inequality). Thus, in view of Riesz's theorem the map $\lambda \in L_{1/\varrho}^2(\omega_{\text{ext}}) \mapsto \varphi_{\lambda} \in L_{\varrho}^2(\omega_{\text{ext}})'$ is an isometry. In other words, $L_{1/\varrho}^2(\omega_{\text{ext}})$ is a realization of the dual $L_{\varrho}^2(\omega_{\text{ext}})'$ with $L^2(\omega_{\text{ext}})$ as a pivot space. Thus, in what follows $L_{\varrho}^2(\omega_{\text{ext}})'$ will be identified to $L_{1/\varrho}^2(\omega_{\text{ext}})$.

Consider now the weighted Sobolev space $W_{\log}^1(\omega_{\text{ext}})$ composed of all the measurable functions on $\omega_{\text{ext}}$ satisfying

$$ (4) \qquad v \in L_{\varrho}^2(\omega_{\text{ext}}), \quad \nabla v \in L^2(\omega_{\text{ext}})^2, $$

and endowed with the norm

$$ \|v\|_{W_{\log}^1(\omega_{\text{ext}})} = \left( \|v\|_{L_{\varrho}^2(\omega_{\text{ext}})}^2 + \|\nabla v\|_{L^2(\omega_{\text{ext}})^2}^2 \right)^{1/2}. $$

We also consider the space $H_{1/\varrho}(\text{div}; \omega_{\text{ext}})$ of all the vector fields $v \in L^2(\omega_{\text{ext}})^2$ satisfying

$$ \varrho^{-1/2} \text{div } v \in L^2(\omega_{\text{ext}}). $$

We also set

$$ \mathcal{D}(\bar{\omega}_{\text{ext}}) = \{ u \bar{\omega}_{\text{ext}} | u \in \mathcal{D}(\mathbb{R}^2) \} $$

where $\mathcal{D}(\mathbb{R}^2)$ is the space of infinitely differentiable functions on $\mathbb{R}^2$ with a compact support. The space $\mathcal{D}(\bar{\omega}_{\text{ext}})$ is dense in $W_{\log}^1(\omega_{\text{ext}})$ and similarly $\mathcal{D}(\bar{\omega}_{\text{ext}})^2$ is dense in $H_{1/\varrho}(\text{div}; \omega_{\text{ext}})$.
It can be proven that when $|x| \to +\infty$, functions of $W_{\log}^1(\omega_{\text{ext}})$ have a decaying behaviour towards 0 or at most a logarithmic growth. More precisely, for $v \in W_{\log}^1(\omega_{\text{ext}})$ we have (see [1] for the proof):

$$ (5) \qquad \lim_{r \to +\infty} \frac{\|v(r, \cdot)\|_{L^2(\mathbb{S}^1)}}{\sqrt{\log r}} = 0, $$

where $\mathbb{S}^1$ is the unit circle of $\mathbb{R}^2$ and

$$ \|v(r, \cdot)\|_{L^2(\mathbb{S}^1)} = \left( \int_{\mathbb{S}^1} v(r, \sigma)^2 d\sigma \right)^{1/2}. $$

In other words,

$$ |v(x)| = \circ((\log|x|)^{1/2}) \text{ when } |x| \to +\infty. $$

We may also observe that functions of $W_{\log}^1(\omega_{\text{ext}})$ and $H_{1/\rho}(\text{div};\omega_{\text{ext}})$ have locally the same nature as functions of the usual Sobolev spaces $H^1$ and $H(\text{div})$ respectively (see, e. g., [11]). Actually, if $K$ stands for an arbitrary compact subset of $\bar{\omega}_{\text{ext}}$, then, for all $v \in W_{\log}^1(\omega_{\text{ext}})$, $v|_K \in H^1(K)$ where $H^1(K)$ denotes the usual Sobolev space defined over $K$. Similarly, for $v \in H_{1/\rho}(\text{div};\omega_{\text{ext}})$, $v|_K \in H(\text{div};K)$, where $H(\text{div};K) = \{v \in L^2(K)^2 \mid \text{div} \, v \in L^2(K)\}$ (see, e. g., [11]). The difference is indeed in the asymptotic behavior when $|x| \to +\infty$. In view of this, and since $\partial\omega$ is bounded, the trace operator $\gamma_0 : v \in W_{\log}^1(\omega_{\text{ext}}) \mapsto v|_{\partial\omega} \in H^{1/2}(\partial\omega)$ is well defined and continuous. Similarly, the normal trace operator

$$ \gamma_N : v \in H_{1/\rho}(\text{div}; \omega_{\text{ext}}) \mapsto v.n \in H^{-1/2}(\partial\omega) $$

is also well defined, continuous, onto and the following Green's formula holds true for $v \in H_{1/\rho}(\text{div};\omega_{\text{ext}})$ and $\varphi \in W_{\log}^1(\omega_{\text{ext}})$:

$$ (6) \qquad \int_{\omega_{\text{ext}}} v.\nabla\varphi dx = - \int_{\omega_{\text{ext}}} \text{div} \, v.\varphi dx + \langle v.n, \gamma_0\varphi \rangle_{H^{-1/2}(\partial\omega), H^{1/2}(\partial\omega)}. $$

Now, we need the following standard assumptions

($\mathcal{H}_1$) $\sigma \in L^\infty(\omega_{\text{ext}})$ and there exists a constant $\sigma_0 > 0$ such that

$$ \sigma(x) \ge \sigma_0 > 0, \text{ ae in } \omega_{\text{ext}}. $$

($\mathcal{H}_2$) $f \in L_{1/\rho}^2(\omega_{\text{ext}})$ and $g \in H^{-1/2}(\partial\omega)$ are such that

$$ (7) \qquad \int_{\omega_{\text{ext}}} f dx + \langle g, 1 \rangle_{\partial\omega} = 0, $$

and we look for solutions of (1)-(2) belonging to $W_{\log}^1(\omega_{\text{ext}})$. However, we may notice that, as in a bounded domain, but not exactly for the same reasons, uniqueness is lost for the Neumann problem (1)-(2). Actually, among polynomial functions only constant ones belong to $W_{\log}^1(\omega_{\text{ext}})$. Under similar asymptotic conditions, this loss of uniqueness occurs also in 3D situations (see, e. g., [8] and [7]), but does not appear in 3D (see, e. g., [4, 3], [5]). For this reason, we complete system (1)-(2) with the following condition

$$ (8) \qquad \int_{\omega_{\text{ext}}} \frac{u}{((|x|^2 + 1)(\log(2 + |x|^2))^2} dx = 0. $$

This condition makes sense since $\rho u \in L^1(\omega_{\text{ext}})$. It also means that $u$ is orthogonal to constant functions with respect to the inner product $(\cdot, \cdot)_{L^2(\omega_{\text{ext}})}$. We have
**Proposition 2.1.** Under assumptions $(\mathcal{H}_1)-(\mathcal{H}_2)$, system (1)-(2)-(8) has one and only one solution $u \in W_{\log}^1(\omega_{\text{ext}})$. Moreover,

$$ (9) \qquad \|u\|_{W_{\log}^1(\omega_{\text{ext}})} \le C(\|f\|_{L_{1/\theta}^2(\omega_{\text{ext}})} + \|g\|_{H^{-1/2}(\partial\omega)}), $$

for some constant $C > 0$ depending only on $\omega$ and the coefficient $\sigma$.

*Proof of Proposition 2.1.* The proof is based on the following variational formulation of the problem:

**Lemma 2.1.** Under the assumptions $(\mathcal{H}_1)-(\mathcal{H}_2)$, $u \in W_{\log}^1(\omega_{\text{ext}})$ is solution of (1)-(2)-(8) iff: for all $v \in W_{\log}^1(\omega_{\text{ext}})$

$$ (10) \qquad \int_{\omega_{\text{ext}}} \sigma \nabla u \cdot \nabla w dx + \left( \int_{\omega_{\text{ext}}} \varrho(x) u dx \right) \left( \int_{\omega_{\text{ext}}} \varrho(x) w dx \right) = \int_{\omega_{\text{ext}}} f w dx + \langle g, w \rangle_{\partial \omega}. $$

*Proof of Lemma 2.1.* If $u$ is solution of (1)-(2)-(8), then multiplying (1) by $v$ and using Green's formula (6) with $v = \nabla u$ and $\varphi = w$ gives (10). Conversely, suppose that $u$ satisfies (10). Choosing $v = 1 \in W_{\log}^1(\omega_{\text{ext}})$ and using (7) gives (8). Then, we prove that $u$ satisfies (1) and (2) in a usual manner. $\square$

Well posedness of (1)-(2)-(8) is a consequence of Lax-Milgram theorem. Actually, continuity of the bilinear form on LHS and the linear form on the RHS of (10) results from assumptions $(\mathcal{H}_1)-(\mathcal{H}_2)$. Coercivity follows from the following Hardy type inequality (see [12]): there exists a constant $C > 0$ such that for all $v \in W_{\log}^1(\omega_{\text{ext}})$

$$ (11) \qquad \inf_{\lambda \in \mathbb{R}} \int_{\omega_{\text{ext}}} \frac{|v - \lambda|^2}{(|x|^2 + 1)(\log(|x|^2 + 2))^2} dx \le C \int_{\omega_{\text{ext}}} |\nabla v|^2 dx. \qquad \square $$

### 3. INVERTED FINITE ELEMENT METHOD. DISCRETIZATION.

In this section, focus is on the discretization of the system (1)-(2)-(8) by means of inverted finite element method (IFEM), which was introduced by one of the authors in [5]. The way we unroll the method here is inspired from [1] where the Dirichlet's problem is considered. Therefore, only the main elements of the method will be introduced here. The reader can refer to [1] for more details (see also [5, 8, 7]).

The starting idea consists to introduce a domain decomposition of $\omega_{\text{ext}}$ into the form

$$ (12) \qquad \overline{\omega_{\text{ext}}} = \bar{\Omega}_0 \cup (\cup_{i=1}^M T_i), $$

where

* $\Omega_0$ is a bounded open subset surrounding $\omega$,

* $T_1, \dots, T_M$ are $M$ (big) infinite triangles satisfying the assumptions (see Figure 1):

- $T_1, \dots, T_M$ have the origin as a common fictitious vertex,

- $\tilde{T}_i \cap \Omega_0 = \emptyset$, for $1 \le i \le M$,

- For $i \ne j$, $T_i \cap T_j$ is either a common infinite edge or empty.

The term *infinite triangle* means a geometric unbounded closed set of the form

$$ T = \{\lambda_0 a_0 + \lambda_1 a_1 + \lambda_2 a_2 \mid \lambda_0 \le 0, \lambda_1 \ge 0, \lambda_2 \ge 0, \lambda_0 + \lambda_1 + \lambda_2 = 1\} $$

where $a_0, a_1$ and $a_2$ are three non-aligned points. Points $a_1$ and $a_2$ are called the real vertices of $T$, while $a_0$ is called the fictitious vertex. We set

$$ \Omega^{\infty} = \overline{T_1 \cup \dots \cup T_M}, $$
and

$$ \Omega_{\star} = \mathbb{R}^2 \setminus \overline{\Omega}_{\infty} = \overline{\omega \cup \Omega_0} $$

The domain $\Omega_{\star}$ is a fictitious *bounded* open set which will play a prominent role in the method.

In the sequel, for $1 \le i \le M$, $S_i$ designates the triangle whose vertices are 0, $a_1^{(i)}$ and $a_2^{(i)}$, where $a_1^{(i)}$ and $a_2^{(i)}$ are the real vertices of $T_i$. Thus,

$$ S_1 \cup \cdots \cup S_M = \overline{\Omega}_{\star}, \cup_{i=1}^{M} (S_i \cup T_i) = \mathbb{R}^2, $$

and

$$ \overline{\Omega}_0 \cap \overline{\Omega}_{\infty} = \overline{\Omega}_{\star} \cap \overline{\Omega}_{\infty} = \cup_{i=1}^{M} (S_i \cap T_i). $$

For each $i \le M$, let $r_i$ be the unique *affine* function defined over $S_i \cup T_i$ satisfying $r_i(0) = 0$ and $r_i(a_1^{(i)}) = r_i(a_2^{(i)}) = 1$. It can be easily seen that $r_i(x) \ge 1$ for $x \in T_i$ and $0 \le r_i(x) \le 1$ for $x \in S_i$. The *global radius* is defined on $\mathbb{R}^2$ as follows

(13) $$ r(x) = r_i(x) \text{ if } x \in S_i \cup T_i, i = 1, \dots, M. $$

This function $x \in \mathbb{R}^2 \mapsto r(x)$ is continuous on $\mathbb{R}^2$ (see [5]). Moreover, there exists two constants $c_1 > 0$ and $c_2 > 0$ such that (see [5])

(14) $$ c_1|x| \le r(x) \le c_2|x|, \text{ for all } x \in \mathbb{R}^2. $$

Define also the polygonal inversion $\Phi : \overline{\Omega}_{\infty} \mapsto \overline{\Omega}_{\star}$ by

$$ \Phi(x) = \frac{x}{r(x)^2}, \text{ for } x \in \overline{\Omega}_{\infty}. $$

The transformation $\Phi$ is bijective and continuous from $\overline{\Omega}_{\infty}$ into $\overline{\Omega}_{\star} \setminus \{0\}$ and

$$ \Phi^{-1}(x_{\star}) = \frac{x_{\star}}{r(x_{\star})^2}, \text{ for } x_{\star} \in \Omega_{\star} \setminus \{0\} $$

(although $\Phi$ can be extended to $\mathbb{R}^2 \setminus \{0\}$ and can be seen as an involution, we prefer to distinguish between $\Phi^{-1}$ and $\Phi$). We also have

$$ r(\Phi(x)) = \frac{1}{r(x)}, \text{ for } x \in \overline{\Omega}_{\infty}, $$

and

* $\Phi(x) = x \text{ iff } x \in \overline{\Omega}_{\infty} \cap \overline{\Omega}_{\star} = \overline{\Omega}_{\infty} \cap \overline{\Omega}_0$.

* $\Phi(T_i) = S_i \setminus \{0\}, i = 1, \dots, M$.

**Remark 3.1.** Note that $M$, the number of infinite triangles, is not a discretization parameter and is not intended to become large or to tend to infinity. In general, the large domains $T_1, \dots, T_M$ remain fixed during mesh refinement ((12) is indeed a domain decomposition). For example, if $\overline{\omega} \subset [-R, R]^2$ for some $R > 0$, one can choose the following decomposition (which will be adopted in the numerical tests):

(15) $$ \omega_{\text{ext}} = \Omega_0 \cup T_1 \cup T_2 \cup T_3 \cup T_4, $$

with

$$ T_1 = \{(x, y) \in \mathbb{R}^2 \mid x \ge \max(R, |y|)\}, $$

$$ T_2 = \{(x, y) \in \mathbb{R}^2 \mid y \ge \max(R, |x|)\}, $$

$$ T_3 = \{(x, y) \in \mathbb{R}^2 \mid x \le -\max(R, |y|)\}, $$

$$ T_4 = \{(x, y) \in \mathbb{R}^2 \mid y \le -\max(R, |x|)\}, $$

and $\Omega_0 = \mathbb{R}^2 \setminus (\overline{\omega} \cup T_1 \cup T_2 \cup T_3 \cup T_4)$. With such a decomposition, we have

$$ r(x, y) = \max(\frac{|x|}{R}, \frac{|y|}{R}). $$
<figure><img src="image_3.png" /><figcaption>FIGURE 1. (left) decomposition of the exterior domain $\omega_{\text{ext}}$ into the union of four infinite triangles and a bounded domain $\Omega_0$. (right) The four triangles $S_1, \dots, S_4$ make up the fictitious bounded domain $\Omega_*$.</figcaption></figure>

Now, we introduce a pair of triangulations $(\mathcal{T}, \mathcal{T}^*)$ satisfying the following requirements:

* $\mathcal{T}$ (respectively $\mathcal{T}^*$) is a conformal triangulation of $\Omega_0$ (respectively, of $\Omega_*$) which is shape regular in the classical sense,

* for all $i \le M$, the restriction of $\mathcal{T}^*$ to $S_i$ is a triangulation of $S_i$,

* the triangulations $\mathcal{T}$ and $\mathcal{T}^*$ have the same trace on $\partial\Omega^\infty \cap \partial\Omega_0 = \partial\Omega_* \cap \partial\Omega_0$.

* $\mathcal{T}^*$ is $\mu$-graded, where $\mu \in (0, 1]$ is a fixed gradation parameter (see [5]). In other words, $\mathcal{T}^*$ satisfies

$$ (16) \qquad \max_{K \in \mathcal{T}^{**}} \frac{h_K}{d_K^{1-\mu}} \le c_1^* h, $$

$$ (17) \qquad \max_{K \in \mathcal{T}^* \setminus \mathcal{T}^{**}} h_K \le c_2^* h^{1/\mu}, $$

$$ (18) \qquad \min_{K_* \in \mathcal{T}^{**}} d_{K_*} \ge c_3^* h^{1/\mu}, $$

where $\mathcal{T}^{**} = \{K_* \in \mathcal{T}^* | 0 \notin K_*\}$ and

$$ h = \max_{K \in \mathcal{T} \cup \mathcal{T}^*} h_K \text{ and } d_K = \text{dist}(0, K) := \inf_{x \in K} |x|. $$

Here $c_1^* > 0$, $c_2^* > 0$ and $c_3^*$ are three constants which are independent of the triangulation.

The construction of triangulations satisfying conditions (16), (17) and (18) is detailed in [5], [8], [7] and [1]. We may observe that triangles of $\mathcal{T}^*$ are of size $h$ when they are adjacent to the interface $\bar{\Omega}_* \cap \bar{\Omega}_\infty$, and of size $h^{1/\mu}$ when they touch the origine.

The last item in unrolling the inverted finite element method is the functional transform $v \to \hat{v}$ which associates to each function $v$ defined over $\bar{\Omega}_\infty$, the function $\hat{v}$ defined over $\bar{\Omega}_*$ as follows:

$$ (19) \qquad \hat{v}(x_*) = r(x_*)^{-\theta+1} v(\Phi^{-1}(x_*)) \text{ for } x_* \in \bar{\Omega}_*. $$

Here $\theta$ is a fixed real parameter. For reasons that will become clear later, we assume that

$$ (20) \qquad \theta > 0. $$

Conversely, $v$ can also be expressed in terms of $\hat{v}$ as follows

$$ (21) \qquad v(x) = r(x)^{-\theta+1} \hat{v}(\Phi(x)), \text{ for } x \in \bar{\Omega}_*. $$
We fix now a natural number $k \ge 1$ (which will be chosen equal in the numerical tests) and define the discrete space

$$W_h = \{v_h \in \mathcal{C}^0(\mathbb{R}^2 \setminus \omega) \mid \forall K \in \mathcal{T}, v_{h|_K} \in \mathbb{P}_k(K), \\ \forall K_* \in \mathcal{T}^*, \hat{v}_{h|_{K_*}} \in \mathbb{P}_k(K_*), \hat{v}_h(0) = 0\}.$$

Function of $W_h$ are piecewise polynomial in the FEM region $\Omega_0$, but not in the IFEM region $\Omega^\infty$. The latter fact is due to the distorsion caused by the inversion.
In addition, observe that for $v \in W_h$,

$$|\hat{v}(x_*)| \le C|x_*| \quad \text{for } x_* \in \Omega_*,$$

and in view of (21) and (14)

$$|v(x)| \le Cr(x)^{-\theta+1}r(\Phi(x)) \le C|x|^{-\theta}, \quad \text{for } x \in \Omega^\infty,$$

where $C$ is a generic constant. Similarly, we also have

$$|\nabla v(x)| \le C|x|^{-\theta-1}, \quad \text{for } x \in \Omega^\infty.$$

Since $\theta > 0$, we deduce that

$$(22) \qquad W_h \subset W_{\log}^1(\omega_{\text{ext}}).$$

Furthermore, $W_h$ is a finite dimensional space. When $k = 1$, its dimension is given by (see Section 4):

$$(23) \qquad \dim W_h = N_i + N_b + N_i^* + N_b^*.$$

where

* $N_i$ is the number of internal nodes of the triangulation $\mathcal{T}$ (that is to say the nodes which are strictly inside $\Omega_0$),

* $N_b$ is the number of nodes of the triangulation $\mathcal{T}$ belonging to the boundary $\partial\omega$,

* $N_i^*$ is the number of internal nodes in the triangulation $\mathcal{T}^*$ of $\Omega_*$ (the node 0 is not counted),

* $N_b^*$ is the number nodes which are common to the triangulations $\mathcal{T}$ and $\mathcal{T}^*$ (that is, nodes belonging to $\bar{\Omega}_* \cap \bar{\Omega}_\infty = \bar{\Omega}_0 \cap \bar{\Omega}_\infty$).

The discrete problem can be written as: find $u_h \in W_h$ such that: for all $w_h \in W_h$

$$(24) \quad \int_{\omega_{\text{ext}}} \sigma \nabla u_h \cdot \nabla w_h dx + \left( \int_{\omega_{\text{ext}}} \varrho(x) u_h dx \right) \left( \int_{\omega_{\text{ext}}} \varrho(x) w_h dx \right) = \int_{\omega_{\text{ext}}} f w_h dx + \langle g, w_h \rangle_{\partial \omega_{\text{ext}}}.$$

We state this

**Proposition 3.1.** Under assumptions $(\mathcal{H}_1)-(\mathcal{H}_2)$, the discrete problem has one and only one solution $u_h \in W_h$. Furthermore, if $u \in H_{\text{loc}}^{k+1}(\omega_{\text{ext}})$ is such that

$$(25) \qquad \forall |\lambda| \le k+1, \quad \|(\partial^\lambda u)(|x|, \sigma)\|_{L^2(\mathbb{S}^1)} \le \frac{C}{|x|^{\theta+|\lambda|}} \quad \text{for } |x| \ge R,$$

for some constants $C > 0$ and $R > 0$ are two constants, then for each $\varepsilon \in ]0, \min(1, \theta)]$, there exists a constant $C_\varepsilon > 0$, not depending on $u$ and $h$ such that

$$(26) \qquad \|u - u_h\|_{W_{\log}^1(\omega_{\text{ext}})} \le C_{\varepsilon}(h^k \|u\|_{H^{k+1}(\Omega_0)} + h^{k \min(\frac{\mu_0}{\mu}, 1)} \|u\|_{W_{k+\theta-\varepsilon}^{k+1}(\Omega^\infty)}),$$

where $\mu_0 = \frac{\theta-\varepsilon}{k}$. If in addition $\mu < \frac{\theta}{k}$, then

$$(27) \qquad \|u - u_h\|_{W_{\log}^1(\omega_{\text{ext}})} \le C^* h^k,$$

for some constant $C^*$ not depending on $h$.
Proof. Well posedness of the finite dimensional problem (24) is straightforward (this is a Galerkin approximation of an elliptic problem). Error estimate (27) follows from Céa's Lemma

$$ (28) \qquad \|u - u_h\|_{W_{\text{log}}^{1}(\omega_{\text{ext}})} \le C \inf_{v_h \in W_h} \|u - v_h\|_{W_{\text{log}}^{1}(\omega_{\text{ext}})}. $$

and estimate of the best approximation error given in [5]. The reader can refer to [1] and [5] for other details. □

## 4. IMPLEMENTATION AND COMPUTATIONAL RESULTS

Here focus is to briefly describe the implementation of the method and to show some numerical results when $k = 1$. Let us first give a basis of $W_h$. Let $(x_i)_{1\le i \in J}$ be the nodes of $\mathcal{T}$ and $(x_i)_{1\le i \in \partial J}$, with $\bar{\Omega}_* \cap \Omega^\infty$, the nodes of $\mathcal{T}$ belonging to the common common boundary (thus, $(x_i)_{1\le i \in J \setminus \partial J}$ are internal nodes of $\mathcal{T}$). Similarly, let $(x_i^*)_{i \in J^*}$ be the nodes of $\mathcal{T}^*$ from which we exclude the origine and $(x_i^*)_{i \in \partial J^*}$ nodes of $\mathcal{T}^*$ belonging to the common boundary $\bar{\Omega}_* \cap \bar{\Omega}_\infty$ (thus, $(x_i^*)_{i \in J^* \setminus \partial J^*}$ are internal nodes of $\mathcal{T}^*$). Since the triangulations $\mathcal{T}$ and $\mathcal{T}^*$ have the same trace on $\bar{\Omega}_* \cap \Omega^\infty$, the sets $\{x_i | i \in \partial J\}$ and $\{x_i^* | i \in \partial J^*\}$ are equal. Define $w_i, i \in J$, as the unique function of $W_h$ satisfying

$$ w_i(x_j) = \delta_{i,j} \text{ for all } j \in J, \quad \hat{w}_i(x_j^*) = 0 \text{ for all } j \in J^* \setminus \partial J^*. $$

Define also $w_i^*, i \in J^*$, as the unique function of $W_h$ satisfying

$$ w_i^*(x_j) = 0 \text{ for all } j \in J \setminus \partial J, \quad \hat{w}_i^*(x_j^*) = \delta_{i,j} \text{ for all } j \in J^*. $$

We may notice that

* $w_i = 0$ in $\bar{\Omega}_\infty$ for all $i \in J \setminus \partial J$.

* $w_i^* = 0$ in $\bar{\Omega}_0$ for $i \in J^* \setminus \partial J^*$.

* If $x_i \in \bar{\Omega}_* \cup \bar{\Omega}_\infty$, $i \in J$, then the support of $w_i$ lies straddles the two domains $\bar{\Omega}_0$ and $\bar{\Omega}_\infty$.

The functions $(w_i)_{i \in J}$ and $(w_i^*)_{i \in J^* \setminus \partial J^*}$ form a basis of $W_h$ and for all $v_h \in V_h$ we can write

$$ (29) \qquad v_h(x) = \sum_{i \in J} v_h(x_i) w_i(x) + \sum_{i \in J^* \setminus \partial J^*} \hat{v}_h(x_i^*) w_i^*(\Phi(x)), \quad \text{for } x \in \omega_{\text{ext}}. $$

Using this basis, one can write the discrete problem (24) into the form $AU = B$ with $A$ the stiffness matrix whose entries are integrals of the form

$$ (30) \qquad (A)_{i,j} := \int_{\omega_{\text{ext}}} \sigma(x) \nabla \psi_i \cdot \nabla \psi_j dx = \int_{\Omega_0} \sigma(x) \nabla \psi_i \cdot \nabla \psi_j dx + \int_{\Omega^\infty} \sigma(x) \nabla \psi_i \cdot \nabla \psi_j dx, $$

where $\psi_i$ and $\psi_j$ are two basis functions. The integral on $\Omega_0$ on the right hand side can be computed in the same way as in the finite element method. The second integral is over the unbounded subdomain $\Omega^\infty$ and can be transformed to a (singular) integral on the bounded domain $\Omega_*$ by means of a change of variables. In [1] it is proven that

$$ \begin{aligned} & \int_{\Omega^\infty} \sigma(x) (\nabla \psi_i)^T \nabla \psi_j \, dx = \int_{\Omega_*} (\nabla \hat{\psi}_i)^T G(x_*) (\nabla \hat{\psi}_j) \, dx_* \\ & + \int_{\Omega_*} k(x_*) \hat{\psi}_i(x_*) \hat{\psi}_j(x_*) \, dx_* \\ & + \int_{\Omega_*} V(x_*) [\hat{\psi}_i(x_*) \nabla \hat{\psi}_j(x_*) + \hat{\psi}_j(x_*) \nabla \hat{\psi}_i(x_*)]^T \, dx_* \end{aligned} $$
with

$$G(x_{\star}) = \sigma(\Phi(x_{\star})) \frac{r(x_{\star})^{2\theta-4}}{|\eta|^2} (|\eta|^2 r(x_{\star})^2 I - 2r(x_{\star}) (\eta x_{\star}^T + x_{\star}\eta^T) + 4x_{\star}x_{\star}^T),$$

$$k(x_{\star}) = -(\theta - 1)^2 \sigma(\Phi(x_{\star})) \frac{r(x_{\star})^{2(\theta-2)}}{|\eta|^2},$$

$$V(x_{\star}) = -(\theta - 1)\sigma(\Phi(x_{\star})) \frac{r(x_{\star})^{2\theta-3}}{|\eta|^2} (\eta - \frac{2x_{\star}}{r(x_{\star})}),$$

where $\eta = \eta(x_{\star})$ is the altitude vector corresponding to $T_1, \dots, T_M$ (see [1] or [5]). In the latter formulas, the vectors $x_{\star}, \nabla\psi_i$ and $\eta$ are considered as column ones.

These integrals are then computed by means of quadrature rules on triangles, as explained in [1] and [5] (see also [14, 13], [10], [16], [17] or [18]). One can also use an affine mapping to transform the integrals on a triangle $K$ to the standart triangle $\hat{K} = \{(x, y) \in \mathbb{R}^2 | x \ge 0, y \ge 0, x+y \le 1\}$, or the map $(\xi, \eta) \in [0, 1]^2 \mapsto (\xi, \eta(1-\xi)) \in \hat{K}$ which transforms the square $[0, 1]^2$ into $\hat{K}$ (allowing the use of tensorized quadrature formulas).

All the computational experiences which follow are carried out with Matlab software. The purpose of these numerical experiences is to verify that the method is efficient with this Neumann problem, when the coefficients are constant or even variable to infinity. The impact of the mesh gradation, already studied many times (see [5], [1], [8] and [7] and [6]), will not be studied again here.

In all the following simulations $\omega$ is the unit ball of $\mathbb{R}^2$ and unless otherwise indicated, we choose:

$$(31) \qquad \theta = 1.01.$$

**Example 1.** We consider the system

$$-\Delta u = f \text{ in } \mathbb{R}^2 \setminus \bar{\omega}, \quad \frac{\partial u}{\partial n} = 0 \text{ on } \partial \omega,$$

with $f$ chosen such that the exact solution is given by

$$(32) \qquad u(x, y) = \frac{x}{\sqrt{x^2 + y^2}} \sin \left( \frac{\pi}{2} \frac{1}{(x^2 + y^2)^2} \right).$$

The long distance behavior of this function is as follows:

$$(33) \qquad |u(x, y)| \sim \frac{c_1}{r^4}, \quad |\nabla u| \sim \frac{c_2}{r^5}, \quad |D^2 u| \sim \frac{c_3}{r^6},$$

with $r = \sqrt{x^2 + y^2}$. It follows that $u \in H_{\text{loc}}^2(\omega_{\text{ext}})$ and the assumption (25) is satisfied with $k = 1$. According to Proposition 3.1, one must have

$$(34) \qquad \|u - u_h\|_{W_{\text{log}}^1(\omega_{\text{ext}})} \le C_1 h,$$

regardless of the value chosen by the parameter $\mu \in (0, 1]$.

Table 1 summarizes the numerical results obtained with this example. In particular, it includes the values of the relative weighted $L^2$-error on the solution $u$ and the $L^2$-error on its gradient and this for different values of the mesh size $h$ and of the gradation parameter $\mu$ ($\mu = 1, 0.75$ or $0.5$). In Figure 2 the curves of these errors are displayed versus the mesh size $h$ in a logarithmic scale. As expected, the two errors decrease with respect to $h$ and we can observe that the value of gradation parameter $\mu$ has no significant effect on the errors (curves are quasi-superposed). In other words, since the solution decreases sufficiently at large distances, that is when $r \to +\infty$, it is not necessary to grade the inverted mesh in order to capture the behavior of the solution at infinity; the method is efficient regardless of the chosen value of the gradation mesh $\mu$. This
in accordance with forecasts of Proposition 3.1. Furthermore, by calculating the slopes of these curves we observe that the weighted $L^2$ error on the solution decreases as $h^2$ while the $L^2$ error on the gradient is of the order of $h$. Although this confirms estimate (34), there is a super-convergence of the 0th order error in the (weighted) $L^2$-norm, as in the finite element method in bounded domains. We have no theoretical proof of this observation which is specific to this example, since, contrary to the method of finite elements (see, e. g., [9]), this superconvergence of the (weighted) $L^2$-error does not occur with other problems (see [5], [1], [8], [7], and [6]) nor in the case of the 3rd example which will follow.

Figure 3 shows the 2D view of the exact solution and the approximate solutions for $\mu = 1$. We can observe that the approximate solutions is close to the exact one; the two solutions are difficult to distinguish with the naked eye.

<table><caption>TABLE 1. (Example 1) Global and local relative weighted $L^2$ errors on $u$ and relative $L^2$ error on $\nabla u$ for several values of $\mu$. Here $\Omega_0 =] -1.5, 1.5]^2 \setminus \bar{\omega}$ (FEM region) and $\Omega^\infty = \mathbb{R}^2 \setminus [-1.5, 1.5]^2$ (IFEM region).</caption><thead><tr><th>μ</th><th>Mesh size h</th><th>0.65</th><th>0.32</th><th>0.21</th><th>0.15</th><th>0.07</th></tr></thead><tbody><tr><td rowspan="5">1</td><td>Wghd. mean value</td><td>-2.3e-5</td><td>4.6e-5</td><td>8.3e-5</td><td>2.4e-5</td><td>4.7e-6</td></tr><tr><td><strong>Rel. L<sup>2</sup><sub>q</sub>-error</strong></td><td><strong>0.166</strong></td><td><strong>0.034</strong></td><td><strong>0.015</strong></td><td><strong>0.012</strong></td><td><strong>0.003</strong></td></tr><tr><td>in fem region</td><td>0.164</td><td>0.034</td><td>0.015</td><td>0.012</td><td>0.003</td></tr><tr><td>in ifem region</td><td>0.233</td><td>0.014</td><td>0.004</td><td>0.003</td><td>0.0008</td></tr><tr><td><strong>Rel. L<sup>2</sup>-error on ∇u</strong></td><td><strong>0.716</strong></td><td><strong>0.238</strong></td><td><strong>0.163</strong></td><td><strong>0.147</strong></td><td><strong>0.074</strong></td></tr><tr><td rowspan="5">0.75</td><td>in fem region</td><td>0.734</td><td>0.243</td><td>0.168</td><td>0.152</td><td>0.076</td></tr><tr><td>in ifem region</td><td>0.298</td><td>0.152</td><td>0.046</td><td>0.061</td><td>0.037</td></tr><tr><td>Wghtd. mean value</td><td>-4.4e-5</td><td>4.6e-5</td><td>8.15e-5</td><td>2.3e-5</td><td>4.6e-6</td></tr><tr><td><strong>Rel. L<sup>2</sup><sub>q</sub>-error</strong></td><td><strong>0.167</strong></td><td><strong>0.033</strong></td><td><strong>0.015</strong></td><td><strong>0.012</strong></td><td><strong>0.003</strong></td></tr><tr><td>in fem region</td><td>0.164</td><td>0.033</td><td>0.014</td><td>0.012</td><td>0.003</td></tr><tr><td rowspan="5">0.5</td><td>in ifem region</td><td>0.249</td><td>0.015</td><td>0.004</td><td>0.002</td><td>0.001</td></tr><tr><td><strong>Rel. L<sup>2</sup>-error on ∇u</strong></td><td><strong>0.716</strong></td><td><strong>0.238</strong></td><td><strong>0.163</strong></td><td><strong>0.147</strong></td><td><strong>0.074</strong></td></tr><tr><td>in fem region</td><td>0.733</td><td>0.242</td><td>0.167</td><td>0.151</td><td>0.075</td></tr><tr><td>in ifem region</td><td>0.312</td><td>0.167</td><td>0.051</td><td>0.068</td><td>0.041</td></tr><tr><td>Wght. Mean value</td><td>-6.9e-5</td><td>4.6e-5</td><td>7.9e-5</td><td>2.3e-5</td><td>4.4e-6</td></tr><tr><td rowspan="5">0.5</td><td><strong>Rel. L<sup>2</sup><sub>q</sub>-error</strong></td><td><strong>0.167</strong></td><td><strong>0.033</strong></td><td><strong>0.014</strong></td><td><strong>0.012</strong></td><td><strong>0.003</strong></td></tr><tr><td>in fem region</td><td>0.164</td><td>0.033</td><td>0.014</td><td>0.012</td><td>0.003</td></tr><tr><td>in ifem region</td><td>0.275</td><td>0.029</td><td>0.005</td><td>0.005</td><td>0.002</td></tr><tr><td><strong>Rel. L<sup>2</sup>-error on ∇u</strong></td><td><strong>0.717</strong></td><td><strong>0.240</strong></td><td><strong>0.163</strong></td><td><strong>0.148</strong></td><td><strong>0.074</strong></td></tr><tr><td>in fem region</td><td>0.733</td><td>0.242</td><td>0.167</td><td>0.151</td><td>0.075</td></tr><tr><td>in ifem region</td><td>0.345</td><td>0.201</td><td>0.064</td><td>0.085</td><td>0.052</td></tr></tbody></table>

*Example 2.* We consider the system (1)-(2) with a coefficient $\sigma$ varying to infinity:

$$ \sigma(x, y) = 1 - \frac{x^2 - y^2}{2(x^2 + y^2)}, $$

and $f$ chosen such that the exact solution is

$$ (35) \qquad u(x, y) = \frac{x}{\sqrt{x^2 + y^2}} \sin \left( \frac{\pi}{2} \frac{1}{(x^2 + y^2)^2} \right). $$

Here also the solution has a behavior of type (33) at large distances, and according to (27), the error estimate (34), remains true, independently of the value taken by $\mu$.
<figure><img src="image_4.png" /><figcaption>FIGURE 2. (example 1) global relative weighted $L^2$ error on $u$ (left figure) and global relative $L^2$ error on $\nabla u$ (right figure).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 3. (example 1) exact solution (left) and approximate solutions (right) with $\mu = 1$.</figcaption></figure>

The numerical results illustrated in Table 3 and Figure 4 are analogous to those obtained in the first example, although the coefficient $\sigma$ is variable here up to infinity, that is to say without being constant at large distances. This again confirms the usability of the method even with variable coefficients, i. e. when one cannot rely on an explicit knowledge of the fundamental solution to make an integral formulation or a series expansion in the farthest region.

## 5. CONCLUSION

As a conclusion, it has been demonstrated that the inverted finite element method easily adapts to second-order problems with Neumann boundary conditions in two-dimensional exterior domains. The numerical results it gives confirm its effectiveness and the behavior of the error is in agreement with the theoretical predictions, even when the coefficient of the equations is varying over large distances and up to infinity. Furthermore, the method provides very good results regardless of the selected value of the parameter $\mu$. Indeed, in the examples we have considered here the gradation does not affect the rapidity of the algorithm whose both relative weighted $L^2$ errors on $u$ and relative $L^2$ error on $\nabla u$ are decreasing with respect to the mesh size
<table><caption>TABLE 2. Global relative weighted $L^2$ error on $u$ (left figure) and global relative $L^2$ error on $\nabla u$ (right figure).</caption><thead><tr><th>μ</th><th>Mesh size h</th><th>0.65</th><th>0.30</th><th>0.21</th><th>0.15</th><th>0.13</th><th>0.07</th></tr></thead><tbody><tr><td rowspan="6">1</td><td>Wght. mean value</td><td>1.7e-2</td><td>1.1e-4</td><td>3.7e-4</td><td>7.8e-5</td><td>2.8e-5</td><td>1.01e-5</td></tr><tr><td><b>Rel. L<sub>q</sub><sup>2</sup>-error</b></td><td><b>0.942</b></td><td><b>0.106</b></td><td><b>0.037</b></td><td><b>0.022</b></td><td><b>0.011</b></td><td><b>0.007</b></td></tr><tr><td>in fem region</td><td>0.890</td><td>0.101</td><td>0.037</td><td>0.021</td><td>0.011</td><td>0.007</td></tr><tr><td>in ifem region</td><td>2.465</td><td>0.230</td><td>0.082</td><td>0.037</td><td>0.022</td><td>0.013</td></tr><tr><td><b>Rel. L<sup>2</sup>-error on ∇u</b></td><td><b>0.984</b></td><td><b>0.247</b></td><td><b>0.165</b></td><td><b>0.148</b></td><td><b>0.096</b></td><td><b>0.074</b></td></tr><tr><td>in fem region</td><td>0.967</td><td>0.250</td><td>0.170</td><td>0.152</td><td>0.099</td><td>0.076</td></tr><tr><td rowspan="6">0.75</td><td>Wght. mean value</td><td>1.7e-2</td><td>1.2e-4</td><td>3.7e-4</td><td>7.7e-5</td><td>2.4e-5</td><td>9.9e-6</td></tr><tr><td><b>Rel. L<sub>q</sub><sup>2</sup>-error</b></td><td><b>0.946</b></td><td><b>0.105</b></td><td><b>0.039</b></td><td><b>0.022</b></td><td><b>0.011</b></td><td><b>0.007</b></td></tr><tr><td>in fem region</td><td>0.893</td><td>0.101</td><td>0.037</td><td>0.021</td><td>0.010</td><td>0.006</td></tr><tr><td>in ifem region</td><td>2.486</td><td>0.224</td><td>0.081</td><td>0.035</td><td>0.022</td><td>0.013</td></tr><tr><td><b>Rel. L<sup>2</sup>-error on ∇u</b></td><td><b>0.984</b></td><td><b>0.247</b></td><td><b>0.165</b></td><td><b>0.148</b></td><td><b>0.096</b></td><td><b>0.074</b></td></tr><tr><td>in fem region</td><td>0.967</td><td>0.250</td><td>0.169</td><td>0.152</td><td>0.099</td><td>0.076</td></tr><tr><td rowspan="6">0.5</td><td>Wght. mean value</td><td>1.7e-2</td><td>1.2e-4</td><td>3.6e-4</td><td>3.6e-4</td><td>2.4e-5</td><td>9.8e-6</td></tr><tr><td><b>Rel. L<sub>q</sub><sup>2</sup>-error</b></td><td><b>0.948</b></td><td><b>0.105</b></td><td><b>0.038</b></td><td><b>0.022</b></td><td><b>0.011</b></td><td><b>0.007</b></td></tr><tr><td>in fem region</td><td>0.895</td><td>0.101</td><td>0.037</td><td>0.021</td><td>0.0105</td><td>0.006</td></tr><tr><td>in ifem region</td><td>2.493</td><td>0.212</td><td>0.080</td><td>0.033</td><td>0.022</td><td>0.012</td></tr><tr><td><b>Rel. L<sup>2</sup>-error on ∇u</b></td><td><b>0.985</b></td><td><b>0.249</b></td><td><b>0.165</b></td><td><b>00.149</b></td><td><b>0.096</b></td><td><b>0.075</b></td></tr><tr><td>in fem region</td><td>0.968</td><td>0.250</td><td>0.169</td><td>0.152</td><td>0.099</td><td>0.076</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>FIGURE 4. (exemple 2) global relative weighted $L^2$ error on $u$ (left figure) and global relative $L^2$ error on $\nabla u$ (right figure) versus $h$.</figcaption></figure>
$h$. However, from a qualitative point, this behavior turns out to be mainly the same as when we impose Dirichlet type conditions. Only the theoretical formulation of the problem is different.

**Acknowledgment.**

N. Kerdid, S.K. Bhowmik and S. Mziou gratefully acknowledge the support of the National Plan for Science, Technology and Information (MAARIFAH), King Abdulaziz City for Science and Technology, KSA, award number 12-MAT2996-08.

**Declarations.**

*Conflicts of interests:* The authors declare no competing interests.

**REFERENCES**

[1] S. Bhowmik, S. Mziou, R. Belbaki, and T. Boulmezaoud. Solving two dimensional second order elliptic equations in exterior domains using the inverted finite elements method. *Computers & Mathematics with Applications*, 72 (9):2315-2333, 2016.

[2] T. Z. Boulmezaoud. Vector potentials in the half-space of $\mathbb{R}^3$. *C. R. Acad. Sci. Paris Sér. I Math.*, 332(8):711-716, 2001.

[3] T. Z. Boulmezaoud. On the Stokes system and on the biharmonic equation in the half-space: an approach via weighted Sobolev spaces. *Math. Methods Appl. Sci.*, 25(5):373-398, 2002.

[4] T. Z. Boulmezaoud. On the Laplace operator and on the vector potential problems in the half-space: an approach using weighted spaces. *Math. Methods Appl. Sci.*, 26(8):633-669, 2003.

[5] T. Z. Boulmezaoud. Inverted finite elements: a new method for solving elliptic problems in unbounded domains. *M2AN Math. Model. Numer. Anal.*, 39(1):109-145, 2005.

[6] T. Z. Boulmezaoud, K. Kaliche, and N. Kerdid. Inverted finite elements for div-curl systems in the whole space. *Advances in computational mathematics*, 43:1469-1489, 2017.

[7] T. Z. Boulmezaoud, S. Mziou, B. Boudjedaa, and M. M. Babatin. Inverted finite elements for degenerate and radial elliptic problems in unbounded domains. *Jpn. J. Ind. Appl. Math.*, 32(1):237-261, 2015.

[8] T. Z. Boulmezaoud, S. Mziou, and T. Boudjedaa. Numerical approximation of second-order elliptic problems in unbounded domains. *J. Sci. Comput.*, 60(2):295-312, 2014.

[9] Ph.-G. Ciarlet. *The finite element method for elliptic problems*. North-Holland Publishing Co., Amsterdam, 1978.

[10] D. A. Dunavant. High degree efficient symmetrical Gaussian quadrature rules for the triangle. *Internat. J. Numer. Methods Engng.*, 21(6):1129-1148, 1985.

[11] V. Girault and P. A. Raviart. *Finite element methods for Navier-Stokes equations*, volume 5 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1986.

[12] J. Giroire. *Etude de quelques problèmes aux limites extérieures et résolution par équations intégrales*. Thèse de Doctorat d'Etat. Université Pierre et Marie Curie, Paris, 1987.

[13] P. C. Hammer, O. J. Marlowe, and A. H. Stroud. Numerical integration over simplexes and cones. *Math. Tables Aids Comput.*, 10:130-137, 1956.

[14] P. C. Hammer and A. H. Stroud. Numerical integration over simplexes. *Math. Tables Aids Comput.*, 10:137-139, 1956.

[15] B. Hanouzet. Espaces de Sobolev avec poids application au problème de Dirichlet dans un demi espace. *Rend. Sem. Mat. Univ. Padova*, 46:227-272, 1971.

[16] P. Silvester. Symmetric quadrature formulae for simplexes. *Math. Comp.*, 24:95-100, 1970.

[17] G. Strang and G. Fix. *An analysis of the finite element method*. Wellesley-Cambridge Press, Wellesley, MA, second edition, 2008.

[18] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu. *The finite element method: its basis and fundamentals*. Elsevier/Butterworth Heinemann, Amsterdam, seventh edition, 2013.