CVEM-BEM Coupling with Decoupled Orders for 2D Exterior Poisson Problems

For the solution of 2D exterior Dirichlet Poisson problems, we propose the coupling of a Curved Virtual Element Method (CVEM) with a Boundary Element Method (BEM), by using decoupled approximation orders. We provide optimal convergence error estimates, in the energy and in the weaker L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textit{L}^\text {2}$$\end{document}-norm, in which the CVEM and BEM contributions to the error are separated. This allows for taking advantage of the high order flexibility of the CVEM to retrieve an accurate discrete solution by using a low order BEM. The numerical results confirm the a priori estimates and show the effectiveness of the proposed approach.


Introduction
In this paper we deal with the following 2D problem where Ω e := R 2 \ Ω 0 is an unbounded domain, exterior to an open bounded one Ω 0 , with Lipschitz boundary Γ 0 .It is known (see [24] and the references therein) that Problem (1) admits a unique solution in the space with ω(x) := , satisfying the asymptotic conditions The constant α represents the asymptotic behaviour of u e at infinity and, here, its value is not fixed in advance.

arXiv:2203.13502v1 [math.NA] 25 Mar 2022
The above problem is of interest in many engineering and physical applications, for example when studying electric and thermal plane fields on infinite domains produced by point sources, or when solving problems of fluid flows around obstacles.Many and various numerical methods have been proposed and analysed for its solution, among which we mention the traditional BEM.This latter is the most natural way to deal with unbounded domains (for a reference, see [28] and the bibliography therein contained).Another common approach is the coupling of a classical variational or finite difference method with a transparent (absorbing or non-reflecting) condition defined on an artificial boundary Γ, properly chosen to delimit a finite computational domain.Among the most commonly used Non-Reflecting Boundary Conditions (NRBCs), those of integral type are exact (i.e.not approximated) and allow treating artificial boundaries of arbitrary, even non-convex, shapes.
The aim of this paper is to propose such a coupling by means of the interior CVEM and the one-equation BEM.Standard VEMs have been applied to a wide variety of interior problems (see the pioneering [3] for the Poisson problem and [1,26,2,6] for more recent applications), but only few papers deal with exterior problems (see [19,20,17,16] for elliptic equations).Among the CVEM approaches till now investigated, we mention those proposed in [10] and [5].
Although the latter deals with local polynomial preserving VEM spaces, we choose the former since it is well-suited for problems characterized by computational domains with prescribed curved boundaries, like ours.
The choice of using VEM, or the more general CVEM, is mainly motivated by the following reasons: it allows us to consider meshes whose elements can be of general shape, and to use local discrete spaces of arbitrarily high order by maintaining the simplicity of implementation independent of it.Moreover, the nature of the VEM allows decoupling the approximation orders and the mesh grids associated with the domain and boundary methods, without the need of using special auxiliary variables (like mortar ones) for the coupling.Indeed, by exploiting the peculiar construction of the VEM, it is possible to add hanging nodes on the edges of the elements that belong to the artificial boundary, without significantly modifying the structure of the interior mesh.
For what concerns the one-equation BEM, we recall that it has been proposed in the well known Johnson & Nédélec Coupling (JNC) (see [14,24]) and it is based on a single Boundary Integral Equation (BIE) that involves the integral operators associated to the fundamental solution (and its normal derivative) of the Laplace equation.
In the recent work [19], in which a similar problem has been studied, the authors consider the Costabel & Han Coupling (CHC) (see [15,22]) combined with an interior VEM.This approach yields to a symmetric and non-positive definite scheme but, involving a BIE of hypersingular type, turns out to be quite onerous from the computational point of view.Even if the CHC has been applied in several contexts, the JNC turns out to be very appealing from the engineering point of view, this latter being cheaper and easier to implement.We remark in addition that, unlikely in [19], we deal with the asymptotic condition (2) that entails Γ λ e (y)dΓ y = 0, where λ e (y) := ∂ue ∂n (y) denotes the normal derivative of u e along the artificial boundary Γ.As a consequence, suitable spaces satisfying identity (3) have to be considered.
For the discretization of our coupled problem we consider a full Galerkin approach based on a CVEM in the interior of the computational domain and on a BEM associated to basis functions chosen in such a way that (3) is satisfied.We study the proposed approach from the theoretical point of view in a quite general framework, and we provide optimal error estimates in the energy and in the weaker L 2 -norm.In particular, since we consider here curved domains, the use of curvilinear elements instead of polygonal ones, allows us to reach the optimal convergence rate for degrees of accuracy higher than 2, avoiding the sub-optimal rate caused by the approximation of the domain.
By a careful study, we show that the source of the approximation error of the discrete solution, both in the energy and in the L 2 -norm, can be split into two contributions: a boundary (BEM) and an interior (CVEM) one.In particular, we show that the boundary contribution behaves like h k ∂ ∂ (h ∂ denoting the maximum edge length of the artificial boundary and k ∂ representing the BEM polynomial degree of accuracy), and the interior one like h k• • (h • being the element diameter and k • the CVEM order degree).Hence, for h ∂ h • and by fixing k ∂ , it results that the bulk error dominates the boundary one up to a certain CVEM order, an aspect that allows obtaining a high accuracy of the global scheme with a low BEM order.
The paper is organized as follows: in the next section we present the model problem for the Poisson equation and its reformulation in a bounded region, obtained by introducing the artificial boundary and its associated one equation Boundary Integral Non Reflecting Boundary Condition (BI-NRBC).In Section 3 we introduce the variational formulation of the problem restricted to the finite computational domain.In Section 4 we apply the Galerkin method and we prove error estimates in the energy and in the L 2 -norm in an abstract framework, provided that suitable hypotheses are assumed.Then we show that these latter are satisfied by the CVEM-BEM approximation spaces introduced in Section 5. Finally, in the last section we detail the choice of the particular basis functions used for the approximation of the normal derivative unknown, and we present some numerical test which confirm the theoretical results.

The model problem
Let Ω e := R 2 \ Ω 0 be an unbounded domain, exterior to an open bounded domain Ω 0 ⊂ R 2 , and denote by Γ 0 := ∂Ω e its Lipschitz boundary having positive Haussdorf measure (see Figure 1 (a)).We consider the exterior Dirichlet Poisson problem (1) in the unknown solution u e , where f ∈ L 2 (Ω e ) represents a source term having a compact support in Ω e .
To determine the solution u e of Problem (1) by means of an interior domain method, we surround the physical obstacle Ω 0 by an artificial boundary Γ; this allows decomposing Ω e into a finite computational domain Ω, bounded internally by Γ 0 and externally by Γ, and an infinite residual one, denoted by Ω ∞ (see Figure 1 (b)).For the theoretical analysis of the numerical approach we propose, we need to assume that Γ 0 consists of a finite number of curves of class C m+1 , with m ≥ 0, and that Γ is a contour of class C ∞ .Denoting by u and u ∞ the restrictions of the solution u e to Ω and Ω ∞ respectively, and by n and n ∞ the unit normal vectors on Γ pointing outside Ω and Ω ∞ (consequently n ∞ = −n), we consider the following compatibility and equilibrium conditions on Γ: In the above relations and in the sequel we omit, for simplicity, the use of the trace operators to indicate the restriction of H 1 functions to the boundary Γ from the exterior or interior.
Assuming, for simplicity, that Γ is chosen such that supp(f ) is a bounded subset of Ω, the following Kirchhoff's formula allows us to represent the solution u ∞ in Ω ∞ .In (5), G and ∂G/∂n ∞,y denote, respectively, the fundamental solution of the 2D Laplace equation and its normal derivative with respect to the unit vector n ∞,y having initial point in y ∈ Γ.

The variational formulation
Let us introduce the bilinear form a : The variational formulation of Problem (7) where I stands for the identity operator and (•, •) L 2 (Ω) denotes the L 2 (Ω)-inner product.It is worth noting that, since we test Equation (7c) with µ ∈ H − 1 /2 0 (Γ), satisfying by definition µ, 1 Γ = 0, the unknown constant α does not appear in the variational formulation (8).Nevertheless, the asymptotic behaviour α is intrinsic to the interior domain problem, and it can be recovered by the numerical scheme when choosing Γ sufficiently far from the obstacle (see Example 7.2).
To reformulate the above problem in operator form, following [24], we introduce the Hilbert space

Then, we define the bilinear form
for û = (u, λ) and v = (v, µ), and the linear continuous operator Hence, we rewrite Problem (8) as follows: find û ∈ V such that whose well-posedness has been proved in [24] (see Lemma 2).
Finally, for the forthcoming analysis, it is convenient to rewrite A = B + K where the the bilinear forms B, K : V × V → R are defined as follows: In the following sections, for the solution of Problem ( 9), we will describe a numerical approach consisting of a CVEM-BEM coupling.This method and the corresponding theoretical analysis is based on that proposed for the Helmholtz problem in [16], to which we refer whenever the theoretical results therein proved hold in our context as well.
It is worth noting that the theoretical analysis for the exterior Poisson problem cannot be obtained as a particular sub-case of the Helmholtz one given in [16], by simply choosing the wave number equal to zero.Indeed, the NRBC associated to the Laplace equation is different from that of the Helmholtz one, both for what concerns the kernel functions appearing in the boundary integral operators and the choice of the discrete function spaces for the approximation of the unknown λ.In fact, in this case, since we do not know a priori the asymptotic value α in (7c), the choice of the space H − 1 /2 0 (Γ) becomes mandatory and, as a consequence, a proper discrete subspace of it is needed.Moreover, another important novelty of the theoretical study, with respect to that of [16], consists in the use of decoupled degrees of approximation for the interior CVEM and the BEM.This allows in particular the application of the CVEM with order higher than that of the BEM, a key aspect for the global scheme since the BEM requires high efforts to efficiently compute the associated system matrices.

The numerical method
To describe the numerical approach we propose to solve (9), we start by introducing a suitable decomposition of the domain Ω, which consists of generic elements and is not limited to the more commonly used triangles.
Let us denote by E a generic "polygon" having at most one curved edge and by h E its diameter; similarly we denote by e a generic "edge", eventually curved, and by h e its length.We introduce a sequence {T h• } h• of unstructured meshes T h• = {E}, which cover the domain Ω, where h • := max E∈T h• h E .We denote by T Γ h ∂ the decomposition of the artificial boundary Γ which, according to the regularity assumption required for Γ, consists of curvilinear parts.The subscript h ∂ denotes the mesh size defined by We suppose there exists a constant > 0 such that for each h • and for each element E ∈ T h• , E is star-shaped with respect to a ball of radius greater than h E and the length of any (eventually curved) edge of E is greater than h E .
For any k ∈ N, let P k (E) be the space of polynomials of degree k defined on E, and The local projector Π ∇,E k can be naturally extended to the global one being the space of piecewise polynomials with respect to the decomposition

By introducing the local bilinear form a
we can write a(u, v) Finally, we introduce the product space H 1 (T h• ) := E∈T h• H 1 (E) and define the associated broken H 1 -norm: To apply the Galerkin method to Problem (9), we introduce the discrete spaces (Γ) associated to the meshes T h• and T Γ h ∂ , respectively, and the product space where are suitable approximations of A, B and L f , respectively.Proceeding analogously as in [16], we introduce sufficient conditions on the discrete spaces, on the bilinear form B h and on the linear operator L f,h to guarantee existence and uniqueness of the solution ûh ∈ V k h and to prove convergence error estimates.
In particular we assume: for any s ≥ 1 In the above assumptions the notation where c is a positive constant that, unless explicitly stated, does not depend on any relevant parameter involved in the definition of According to the definition of the • V norm, (H1.a) and (H1.b) ensure the following approximation property for the product space V k h : Recalling that the evaluation of the bilinear form B on elements of V k h is well defined provided that a(•, •) is split into the sum of the local contributions a E (•, •), and assuming that the approximated bilinear form B h is well defined on the space H 1 (T h• ), we further assume: In the following theorem we show the validity of the inf-sup condition for the discrete bilinear form A h .
In what follows we provide the error estimate in the weaker W -norm, where W := L 2 (Ω) × H − 3 /2 (Γ).To this aim, we start by assuming the following property: (H3.b) consistency: for all q ∈ P 1 (T Such assumption, together with some of those previously introduced, allows us to prove the following approximation error estimate for the operator L f . Proof.Using (H3.b) and (H3.a), we can write From (H1.a) and by using standard polynomial approximation estimates (see, for example, Lemma 4.3.8 in [13]), we get which, combined with the previous estimate, leads to the claim.
To prove the error estimate in the weaker W -norm, we introduce the dual space W := L 2 (Ω) × H 3 /2 (Γ) and denote by [•, •] : W × W → R the associated duality pairing.Further, we define the adjoint operator A * : V → V as which, by Lemma 3 in [24], is an isomorphism, whose inverse  9), satisfies û ∈ H s+1 (Ω) × H s− 1 /2 (Γ), for s ≥ 1, the following estimate Proof.Let ŵ ∈ W and v := A * −1 ŵ the unique element that, according to the above mentioned property of and Now, choosing ẑ = û − ûh in (21), where û and ûh are the solutions of ( 9) and ( 12) respectively, and denoting by the last inequality following from the continuity of A and Lemma 4.3.By applying Theorem 4.2 and the interpolation property (13), we estimate and, since v To estimate III in (22), we add and subtract the terms B h ((Π ∇ k• u, λ), vI h ) and B((Π ∇ k• u, λ), vI h ) which, for (H2.a), are equal.Hence we get Similarly, adding and subtracting the two equal terms (see (H3.b)) Using the continuity of B h (see (H2.b)) and of B, we have . The first factor of the above inequality is estimated, by using Theorem 4.2 and standard polynomial approximations, as follows Then, using (20), we obtain Finally, the assertion easily follows combining ( 22) with ( 23), ( 24) and ( 25).

The CVEM-BEM method
In this section we describe the discrete CVEM-BEM coupling procedure for the solution of Problem (9).In particular, we will show that all the assumptions, introduced in Section 4 and used to prove Theorems 4.2 and 4.4, are satisfied.
Referring to the notation introduced in Section 4, and denoting where e 1 , . . ., e n E denote the edges of the boundary of E, whose first element e 1 is assumed to be curved and parametrized by a local map γ E : I E → e 1 , and We omit here, for brevity, the complete description of such space and we refer to [3,10] for a deeper presentation.Further, since we will use some of the theoretical results proved in [16], we also refer to this latter, in particular for what concerns the notation and the implementation details.
On the basis of the definition of the local virtual space The validity of Assumption (H1.a) is based on the proof of Lemma 5.2 in [16], in which the results hold both for the space Q k• h• and for a suitable enhanced space associated to it.For each E ∈ T h• , in the spirit of the virtual element method, we define the approximation a E h• of the bilinear form a E (see definition (11)), as follows: where s E is the standard "dofi-dofi" stabilization term (see Eq. (3.22) of [4]).The global approximate bilinear form The boundary element space X k ∂ h ∂ associated to the artificial boundary Γ is defined as follows: and we refer to [25] for the validity of the associated Assumption (H1.b).We then define the approximate bilinear form For these choices, from [16] (see Section 5.2), Assumptions (H2.a)-(H2.c)are satisfied .

Algebraic details and computational issues
In this section we briefly detail the construction of the final linear system associated with the CVEM-BEM scheme, and we give some implementation issues concerning the BEM matrices.
We start by re-ordering and splitting the complete index set S of the basis functions {Φ j } j∈S of Q k• h• as S = S I ∪ S Γ , where S I and S Γ denote the sets of indices related to the internal degrees of freedom and to those lying on Γ, respectively.Moreover, we denote by {ϕ j } j∈G the basis functions of X k ∂ h ∂ , G being the corresponding index set.In order to write the linear system associated with the discrete problem (12), we expand the unknown function Hence, using the basis functions of To write the matrix form of the above linear system, we introduce the stiffness matrix A and the matrix Q whose entries are respectively defined by In accordance with the splitting of the set of the degrees of freedom, we consider the block partitioned representation of the above matrices and vectors (with obvious meaning of the notation), and we rewrite equations ( 27) as follows: For what concerns the discretization of the BI-NRBC, by inserting ( 26) in (8b) and testing with the functions ϕ i , i ∈ G, we obtain that in matrix form reads where By combining (28) with (30) we obtain the final linear system It is worth to point out that, since in the theoretical analysis we have assumed that the boundary integral operators are not approximated, it is crucial to compute the integrals defining the BEM entries of V and K with a high accuracy.Hence, to retrieve their approximation without affecting the overall accuracy of the coupled CVEM-BEM scheme, suitable efficient quadrature formulas must be considered.In [17] we have proposed and successfully applied a smoothing technique to weaken the logsingularity of the single layer operator and to compute the corresponding entries with high accuracy by using the Gauss-Legendre product quadrature rule with few nodes.Such a strategy has been tuned for the standard nodal linear and quadratic basis functions and could be, in principle, adopted also in this context for the basis functions satisfying property (3).However, since the strategy associated to the standard Lagrangian basis is a well-established task, we take advantage of it by applying a computational trick.To describe this latter, we introduce the space , whose Lagrangian basis functions are denoted by ϕ j .
Further, we denote by V, K and Q the matrices associated with the choice of the space X k ∂ h ∂ , which differ from V, K and Q by the presence of the functions ϕ i instead of ϕ i .In the forthcoming Remark 6.1 we detail the quadrature adopted to efficiently compute V.Here we describe how to retrieve the matrices V, K and Q from the corresponding V, K and Q.To this aim it is sufficient to define the functions ϕ i as a suitable linear combination of the standard ones ϕ i and, hence, to combine accordingly the rows and/or the columns of V, K and Q.Such a combination depends on the order k ∂ , on the shape of the artificial boundary Γ and on the associated mesh T Γ h ∂ .In particular, for k ∂ = 2, we define ϕ i := c i ϕ i + ϕ i+1 , where ϕ i and ϕ i+1 are two consecutive piece-wise linear nodal basis functions.For k ∂ = 3, we distinguish the following two cases: a) ϕ 2i := c 2i ϕ 2i + ϕ 2i+1 ; b) ϕ 2i+1 := ϕ 2i+1 + c 2i+1 ϕ 2i+2 .The coefficients c are chosen such that the relation Γ ϕ = 0 is satisfied and are retrieved by applying a ν-point Gauss-Legendre quadrature formula, with ν chosen such that the integral over Γ of the associated ϕ functions is accurately computed.It is worth noting that dim(X In Figures 2 and 3, we show the basis functions of the spaces  Remark 6.1.We recall that the numerical integration difficulties in the computation of the V entries spring from the log-singularity of G(r) near the origin, the latter being the kernel of the single layer operator V. To compute the corresponding integrals with high accuracy by few nodes, we have used the very simple and efficient polynomial smoothing technique proposed in [27], referred as the q-smoothing technique.It is worth noting that such technique is applied only when the distance r approaches zero.This case corresponds to the matrix entries belonging to the main diagonal and to the co-diagonals for which the supports of the basis functions overlap or are contiguous.After having introduced the q-smoothing transformation, with q = 3, we have applied the n-point Gauss-Legendre quadrature rule with n = 9 for the outer integrals, and n = 8 for the inner ones (see [18] and Remark 3 in [17] for further details).
For the computation of all the other integrals, we have applied a 9 × 8-point Gauss-Legendre product quadrature rule.Incidentally, we point out that the integrals involving the kernel function ∂ n G, appearing in the double layer operator K, do not require a particular quadrature strategy, since its singularity 1/r is factored out by the same behaviour of ∂ n r near the origin.Hence, for the computation of the entries of the matrix K, we have directly applied a 9 × 8-point Gauss-Legendre product quadrature rule.The quadrature strategy described above guarantees the computation of all the mentioned integrals with a full precision accuracy (16-digit double precision arithmetic) for both k ∂ = 2 and k ∂ = 3.

Numerical results
In this section, we present some numerical test to validate the theoretical results and to show the effectiveness of the proposed method.
For the generation of the partitioning T h• of the computational domain Ω, we have used the GMSH software to construct unstructured conforming meshes consisting of quadrilaterals (see [21]).If a polygon E has a (straight) edge bordering with the interior boundary or with the artificial one, we transform it into a curved boundary edge by means of a suitable parametrization.We remark that, even if in principle it is possible to fully decouple the interior and boundary meshes, we consider here for simplicity the boundary mesh inherited by the interior one, for which it turns out h ∂ ≤ h • .Furthermore, we point out that in all the numerical test we have considered k ∂ = 2, 3; larger values than those considered would require a tailored quadrature technique for the accurate computation of the BEM matrices that we have not performed yet.Since this usually is considered the bottleneck of the BEM, the use of decoupled approximation orders allows us to exploit the flexibility of the CVEM to retrieve high accuracy by increasing only the approximation order k • .In Example 7.1 we will investigate this aspect.

Example 1
Let us consider the unbounded region Ω e , external to the unitary disk We consider Problem (1) with f = 0 and g(x) = x 1 + 2 prescribed on the boundary Γ 0 = ∂Ω 0 .In this case, the exact solution u(x) is known and its expression is given by u We choose as artificial boundary Γ the circumference of radius 2, so that the finite computational domain Ω is the annulus bounded internally by Γ 0 and externally by Γ.
To develop a convergence analysis, we start by considering a coarse mesh, associated to the level of refinement zero (lev.0), and all the successive refinements are obtained by halving each side of its elements.In Table 1, we report the total number of the degrees of freedom associated to the CVEM space, corresponding to each decomposition level of the computational domain, and the approximation orders k • = k ∂ = 2, 3 (see Figure 4 for the meshes corresponding to level 0 (left plot) and level 2 (right plot)).
To test our numerical approach and to validate the theoretical analysis, the order k • of the approximation spaces is chosen equal to 2 (quadratic) and 3 (cubic), and k ∂ = k • .Moreover, recalling that the approximate solution u h• is not known inside the polygons, as suggested in [10] we compute the H 1 -seminorm and L 2 -norm relative errors, and the corresponding EOC, by means of the following formulas: In the above formulas the superscript k • = 2, 3 refers to the approximation order of u, and the subscript lev refers to the refinement level.For what concerns the evaluation of these errors, to compute the associated integrals over polygons we have used the n-point quadrature formulas proposed in [29] and [30], which are exact for polynomials of degree at most 2n.For curved polygons, we have applied the generalization of these formulas suggested in [10] (see Section 4.3).In both cases, we have chosen n = 8.
In Table 2 we report ε ∇,k• lev and ε 0,k• lev and the corresponding EOC by varying the refinement level from 0 to 5. As we can see the H 1 -seminorm error and the L 2 -norm error estimates confirm the expected convergence order of the method.For this example, we further investigate the possibility of choosing different approximation orders.In particular, since the meshes we have considered to generate Table 2 have the property h • = 2h ∂ , we analyse the convergence of the scheme by fixing k ∂ and varying k • .In Figures 5 and 6 we report the behaviour of the H 1 -seminorm and L 2 -norm relative

Example 2
We consider the example proposed in [25] (and in [11]), for which Γ 0 is the boundary of the unit disk, centered at the origin of the cartesian axis, f = 0 and the datum g on Γ 0 is defined as g(x) = x 4 1 x 1 ≥ 0, 0 x 1 < 0. Solving the Dirichlet Laplace problem in polar coordinates, and expanding the solution in terms of the eigenvectors of the associated Sturm Liouville system, the solution in polar coordinates reads u(ρ, θ) =

Conclusions
We have proposed and analysed the coupling of a Curved Virtual Element Method with the one-equation Boundary Element Method to solve 2D exterior Poisson problems.The peculiarity of the scheme consists in the use of decoupled approximation orders for the interior CVEM and the boundary integral NRBC.This strategy has allowed us to exploit the well-known flexibility of the CVEM to retrieve an accurate solution by a low order approximation for the BEM.Since high order BEMs require non-trivial computational efforts to efficiently evaluate the matrix entries of the associated integral operators, the advantage of using a low order BEM turns out to be a key aspect to achieve a good accuracy and convergence rate weighted against computational costs.
The good performances obtained by applying the proposed scheme to elliptic problems, encourage us to consider it within other contexts, such as time dependent exterior problems, for which both the pure BEM and its coupling with standard interior domain methods could become prohibitive.

3
respectively, associated with the uniformly partitioned parametrization interval [0, 2π) of the particular choice of a circumference.For this choice, it is immediate to get c i = −1 for k ∂ = 2, and c 2i = c 2i+1 = −2 for k ∂ = 3.

(− 1 )
(n−1) /2 ρ −n n 5 − 20n 3 + 64n cos (nθ), from which we deduce that the asymptotic behaviour is characterized by the constant α = 3 /16 = 0.1875.We choose as artificial boundary the ellipse of semi-axes 50 and 15, so that the values of the numerical solution at the points (−50, 0) and (50, 0) can be considered good approximations of α.In Figure7we compare the behaviour of the exact and numerical solutions in the intervals [−50, −1] (left plot) and [1, 50] (right plot) for a fixed mesh of the computational domain and for different choices of the approximation orders.Besides noting a very good agreement of the solutions, we report that the corresponding absolute errors at (−50, 0) and (50, 0) are approximately equal to 4.0e − 04 for k • = 1 and k ∂ = 2, 5.0e − 05 for k • = k ∂ = 2 and 1.0e − 08 for k • = k ∂ = 3.

Table 1 :
Number of the degrees of freedom associated to the CVEM space