Lise-Marie Imbert-G Érard

Bruno Despr És
email: despres@ann.jussieu.fr

A generalized plane wave numerical method for smooth non constant coefficients

Keywords: generalized plane wave, wave equations, high order

We propose an original method based on generalized plane waves and approximated coefficients for the numerical approximation of the Helmholtz equation with a smooth constant coefficient. This is justified by a high order convergence estimate rate. Our motivation stems from the Maxwell's equations with Hermitian dielectric tensor ε which are used to model reflectometry in fusion plasma. Simplified models split them into two different propagation modes. Some numerical results are presented in dimension one and two.

Introduction

Our aim is to describe a new numerical method with generalized plane waves for the numerical approximation of time harmonic wave equations with smooth non constant coefficients. Our model problem is the Helmholtz problem with a smooth non constant coefficient

-∆ u + αu = f , x ∈ Ω , (∂ ν + iγ) u = Q (-∂ ν + iγ) u + g, x ∈ Γ .
(1.1)

Here the smooth coefficient is real, α ∈ R, but everything could be adapted replacing the homogeneous problem by the adjoint problem. The γ function can be a variable physical parameter satisfying 0 < γ m γ γ M , but for the sake of simplicity we will consider it constant and positive. The sign of α may change even if we restrict the presentation to real coefficients for simplicity. The method can be used for complex valued coefficients as well. The unknown u(x) ∈ C is sought in the space of complex valued functions.

Physical motivations

Our motivation comes from the need of efficient numerical methods for certain Maxwell's harmonic equations appearing in plasma physics. These equations read curl (curlE) -ω 2 c 2 ε(x)E = 0, x = (x, y, z), (1.2)

where E denotes the electric field, ω is the pulsation, c the sound speed and ε the dielectric tensor. The dielectric tensor represents the electromagnetic behavior of the media. The cold plasma theory Swanson (27) yields the already simplified dielectric tensor is

ε(x) =   1 -a(x) iba(x) 0 -iba(x) 1 -a(x) 0 0 0 1 -ca(x)   , i 2 = -1,
where b < 1 and c = 1b 2 . Typically the coefficient a satisfies a = x/x 0 + ã, where the perturbation is described by ã = ã0 exp -

(x-x f ) 2 ω 2 x exp - (y-y f ) 2 ω 2
y cos(q(yy f )), see Gusakov et al. (12). This is completed with boundary conditions of metallic or absorbing type. We refer to Monk (23) for the general theory of Maxwell's equations and to Cessenat & Després (2); Hiptmair et al. (14); Huttunen et al. (16) for the use of specific plane wave methods for the numerical approximation of the solutions of such problems. Two models for different propagation modes are often considered. Both are obtained from equation (1.2) under convenient assumptions on the direction and polarization of the electric field.

The two dimensional equation for what is called the O-mode reduces to

-∆ E z - ω 2 c 2 ε z (x, y)E z = 0, ∆ = ∂ xx + ∂ yy , (1.3)
on the domain Ω and can be completed by the following boundary condition

(∂ ν + iγ)E z = Q(-∂ ν + iγ)E z + g
on the boundary domain Γ . Here ∂ ν denotes the normal derivative, γ > 0 is a smooth positive function and g is for instance an L 2 -function on the boundary. Q is a smooth function allowing to fit the condition : if Q = -1 it gives a Dirichlet condition, if Q = 1 a Neumann condition or if Q = 0 a Robin condition. This O-mode (named for Ordinary mode) presents one cutoff : when ε z is negative or positive the nature of the equation (1.3) is either elliptic coercive or elliptic propagative. This coefficient ε z ∈ R is a real continuous function. It depends on the local density of electrons and on the exterior frozen magnetic field. Since the electron density is continuous, it explains why the coefficient of the equation is also a continuous function. A more general setting of the physical problem could be to introduce some dissipation with a complex valued coefficient. This case is not considered hereafter. A further simplified one dimensional model reads

- d 2 dx 2 E z + xE z = 0.
(1.4)

Equations (1.3) and (1.4) are particular cases of our model problem (1.1). The fundamental solutions are the two Airy functions Ai and Bi. The first Airy function Ai displays important properties which are fundamentally related to the physics of the problem. It will be used for validation in our numerical tests.

More generally a challenge is to adapt advanced numerical methods which are at the frontier of what is used in classical engineering so as to obtain efficient algorithms which can be used in the context of the numerical modeling of Fusion plasmas. We are in particular interested by reflectometry which is a diagnostic method to measure density in fusion plasmas, based on the reflection of the probing wave at a plasma cut-off. The local density at the cut-off layer can be deduced from the reflected signal. Therefore the mathematical model of the cut-off is crucial for the application. Our motivation is that the equation (1.3) models the cut-off as the zero of the coefficient ε z , it is then natural to take into account the smoothness feature of the coefficient in the design of the numerical method. In fact, classical numerical methods described in paragraph 1.2 consider only piecewise coefficients which would damage the description of the cut-off. We hope this work can be considered as a first step in this direction.

Plane wave methods

The numerical method that we propose is an extension of plane waves methods, such as the ultra weak variational formulation (UWVF) Després (4); Cessenat & Després (2,3); Gittelson et al. (10); Huttunen et al. (17), to problems with smooth non constant coefficients. Indeed the standard UWVF uses constant coefficients per cell. This is optimal when the physical domain can be split into sub-domains in which the coefficients are constant. But if the coefficients of the problem to solve are non constant and smooth, such a procedure introduces a priori an important error. Our aim is to propose and analyze an extension of UWVF which uses original basis functions based on the generalized plane waves Melenk (20).

We think that the approach proposed in this work is not restricted to UWVF, and can be generalized to different plane wave methods that we describe here. PUFEM Melenk (19); Melenk & Babuska (21) falls in the same class of method Perrey-Debain et al. ( 24); Pluymers et al. (25). It has also been shown that UWVF can be interpreted as a special Discontinuous Galerkin procedure Gittelson et al. (10); Hiptmair et al. (14);Farhat et al. (6,8). The analysis of the classical ultra weak variational formulation method described as a discontinuous Galerkin method is performed in Huttunen et al. (16), as well as the corresponding h and p convergence theory. It has been proved that the analysis of h-convergence takes great advantage of this fact in Buffa & Monk (1); Gittelson et al. (10). The analysis of p convergence is treated in Hiptmair et al. (13). Comparisons between these methods is investigated in Gabard et al. (9); Huttunen et al. (15); Strouboulis et al. (26); Wang et al. (28). Analysis with respect to the wavenumber k is performed in Melenk & Sauter (22). The new family of generalized plane waves described in this work generates a high order method with respect to the basis functions and the coefficients of the problem: in this direction we refer also to the Enrichment method Kalashnikova et al. (18); Farhat et al. (7) which could provide an alternative to our method in order to enrich a more conventional polynomial basis with our generalized plane wave basis. Ultimately there is no opposition between all these approaches in the sense that it is possible at the level of principles to mix polynomial basis functions and generalized plane waves to obtain a method with improved approximation properties. However it is far beyond the scope of this paper and will not be considered.

The performance of the new method relies on a further investigation of the ultra weak variational formulations method presented in Després (4). In fact, the ultra weak variational formulation uses basis function adapted to the original problem, in the sense that they are exact solution to the problem. It makes the approximation of the exact solution more relevant for a given number of elements in the mesh and a given number of basis functions per element than the approximation obtained with classical finite elements methods. Then considering smooth non constant instead of piecewise constant coefficients, it's coherent to look for more general basis functions since classical plane waves can no more be exact solutions of the problem. So the idea is to construct, in the vicinity V 0 of a point x 0 ∈ Ω basis functions

• that are generalized plane waves, say ϕ = e P(x) for x ∈ V 0 and P ∈ C[X],

• that are solution to a modified problem : -∆ ϕ + αϕ = 0, such that α satisfies the approximation property αα L ∞ (V 0 ) Ch q , where h denotes the size of V 0 , C is a constant and q 1 is a given entire number.

The aim of this paper is first to design the new adapted basis functions, second to construct the adapted tools to fit with the frame of non confirming finite element methods in order to follow the steps of the second Strang lemma for the estimation of the convergence rate and third to illustrate with basic numerical tests for the Airy equation.

Plan

This work is organized as follows. We present a family of generalized plane waves in the section 2 : these functions have been designed to be the standard plane waves in case the coefficients of the problem are constant in space and non positive. In section 3 we present the general principle of UWVF and adapt it to smooth coefficients. The next section 4 is devoted to the numerical analysis of the method. Our main theoretical result is a proof of convergence in dimension one, using the second Strang's lemma and some uniform coercivity estimates. This is probably the most original theoretical result in our work.

To our knowledge it is the first time that it is introduced and analyzed in the context of generalized plane wave methods. Numerical results are provided in section 5 to illustrate the theoretical results.

In particular we display experimental convergence estimates in dimension two. The numerical results suggest that a different normalization of the generalized plane waves may increase the accuracy, which is indeed what is observed. Additional technical material is provided in the appendix.

Generalized plane waves

Unlike the classical variational formulation used for instance by finite element methods, here the variational formulation requires meshing the domain as a preliminary task. This feature is shared by Ultra Weak Variational Formulations, Discontinuous Galerkin Methods, Enriched methods and other plane wave methods. The coupling strategy between the cells differs of course. But up to this fact the generalized plane waves can be used in principle for all such algorithms. We begin with some notations. The mesh of the domain Ω is denoted

T h = {Ω k } k∈[[1,N h ]]
, such that :

Ω = ∪Ω k , Ω k ∩ Ω j = / 0, ∀k = j, Γ k = Ω k ∩ Γ Σ k j = Ω k ∩ Ω j , oriented from Ω k to Ω j , ∂ Ω k = (∪ j Σ k j ) ∪ Γ k . Ω k Ω j Σ k j Γ j FIG.
1: Example of a meshed square domain Ω , with elements Ω k , edges Σ k j and Γ j respectively oriented toward Ω j and the exterior of the domain.

A general method

Shape or basis functions ϕ are sought as solutions of the homogeneous equation

(-∆ + α)ϕ = 0. (2.1)
If the coefficient α is constant in the cell and negative, it is sufficient to use plane waves, that is in dimension two

x = (x 1 , x 2 ) ϕ(x, y) = e √ α(d 1 x+d 2 y) with d = (d 1 , d 2 ) and (d, d) = 1. (2.2)
If the vector d is real, it is simply the direction of the plane wave. This is the basic idea of all plane wave methods. However if α is non constant in the cell, then we do not know of any simple and general analytical formula for ϕ. For example if α = x is linear with respect to the first variable, it is possible to construct ϕ from the Airy functions Ai and Bi. But the Airy functions are highly transcendental, they are not that evident to manipulate. Our main goal is to describe a method of approximation which can be used for any function α. By comparison with (2.2) it is natural to to generalize the plane formula and to consider generalized plane waves such as ϕ(x, y) = e P(x,y) P a polynomial.

One gets that ϕ is solution of (2.1) if and only if

∂ 2 ∂ x 2 P + ∂ ∂ x P 2 + ∂ 2 ∂ y 2 P + ∂ ∂ y P 2 = α(x, y).
(2.3) However many tries showed that such a representation is not sufficient. The explanation is simple: as P can be expanded as a finite series of monomial x n y m , the result is a finite series of term

∂ 2 ∂ x 2 x n y m + ∂ ∂ x x n y m 2 + ∂ 2 ∂ y 2 x n y m + ∂ ∂ y x n y m 2 = n(n -1)x n-2 + nx 2n-2 y m + x n m(m -1)y m-2 + m 2 y 2m-2 .
(2.4)

For example let us consider the case α(x, y) = x and let us look for a polynomial

P = ∑ n K ∑ m L
a nm x n y m , a KL = 0, solution of (2.3). If K 2 of L 2, the maximal degree of (2.4) cannot decrease which is contradictory with the fact that α = x is a polynomial of degree one. So K 1 and L 1. In this case the degree of (2.4) with respect to x is 0, and the same for the degree with respect to y. In summary solutions of the functional equation (2.3) cannot be polynomial in the general case. Therefore a modification is needed. Instead of considering that α is given and looking for solution of (2.3), we look for approximate solutions. That is we consider the approximate equation

-∆ + α l k ϕ l k = 0 in Ω k
where α l k is an approximation of α in Ω k .

Definition 2.1. Generalized plane waves A generalized plane wave will be understood as any function (with support in cell Ω k ) of the form ϕ l k = e P l k (x,y)

where P l k is a polynomial solution of

∂ 2 ∂ x 2 P l k + ∂ ∂ x P l k 2 + ∂ 2 ∂ y 2 P l k + ∂ ∂ y P l k 2 = α l k (x, y) (2.5)
and α l k is an approximation of α.

Design of the basis functions in dimension one

The one dimensional case is a first step to construct the coefficients α l k and the generalized plane wave functions ϕ l k . Therefore we will suppose in this section that

Ω =]a, b[⊂ R and that Ω = ∪ k∈[[1,N h ]] [x k , x k+1 ], with x k < x k+1 . The middle of the open interval Ω h =]x k , x k+1 [ is denoted by x k+1/2 = x k +x k+1 2
. Apart from providing the technical details of the construction of the basis functions, the central result of this section is an explanation why it is necessary to use different approximations α l k of the function α in the same cell [x k , x k+1 ] in order to avoid a singularity in the construction.

2.2.1 Design principle. We want here to set our choice of basis functions : in order to generalize plane wave methods, we will consider exponential of polynomials ϕ(x) = e P(x) .

Notice that we only need two basis functions per element of the mesh in dimension one: this is a common property of plane wave methods in dimension one; the reason is the number of elementary solutions of a second order differential equation which is two. Plugging the previous representation formula into the homogeneous equation -ϕ ′′ + αϕ = 0 we find the functional equation

P ′′ (x) + P ′ (x) 2 = α(x), x ∈ [x k , x k+1 ].
This equation is non linear and no simple solution is available for general right hand side α. However if α is locally constant, that is

α(x) = α(x k+1/2 ) ∈ R, x ∈ [x k , x k+1 ], then P ± k (x) = ± α(x k+1/2 )x
are two natural solutions which correspond to the two local plane waves ϕ ± k (x) = e P ± k (x) in the case α(x k+1/2 ) < 0.

Local approximation.

To ensure the local approximation of the α, one has to fit the polynomials' coefficients to approximate the Taylor expansion of the equation's coefficient α, which is performed with respect to the parameter h which represents the length of the mesh

h = max k (x k+1 -x k ).
A first idea is to approximate the Taylor expansion of α α = α ± + O(h q )

(2.6) and to look for polynomials solutions of

P ′′ ± + (P ′ ± ) 2 = α ± , x ∈ [x k , x k+1 ].
Without restriction we assume that α admits a local infinite expansion

α = ∞ ∑ i=0 d i α dx i (x k+1/2 ) x -x k+ 1 2 i , x ∈ [x k , x k+1 ].
Using the finite expansion P ± = ∑ i I β ± i y i where y = xx k+ 1 2 , one obtains

α ± = P ′′ ± + (P ′ ± ) 2 = ∑ i I β ± i y i ′′ + ∑ i I β ± i y i ′ 2
.

In order to satisfy (2.6) we have to chose I ∈ N and (β i ) 0 i I such that

∑ i I β ± i y i ′′ + ∑ i I β ± i y i ′ 2 = q-1 ∑ i=0 d i α dx i (x k+1/2 ) y i + O(h q ).
(2.7)

Identifying the coefficients in the polynomial part of the previous equation leads to a system of q equations with I unknowns. Then choosing I high enough ensures that the system is easy to solve. At the same time it is reasonable to choose I as small as possible to minimize the amount of computations. The main question is therefore to determine the optimal value of the degree of the polynomials, parameter I, with respect to the order of approximation, parameter q. Some remarks and examples follow.

• Normalization : β 0 = 0. It is always possible to take β 0 = 0 since β 0 does not show up in (2.7). It implies that the amplitude of the corresponding basis function is normalized in the cell since

e P ± x k+ 1 2 = e 0 = 1.
• Trivial case : q = I = 1. From (2.7) one obtains the equation

β 2 1 = α x k+ 1 2 . One recovers from this procedure β 1 = ± α x k+ 1 2 so P ± (x) = ± α x k+ 1 2 x -x k+ 1 2 .
In the case where α x k+ 1 2 < 0, it yields two plane waves with opposite directions. This case is the trivial one.

• Counter-example : q = I = 2. The discrete equations are obtained from the first two terms in (2.7)

   2β 2 + β 2 1 = α x k+ 1 2 ≡ a, 4β 1 β 2 = α ′ x k+ 1 2 ≡ b. (2.8) Elimination of β 2 yields -2β 3 1 + 2aβ 1 = b. (2.9)
It is of course possible in principle to compute β 1 as any root of this polynomial, β 2 will then be computed as a ratio, i.e.

β 2 = b 4β 1
. So in principle this method has the ability to generate at least two different polynomials P ± . However there is a possibility for β 1 to vanish for some value of a and b. In such a case β 2 would be singular. It must be noticed that we have used such a method in our first numerical tests: indeed it revealed a singularity near α(x) ≈ 0. Another problem is the generalization to high order : indeed this procedure requires to compute exactly the roots of a high order polynomial which generalizes (2.9); this is not possible for orders 5. This is why we do not use this method to compute the coefficients β 1 and β 2 .

• Example : q = 2 and I = 3. Since one needs at least one more degree of freedom in the system to be solved we modify (2.8) and take into account β 3 . The system becomes

2β 2 + β 2 1 = a, 6β 3 + 4β 1 β 2 = b.
(2.10)

This system has 3 unknowns and 2 equations. So it has a priori an infinite number of solutions.

Very fortunately a natural normalization condition arises, by considering that the two basis function should be linearly independent. To insure this we impose that d dx e

P + x k+ 1 2 = 0 ⇐⇒ P ′ + x k+ 1 2 = 0 and d dx e P + x k+ 1 2 = 1 ⇐⇒ P ′ + x k+ 1 2 = 1.
The first case corresponds to β 1 = 0 and the second one to β 1 = 1. With this second normalization it is evident that β 2 and β 3 can be computed explicitly from (2.10) and that the resulting formulas are just polynomial expressions with respect to all coefficients. One obtains two sets of coefficients which are

β + 1 = 1, β + 2 = 1 2 (a -1), β + 3 = 1 3 (b -2a + 2) and β - 1 = 0, β - 2 = a 2 , β - 3 = 1 3 (b -2a). Notice that α + = α - since β + 1 = β - 1 .
We use the method described in the last example at any order. The first thing is to chose a convenient degree I for any order q. In order to obtain an invertible system, one has to consider the first q terms in the left hand side of (2.7). As long as one considers the terms of degree less than or equal to I -2, the index of the coefficient from P ′′ is higher than the indexes of all the coefficients arising from (P ′ ) 2 terms. So that -as long as I -2 q -1 -the computation of the I -1 coefficients {β j } 2 j I is straightforward, given the coefficient β 1 . Moreover if I < q + 1 the terms of degree higher than I -2 in (2.7) will give an overdetermined system. For this reason the choice of Ps degree is set to be I = q + 1.

In fact, in this case we solve the system of q equations with q + 1 unknowns obtained identifying the first q coefficients in both parts of the expansion (2.7) with the normalization

β + 1 = 0 which corresponds to P ′ + x k+ 1 2 = 0 and β - 1 = 1 which corresponds to P ′ -x k+ 1 2 = 1.
The coefficients β + 2 , β + 3 , β + 4 , . . . , and β - 2 , β - 3 , β - 4 , . . . , are calculated one after the other using the formula deduced from the definition (2.7) for all n q -1

(n + 2)(n + 1)β ± n+2 = d n dx n α x k+ 1 2 - j+ j ′ =n ∑ 0 j, j ′ I-1 ( j + 1)( j ′ + 1)β ± j+1 β ± j ′ +1 .
By construction ϕ + = e P + and ϕ -= e P -are linearly independent functions. Once the polynomials P + and P -have been constructed up to order q, we set

α + = P ′′ + + (P ′ + ) 2 and α -= P ′′ -+ (P ′ -) 2 .
(2.11)

By construction the first q coefficients of these polynomials coincide. But of course all other coefficients have no reason to be equal, so α + = α -in the general case.

One can summarize as follows.

Lemma 2.1. The functions α ± defined in (2.11) satisfy the following statements.

1. They are bounded independently from the cell number k, as well as all there derivatives.

2. If α is constant in the cell, then α ± = α, and the basis functions are classical plane waves.

3. By construction there exists a constant C q such that

α ± -α L ∞ (Ω k ) C q h q
where h = |Ω k | and q is the order of the approximation.

This construction is valid even if the sign of α changes.

The last point is critical to be able to address the numerical approximation of the Airy equation.

Remark 2.1. It is also possible to choose another normalization such as

β ± 1 = ± √ x k+1/2
. This choice will be illustrated as a numerical example in section 5.

Remark 2.2. Property 3 of lemma (2.1) establishes a property of approximation with respect to h. One of our numerical tests shows that a similar property of convergence holds with respect to the order parameter q. In practice q-convergence is a highly desirable property since it allows to use big cells.

Design of the basis functions in dimension two

The fundamental equation that defines a generalized plane wave is (2.5). Description of all solutions of this equation and of effective procedures for the computation of such solutions seems difficult. Nevertheless one can rely on a simple linearization procedure in order to define a set of generalized plane waves with an order of approximation q > 1.

Linear coefficients+rotation.

A first remark is that a simple procedure exists in the case of a linear coefficient α = a + bx + cy.

Up to a local rotation it is always possible to assume that c = 0. Assuming the local form

P(x, y) = p(x) + λ y
one ends up with the equation

p ′′ (x) + p ′ (x) 2 = α(x) -λ 2
for which the procedure described in the previous section is well adapted for the construction of a discrete space of approximation. Some details about the choice of λ will be provided in the numerical section 5.3. Our numerical tests (in the numerical section) show this procedure yields a high order method at any order.

2.3.2 Linearization. For a smooth coefficient α which is not necessarily linear, it is always possible to approximate it by a linear function, that is

α(x, y) = a + b(x -x G ) + c(y -y G ) + O(h 2 ) where a = α(x G , y G ), b = ∂ x α(x G , y G ), c = ∂ y α(x G , y G ) and |x -x G | + |y -y G | = h. We can write α(x, y) = α G (x, y) + O(h 2 ), α G (x, y) = a + b(x -x G ) + c(y -y G ).
The approximation of α G with the method described above yields a procedure a approximation of α by generalized plane waves.

UWVF and generalized plane waves

We now described the introduction of generalized plane waves in the UWVF method.

Notation

The function space for the UWV formulation is denoted V as

V = ∏ k∈[[1,N h ]] L 2 (∂ Ω k ), equipped with the Hermitian product (x, y) = ∑ k ∂ Ω k 1 γ x k y k . It defines a norm: x = (x, x).
In particular for any operator A ∈ L (V ), the norm is

A = sup x =0 Ax x .
Remark 3.1. The space V depends on the mesh.

Moreover: if Ω ⊂ R the dimension of V is finite; if Ω ⊂ R d with d 2, the dimension of V is infinite.

A standard ultra weak variational formulation

The ultra weak variational formulation is a convenient reformulation of the initial problem. We need to define

H k (α) = v k ∈ H 1 (Ω k ), (-∆ + α)v k = 0, (Ω k ), ((-∂ ν + iγ)v k ) |∂ Ω k ∈ L 2 (∂ Ω k ) (3.1)
and

H = N h ∏ k=1 H k (α). Theorem 3.1. Let u ∈ H 1 (Ω ) be a solution of problem (1.1) such that ∂ ν k u ∈ L 2 (∂ Ω k ) for any k. Let γ > 0 be a given real number. Then x ∈ V defined by x |∂ Ω k = x k with x k = ((-∂ ν + iγ)u |Ω k ) |∂ Ω k satisfies ∑ k ∂ Ω k 1 γ x k (-∂ ν + iγ)e k -∑ j, j =k Σ k j 1 γ x j (∂ ν + iγ)e k -∑ k,Γ k = / 0 Γ k Q γ x k (∂ ν + iγ)e k = -2i ∑ k ∂ Ω k f e + ∑ k Γ k 1 γ g(∂ ν + iγ)e k , (3.2)
for any e = (e

k ) k∈[[1,N h ]] ∈ H. Conversely, if x ∈ V is solution of (3.2) then the function u defined locally by    u |Ω k = u k ∈ H 1 (Ω k ), (-∆ + α)u k = f |Ω k , (-∂ ν k + iγ)u k = x k , (3.3)
is the unique solution of the problem (1.1).

This result is classical in the context of UWVF. We refer to Cessenat & Després (3); Buffa & Monk (1); Hiptmair et al. (13); Huttunen et al. (16,17); (author?) (Imbert-Gérard & Després). Our main task is to adapt this formulation to the generalized plane waves developed previously. In order to give a more compact formulation useful for further developments, some definitions are required. Definition 3.1. For any f ∈ L 2 (Ω ), let E f be the extension mapping defined by :

E f : V → H, z → e = (e k ) k∈[[1,N h ]] ,
where e is defined ∀k ∈ [[1, N h ]] by the unique solution of the following problem :

(-∆ + α)e k = f (Ω k ), (-∂ ν k + iγ)e k = z k (∂ Ω k ).
Also define E which is the homogeneous extension operator with vanishing right hand side, namely E = E 0 .

Notice that E f is well defined thanks to theorem A.1. Definition 3.2. Let F be the mapping defined by

F : V → V, z → (∂ ν + iγ)E(z) |∂ Ω k k∈[[1,N h ]] .
This operator relates the outgoing and ingoing traces on the boundaries ∂ Ω k .

Definition 3.3. Let Π be the mapping defined by

Π :    V → V, z |Σ k j → z |Σ jk , z |Γ k → Qz |Γ k .
Definition 3.4. If F * denotes the adjoint operator of the operator F, let A be the operator F * Π .

With these notations the problem (3.2) is equivalent Cessenat & Després (3) to Find x ∈ V such that ∀y ∈ V (x, y) -(Π x, Fy) = (b, y), (3.4)
where the right hand side b ∈ V is given by the Riesz theorem

(b, y) = -2i Ω f E(y) + Γ 1 γ gF(y), ∀y ∈ V.
More precisely

• If u is solution of the initial problem (1.1) such that (-∂ ν + iγ)u |∂ Ω k k∈[[1,N h ]] ∈ V , then x = (-∂ ν + iγ)u |∂ Ω k k∈[[1,N h ]] is solution in V of (3.4).
• Conversely if x is solution of (3.4) then u = E f (x) is the unique solution of (3.2). The problem

(3.4) is equivalent to For b ∈ V, find x ∈ V (I -A)x = b. (3.5)
We now give some properties of the operators defined previously. They will be useful for the theoretical study of the method.

Lemma 3.1. The operator Π obviously satisfies Π 1 for any complex function Q such that |Q| 1.

Lemma 3.2. The operator F is an isometry.

Proof. For any y ∈ V , let e ∈ H be E(y). Then

Fy 2 = ∑ k∈[[1,N h ]] ∂ Ω k 1 γ |(∂ ν + iγ)e k | 2 = ∑ k∈[[1,N h ]] ∂ Ω k 1 γ |∂ ν e k | 2 + γ|e k | 2 + 2ℑ(∂ ν e k • e k ),
y 2 = ∑ k∈[[1,N h ]] ∂ Ω k 1 γ |(-∂ ν + iγ)e k | 2 = ∑ k∈[[1,N h ]] ∂ Ω k 1 γ |∂ ν e k | 2 + γ|e k | 2 -2ℑ(∂ ν e k • e k ). Since ∂ Ω k ∂ ν e k • e k = Ω k |∇e k | 2 + α|e k | 2 ∈
R, one gets that Fy 2 = y 2 . This implies the result.

Proposition 3.1. The operator A satisfies A 1.

This operator also satisfies the following property.

Proposition 3.2. The operator I -A is injective.

Proof. Let x ∈ V such that (I -A)x = 0, which means x = F * Π x. Define z ∈ V such that z = Π x, then F * z = x so that Π F * z = z. Then define u ∈ H such that for all k ∈ [[1, N h ]] -∆ u + αu = 0, (Ω k ), (∂ ν + iγ)u = z |∂ Ω k , (∂ Ω k ).
(3.6)

In order to identify

F * z, define y ∈ V such that ∀k ∈ [[1, N h ]], y k = (-∂ ν + iγ)u |Ω k . It is known that ∀v ∈ V , there exists w ∈ H such that w = E(v), which means w satisfies -∆ w + αw = 0, (Ω k ), (-∂ ν + iγ)w = v |∂ Ω k , (∂ Ω k ).
(3.7)

Then (y, v) = ∑ k∈[[1,N h ]] ∂ Ω k 1 γ (-∂ ν + iγ)u |Ω k • (-∂ ν + iγ)w |∂ Ω k , = ∑ k∈[[1,N h ]] ∂ Ω k 1 γ ∂ ν u • ∂ ν w + γu • w + i∂ ν u • w -iu • ∂ ν w, (z, Fv) = ∑ k∈[[1,N h ]] ∂ Ω k 1 γ (∂ ν + iγ)u |Ω k • (∂ ν + iγ)w |∂ Ω k , = ∑ k∈[[1,N h ]] ∂ Ω k 1 γ ∂ ν u • ∂ ν w + γu • w -i∂ ν u • w + iu • ∂ ν w.
On the other hand, from (3.6) and (3.

7) for all k ∈ [[1, N h ]] ∂ Ω k ∂ ν u • w = Ω k ∇u • ∇w + Ω k αu • w, ∂ Ω k u • ∂ ν w = Ω k ∇u • ∇w + Ω k αu • w, so that ∂ Ω k -∂ ν u • w + u • ∂ ν w = 0. As a consequence ∀v ∈ V, (y, v) = (z, Fv), which exactly means that y = F * z. Since Π F * z = z, it leads to Π y = z.
To conclude let's read this last equation in terms of the function u defined in (3.6).

∀(k, j) ∈ [[1, N h ]] 2 , (-∂ ν + iγ)u |Σ jk = (∂ ν + iγ)u |Σ k j , Q(-∂ ν + iγ)u |Γ k = (∂ ν + iγ)u |Γ k ,
so that both u and ∂ ν u are continuous along every interface Σ k j , and now

-∆ u + αu = 0, (Ω ), (∂ ν + iγ)u = Q(-∂ ν + iγ)u, (∂ Ω ).
Thanks to the preliminary result, u is the unique solution of the corresponding (1.1) problem : it is the 0 solution. Then z = 0, and so x = 0. The proof is ended.

An abstract discretization procedure

The next step consists in the discretization of equation (3.4). This could be treated thanks to a standard Galerkin method which is presented below. That is we consider a subspace V h ⊂ V with finite dimension. We seek the discrete solution x h ∈ V h such that

∀y h ∈ V h , (x h , y h ) -(Π x h , Fy h ) = (b, y h ).
(3.8)

The definition of the operator F, through the operator E, is linked to the functional space H ; this fact means that solutions of the homogeneous equation are needed then. In other words this abstract Galerkin procedure needs a companion constructive procedure to design the basis functions to generate V h . Before describing in the next section what is our proposition to make such a Galerkin method effective, we explain below why such a Galerkin approach (3.8) yields a well posed discrete problem. We provide here an analysis of this well known fact which is slightly different from what can be found in the literature Després ( 4 Proof. One has by definition

|||x||| 2 = x 2 + Ax 2 -2R(x, Ax). Since A 1 then |||x||| 2 2 x 2 -R(x, Ax) V = 2R ((I -A) x, x) V = 2Ra(x, x).
The coercivity is proved. The skewed bicontinuity is evident from Cauchy-Schwarz inequality applied to a(x, y) = ((I -A) x, y).

Proposition 3.3. Assume there exists x solution of the problem (3.5). Then any discrete solution x h satisfies the inequality |||x -

x h ||| 2 inf z h ∈V h x -z h . (3.9) Proof. By construction a(x -x h , y h ) = 0 ∀y h ∈ V h . So a(x -x h , x -x h ) = a(x -x h , x -z h ) with z h = y y -x h .
It ends the proof with the coercivity and skewed bicontinuity of lemma 3.3.

Lemma 3.4. For all b ∈ V , the discrete solution x h exists and is unique.

Proof. If x h exists, it is solution of a linear system, the dimension of the system being the dimension of the discrete subspace V h . Therefore it is sufficient to check that if a(x h , y h ) = 0 for all y h ∈ V h , then x h = 0. Apply the inequality (3.9) with the choice x = b = 0. It yields x h 2 inf z h ∈V h z h = 0.

The new method

We consider that the generalized plane waves ϕ l k have been constructed with the procedure described previously. The local discrete space is

W k = Span (-∂ ν + iγ)ϕ l k 1 l p(k) ⊂ L 2 (∂ Ω k ).
The global discrete space V q ⊂ V is defined by : V q = ∏ 1 k N h W k . Regarding these definitions, one sees that the basis functions are defined on the boundaries of the mesh, and that they have compact support. That is the shape function defined from

ϕ l k has support in L 2 (∂ Ω k ) and vanishes in L 2 (∂ Ω k ′ ) for k ′ = k. It is therefore convenient to define the trace v l k ∈ V by v l k = (-∂ ν + iγ)ϕ l k on L 2 (∂ Ω k ), and v l k = 0 on L 2 (∂ Ω k ′ ) k ′ = k.
An equivalent way to define W k and V q could be

W k = Span(v l k ) 1 l p(k) and V q = Span(v l k ) 1 k p(k), 1 p N h .
Next define what are the generalizations of operators E and F in this context. Let E q ∈ L (V q , H) be the discrete mapping defined ∀k ∈

[[1, N h ]] and ∀l ∈ [[1, p(k)]] by E q (v l k ) = ϕ l k on H 1 (Ω k ), and v l k = 0 on H 1 (Ω k ′ ) k ′ = k. (3.10) Similarly define F q ∈ L (V q ,V ), ∀k ∈ [[1, N h ]] and ∀l ∈ [[1, p(k)]], by F q (v l k ) = (∂ ν + iγ)(ϕ l k ) on L 2 (∂ Ω k ), and v l k = 0 on L 2 (∂ Ω k ′ ) k ′ = k.
With these notations and definitions, the abstract UWVF with generalized plane waves is defined as follows.

Definition 3.6. (UWVF method with generalized plane waves) Find x h ∈ V q such that

∀y h ∈ V q , (x h , y h ) V -(Π x h , F q y h ) V = (b q , y h ) V (3.11)
with the right hand side given by

(b q , y h ) V = -2i Ω f E q (y h ) + Γ 1 γ gF q (y h ), ∀y h ∈ V q .
(3.12)

Numerical analysis of the method

In this section we desire to provide tools for the proof of the convergence of the discrete solution defined by (3.11) to the exact solution. Since the discrete method (3.11) can be viewed as a convenient modification of the bilinear form (3.8), it is not surprising that that the convergence analysis strongly relies on the second Strang's lemma as it is the case for non conforming finite element methods Farhat et al. (6). However the technicalities attached to ultra weak formulations are such that the convergence proof will be completed only in dimension one. This is due to the fact that some uniform coercivity properties which are part of the second Strang's lemma are easy to prove in dimension one, see proposition 4.2, but are open problems in higher dimension.

Simplified notations in dimension one

Let the order of approximation q be a given number. Assume that one has two polynomials P k,1 and P k,2 for all k ∈ [[1, N h ]]. The corresponding basis functions and coefficients are denoted ϕ k,1 , α k,1 and ϕ k,2 , α k,2 . For the sake of simplicity, the basis functions space will be now denoted by {ϕ j } j∈ [[1,2N h ]] and the corresponding coefficients

D = {α j } j∈[[1,2N h ]] ; {z j } j∈[[1,2N h ]]
will denote the corresponding traces, i.e.

∀ j ∈ [[1, 2N h ]], z j = {(-∂ ν + iγ)ϕ j|∂ Ω k } k∈[[1,N h ]] .
The family {z j } j∈ [[1,2N h ]] is a basis of the functional space V q . A fundamental property is that

V q = V only in dimension one.
This will greatly reduce the technicalities of the proof. In fact, the lemmas 4.1 and 4.3 rely on the fact that in dimension one dimV q = 2

Preliminary results

For the sake of completeness, here are classical results useful for the study of this new method. The proofs are postponed to the appendix. Then there exists two constants h 0 and C which depend of α L ∞ (O) and γ such that ∀h < h 0

w L 2 (O) C √ h g L 2 (∂ O) , (4.2)
Remark that the existence and uniqueness of the solution is given by theorem A.1. We will also need a result on the approximation error between the problem

-∆ w + αw = f , (O), (-∂ ν + iγ) w = g, (∂ O), (4.3)
and the modified problem

-∆ w + α h w = f , (O), (-∂ ν + iγ) w = g, (∂ O), (

The discrete problem

This paragraph is devoted to showing that the operator F q described in section 3.4 is an approximation of the operator F up to the order q + 1 in h. Consider the following problem Find x h ∈ V q such that (I -A q )x h = b, (4.6)

where A q = (F q ) * Π . Here h and q are given. This result relies on a preliminary lemma.

Lemma 4.1. Let q 2. Suppose h is small enough and basis functions are constructed as described in paragraph 2.2.2. For all k ∈ [[1, N h ]], there exists a constant C independent k such that ∀z ∈ V q and ∀k ∈

[[1, N h ]] ∑ j∈{1,2} |x j | z j L 2 (∂ Ω k ) C ∑ j∈{1,2} x j z j L 2 (∂ Ω k ) . Proof. Set k ∈ [[1, N h ]] and z = x 1 z 1 + x 2 z 2 .
First x j can be written as a function of z. This is a priori possible using {w j } j∈{1,2} which is the dual basis of {z j } j∈{1,2} . For all ( j, l) ∈ {1, 2} 2 , the dual function w j is defined by

(w j , z l ) V = δ jl , (4.7)
where δ denotes the Kronecker symbol. The proof proceeds in several steps.

First step. One has that x j = (z, w j ) V , therefore

∑ j∈{1,2} |x j | z j ∑ j∈{1,2}
z j w j z .

So the claim is proved provided the term between parentheses can be estimated.

Second step: estimation of ∑ j∈{1,2} z j w j . From (4.7) it turns out that

w 1 = -z 2 2 |(z 1 , z 2 )| 2 -z 1 2 z 2 2 z 1 + (z 1 , z 2 ) |(z 1 , z 2 )| 2 -z 1 2 z 2 2 z 2 , w 2 = (z 1 , z 2 ) |(z 1 , z 2 )| 2 -z 1 2 z 2 2 z 1 - z 1 2 |(z 1 , z 2 )| 2 -z 1 2 z 2 2 z 2 , so that ∑ j∈{1,2} z j w j 2 z 1 2 z 2 2 z 1 2 z 2 2 -|(z 1 , z 2 ) 2 | .

Let us set for convenience

A = |(z 1 ,z 2 )| z 1 z 2 so that ∑ j∈{1,2} z j w j 2 1 1-A 2 .
It means that the whole proof relies on an upper bound for A.

Third step: end of the proof. By definition

(z j ) |∂ Ω k = (-∂ ν + iγ)e P j |∂ Ω k
. By construction P j (x k+1/2 ) = 0 for j = 1, 2, P ′ 1 (x k+1/2 ) = 0 and P ′ 2 (x k+1/2 ) = 1. Since by construction all derivatives of P 1 and P 2 are uniformly bounded, one has P j (x) = O(h) for j = 1, 2, P ′ 1 (x) = O(h) and P ′ 2 (x) = 1 + O(h) when h goes to 0 and for all x ∈ [x k , x k+1 ].

So one can estimate

||z 1 || 2 = 1 γ -P ′ 1 (x k+1 ) + iγP 1 (x k+1 ) 2 |exp (P 1 (x k+1 ))| 2 + 1 γ P ′ 1 (x k ) + iγP 1 (x k ) 2 |exp (P 1 (x k ))| 2 = 1 γ -P ′ 1 (x k+ 1 2 ) + iγP 1 (x k+ 1 2 ) 2 + 1 γ P ′ 1 (x k+ 1 2 ) + iγP 1 (x k+ 1 2 ) 2 + O(h), that is ||z 1 || 2 = 2γ + O(h).
With the same method one obtains

||z 2 || 2 = 2 1 + γ 2 γ + O(h) = 1 + γ 2 γ 2 2γ + O(h),
and

(z 1 , z 2 ) = 1 γ (-P ′ 1 (x k+1/2 ) + iγ)(-P ′ 2 (x k+1/2 ) + iγ) + 1 γ (P ′ 1 (x k+1/2 ) + iγ)(P ′ 2 (x k+1/2 ) + iγ) + O(h) that is (z 1 , z 2 ) = 2γ + O(h). Therefore A 2 = γ 2 1+γ 2 + O(h).
It proves the claim for h sufficiently small. Final comment. By construction the polynomials designed in dimension one in section 2.2.2 by the approximation of the Taylor expansion (2.7) are such that all their coefficients are uniformly bounded up to order q for all cells in the domain. This is why the error O(h) in the above analysis is uniform with respect to the cell index k, which is therefore not indicated. This is not true if one constructs the polynomials with the method constructed in the counter example (2.8).

Lemma 4.2. For small h and considering the basis functions constructed as described in paragraph 2.2.2, there exists a constant C F q -F Ch q+1 (4.8)

Proof. For all j ∈ [[1, 2N h ]],
the function ϕ j is by construction ϕ j = E q (z j ) such that

ϕ j ∈ {ϕ l } l∈[[1,2N h ]] satisfies ∀k ∈ [[1, N h ]] z j = (-∂ ν + iγ)ϕ j , (∂ Ω k ), -d 2 dx 2 + α j ϕ j = 0, (Ω k ).
We also define ψ j = E(z j ) such that and the equation with the exact coefficient α

ψ j ∈ H satisfies ∀k ∈ [[1, N h ]] z j = (-∂ ν + iγ)ψ j , (∂ Ω k ), -d 2 dx 2 + α ψ j = 0, (Ω k ). Then |(F q -F)z j | 2 = |(∂ ν + iγ)(ϕ j -ψ j )| 2 , = |(-∂ ν + iγ)(ϕ j -ψ j )| 2 + 2ℜ(iγ(ϕ j -ψ j )∂ ν (ϕ j -ψ j )), = -2γℑ((ϕ j -ψ j )∂ ν (ϕ j -ψ j )),
since ϕ j and ψ j satisfy the same boundary condition:

(-∂ ν + iγ)(ϕ j -ψ j ) = 0.
Then on the only one element where z j is non zero numbered k = k( j)

∂ Ω k 1 γ |(F q -F)z j | 2 = -2ℑ ∂ Ω k (ϕ j -ψ j )∂ ν (ϕ j -ψ j ), = -2ℑ Ω k (ϕ j -ψ j ) d 2 dx 2 (ϕ j -ψ j ) -2ℑ Ω k d dx (ϕ j -ψ j ) 2 , -2ℑ Ω k (ϕ j -ψ j )(α j ϕ j -αψ j ),
since both ϕ j and ψ j satisfy homogeneous equations. Then

∂ Ω k 1 γ |(F q -F)z j | 2 -ℑ Ω k (α j + α)|ϕ j -ψ j | 2 + Ω k (α j -α)(ϕ j -ψ j )(ϕ j + ψ j ) , α j + α L ∞ (Ω k ) ϕ j -ψ j 2 L 2 (Ω k ) + α j -α L ∞ (Ω k ) ϕ j -ψ j L 2 (Ω k ) ϕ j L 2 (Ω k ) + ψ j L 2 (Ω k ) ,
thanks to Cauchy-Schwarz inequality. On the other hand, from (4.2) and (4.5) for small hs

ϕ j -ψ j L 2 (Ω k ) Ch 3 2 z j L 2 (∂ Ω k ) α -α j L ∞ (Ω k ) , ϕ j L 2 (Ω k ) C √ h z j L 2 (∂ Ω k ) , ψ j L 2 (Ω k ) C √ h z j L 2 (∂ Ω k ) ,
and α j + α L ∞ (Ω k ) is bounded as noticed in remark 2.1. So for small h

(F q -F)z j 2 L 2 (∂ Ω k ) C ′ h 2 α j -α 2 L ∞ (Ω k ) z j 2 L 2 (∂ Ω k ) ,
where still k denotes k( j). Now for all k

∈ [[1, N h ]] let L(k) be the set of indexes j ∈ [[1, 2N h ]] such that Ω k is the support of z j .
Hence, for all z ∈ V q then z |∂ Ω k = ∑ l∈L(k) x l z l where both z l s vanish on ∂ Ω j for all j = k, it yields

(F q -F)z L 2 (∂ Ω k ) ∑ l∈L(k) |x l | (F q -F)z l L 2 (∂ Ω k ) Ch max l∈L(k) α l k -α L ∞ (Ω k ) ∑ l∈{1,2} |x l | z l L 2 (∂ Ω k ) .
Thanks to lemma 4.1 it means that

(F q -F)z L 2 (∂ Ω k ) √ C ′ h max l∈L(k) α l k -α L ∞ (Ω k ) z L 2 (∂ Ω k ) .
Going back to the definition of the V norm for all z ∈ V

(F q -F)z Ch max j∈[[1,2N h ]] α j -α L ∞ (Ω k ) z , which exactly means F q -F Ch max j∈[[1,2N h ]] α -α j L ∞ (Ω k ) .
The result then comes from equation (2.6) ensured by the construction of approximated coefficients α j s.

Some norms

The whole point of this paragraph is to define a useful norm to adapt the second Strang lemma.

Lemma 4.3. There exists a constant C such that for all x ∈ V : Ch 3/2 x (I -A)x .

Remark that, in dimension one, the dimension of the space V is finite, so all the norms are equivalent; but the constants in the continuity inequalities depend on h, and this lemma specifies the dependence in this mesh parameter. Proof.

First step Take x ∈ V , and define b = (I -A)x. In order to interpret this equality in V , define u = E(x)

and w = E(b), so that (u, w) ∈ H × H and ∀k ∈ [[1, N h ]] -d 2 dx 2 + α u = 0, (Ω k ), (-∂ ν + iγ)u = x k , (∂ Ω k ), ∀k ∈ [[1, N h ]] -d 2 dx 2 + α w = 0, (Ω k ), (-∂ ν + iγ)w = b k , (∂ Ω k ). Since F is an isometry one has Fx -Π x = Fb.
It means on every interface

∀k ∈ [[1, N h ]], (-∂ ν + iγ)u |Ω k (x k ) -1 k =1 (-∂ ν + iγ)u |Ω k-1 (x k ) = (-∂ ν + iγ)w |Ω k (x k ), (∂ ν + iγ)u |Ω k (x k+1 ) -1 k =N h (∂ ν + iγ)u |Ω k+1 (x k+1 ) = (∂ ν + iγ)w |Ω k (x k+1 ).
This leads to a system of jump conditions on the interfaces

               (-∂ ν + iγ)u |Ω 1 (x 1 ) = (-∂ ν + iγ)w |Ω 1 (x 1 ), ∀k ∈ [[2, N h ]], d dx u |Ω k-1 - d dx u |Ω k (x k ) = 1 2 (-∂ ν + iγ)w |Ω k -(∂ ν + iγ)w |Ω k-1 (x k ), u |Ω k -u |Ω k-1 (x k ) = 1 2iγ (-∂ ν + iγ)w |Ω k -(∂ ν + iγ)w |Ω k-1 (x k ), (∂ ν + iγ)u |Ω N h (x N h +1 ) = (∂ ν + iγ)w |Ω N h (x N h +1 ).
(4.9) Considering U 0 and U 1 the two fundamental solutions of the homogeneous equation such that

-d 2 dx 2 + α u = 0 on Ω , then u satisfies ∀k ∈ [[1, N h ]], u |Ω k = δ k 0 U 0 + δ k 1 U 1 , (4.10)
where

(δ k 0 , δ k 1 ) k∈[[1,N h ]]
completely determine u ∈ H. Plugging (4.10) in (4.9), and defining

             λ 0 = (-∂ ν + iγ)w |Ω 1 (x 1 ), ∀k ∈ [[2, N h ]], λ k-1 = 1 2 (-∂ ν + iγ)w |Ω k -(∂ ν + iγ)w |Ω k-1 (x k ), µ k-1 = 1 2iγ (-∂ ν + iγ)w |Ω k -(∂ ν + iγ)w |Ω k-1 (x k ), µ N h = (∂ ν + iγ)w |Ω N h (x N h +1 ), then            (-∂ ν + iγ)U 0 (x 1 )δ 1 0 + (-∂ ν + iγ)U 1 (x 1 )δ 1 1 = λ 0 , ∀k ∈ [[2, N h ]], d dx U 0 (x k )(δ k-1 0 -δ k 0 ) + d dx U 1 (x k )(δ k-1 1 -δ k 1 ) = λ k-1 , U 0 (x k )(δ k-1 0 -δ k 0 ) + U 1 (x k )(δ k-1 1 -δ k 1 ) = µ k-1 , (∂ ν + iγ)U 0 (x N h +1 )δ N h 0 + (∂ ν + iγ)U 1 (x N h +1 )δ N h 1 = µ N h . (4.11) Given the change of variable ∀k ∈ [[1, N h -1]] D k 0 = δ k 0 -δ k+1 0 , D k 1 = δ k 1 -δ k+1 1 , (4.12)
the system (4.11) gives a linear system with unknowns

(D k 0 , D k 1 ) k∈[[1,N h -1]] . Defining the Wron- skian W 0 = U 1 d dx U 0 -U 0 d dx U 1 -which is non zero -the solution is ∀k ∈ [[1, N h -1]]        D k 0 = 1 W 0 λ k U 1 (x k+1 ) -µ k d dx U 1 (x k+1 ) , D k 1 = 1 W 0 µ k d dx U 0 (x k+1 ) -λ k U 0 (x k+1 ) .
Then the structure of the system (4.11) is

       αδ 1 0 + β δ 1 1 = λ 0 , δ k 0 -δ k+1 0 = D k 0 , ∀k ∈ [[1, N h -1]], δ k 1 -δ k+1 1 = D k 1 , ∀k ∈ [[1, N h -1]], γδ N h 0 + ηδ N h 1 = µ N h . Eliminating (δ k 0 , δ k 1 ) k∈[[1,N h -1]] it yields δ 1 0 = ∑ N h -1 k=1 D k 0 + δ N h 0 , δ 1 1 = ∑ N h -1 k=1 D k 1 + δ N h 1 ,
and αδ

N h 0 + β δ N h 1 = L, γδ N h 0 + ηδ N h 1 = µ N h , (4.13) with L = λ 0 -(-∂ ν + iγ)U 0 (a) N h -1 ∑ k=1 D k 0 -(-∂ ν + iγ)U 1 (a) N h -1 ∑ k=1 D k 1 . (4.14)
The determinant of the system (4.13) is

W 1 = (-∂ ν +iγ)U 0 (a)(∂ ν +iγ)U 1 (b)-(∂ ν +iγ)U 0 (b)(-∂ ν + iγ)U 1 (a)
. If it were zero, then its columns would be linearly dependent, say

a 0 C 1 + a 1 C 2 = 0 ; this would mean (∂ ν + iγ)(a 0 U 0 + a 1 U 1 )(x 1 ) = 0 and (∂ ν + iγ)(a 0 U 0 + a 1 U 1 )(x N h ) = 0 so that u = a 0 U 0 + a 1 U 0 would satisfy -u ′′ + αu = 0, (∂ ν + iγ)u = 0.
Then u would be the unique solution (zero) of this last system, which is not possible since U 0 and U 1 are independent. Then W 1 is non zero. One finally obtains that

                         δ N h 0 = 1 W 1 L(∂ ν + iγ)U 1 (b) -µ N h (-∂ ν + iγ)U 1 (a) , δ N h 1 = 1 W 1 µ N h (-∂ ν + iγ)U 0 (a) -L(-∂ ν + iγ)U 1 (a) , ∀k ∈ [[1, N h -1]] δ k 0 = δ N h 0 + N h -1 ∑ j=k D j 0 , δ k 1 = δ N h 1 + N h -1 ∑ j=k D j 1 .
(4.15)

Now u is completely known.

Second step The next step is the estimation of the coefficients

(δ k 0 , δ k 1 ) k∈[[1,N h ]] using (4.15
). Since F is an isometry, and λ k and µ k are linear combinations of the components of Fb

   ∀k ∈ [[0, N h -1]], |λ k | √ γ b , ∀k ∈ [[1, N h ]], |µ k | 1 √ γ b . (4.16)
Thus from (4.12) and (4.16), with C depending on U 0 , U 1 , γ and W 0 ,

N h -1 ∑ k=1 D k 0 CN h b , , N h -1 ∑ k=1 D k 1 CN h b .
From (4.14), |L| CN h b , and since |µ N h | C b one has from (4.15)

δ N h i CN h b , ∀i ∈ {0, 1}, and next for k ∈ [[1, N h -1]]: δ k i δ N h i + N h -1 ∑ k=1 D k i CN h b . Then all δ terms satisfy δ k i CN h b for i ∈ {0, 1} and k ∈ [[1, N h ]].
End of the proof A last calculus leads to the following inequalities

x 2 = ∑ k∈[[1,N h ]] δ k 0 (-∂ ν + iγ)U 0 + δ k 1 (-∂ ν + iγ)U 1 2 L 2 (∂ Ω k ) ∑ k∈[[1,N h ]] 2C(|δ k 0 | + |δ k 1 |) 2 C ∑ k∈[[1,N h ]] N 2 h b 2 C b 2 N 3 h , so that x Ch -3/2 b .
Definition 4.1. Let us define |||x||| q = (I -A q )x for all x ∈ V . This is a norm under the condition of the next proposition.

Proposition 4.1. Let q 2 be given and let h be small enough. There exists a constant C > 0 such that Ch 3/2 x |||x||| q , ∀x ∈ V.

(4.17)

Proof. One has ∀x ∈ V, (I -A)x (I -A q )x + (A q -A)x (I -A q )x +Ch q+1 x . So (I -A)x V -Ch q+1 x (I -A q )
x , ∀x ∈ V . Then lemma 4.3 concludes the proof since h q+1 < h 3/2 for h small enough.

Proposition 4.2. There exists a constant h 1 > 0 such that the bilinear form a q (x, y) = ((

I -A q )x, y) is uniformly coercive, i.e. ∀h h 1 |||x||| 2 q 3R (a q (x, x)) , ∀x ∈ V.
Proof.

One has |||x||| 2 q x 2 -2R (A q x, x) + A q x 2 . Since

A q x Ax + (A q -A)x (1 +Ch q+1 ) x (4.18)
there exists another constant denoted as

C ′ > 0 such that A q x 2 (1 +C ′ h q+1 ) x 2 . Therefore |||x||| 2 q 2 x 2 +C ′ h q+1 x 2 -2R (A q x, x) ,
that is |||x||| 2 q -C ′ h q+1 x 2 2R (a q (x, x)). For small h since q > 3/2 and due to the proposition 4.1 one has

C ′ h q x 2 C ′ h 1/3 (h 3/2 -h q ) x 2 h 1/3 |||x||| 2 q , then 2 3 |||x||| 2 q |||x||| 2 q -C ′ h q x 2 .
Combined with the previous inequality it proves the claim.

Convergence

The main convergence result is an adapted version of Strang second lemma with the ||| • ||| q norm. Theorem 4.3. Suppose that q 2 and h min(h 0 , h 1 ). Denote x ∈ V the solution of the exact problem (3.5) in dimension one and x h ∈ V the solution of the discrete problem (4.6). Then there exists a constant C > 0 such that

|||x -x h ||| q Ch -3/2 inf y h ∈V |||x -y h ||| q + sup w h ∈V -{0} |a q (x, w h ) -f q (w h )| w h , (4.19)
where f q (y) = (b q , y) V .

The proof relies on the following intermediate result already proved in (4.18).

Lemma 4.4. The operator A q satisfies A q 1 +Ch q+1 .

Proof. Of theorem 4.3

• The first remark is the uniform coercivity with respect to |||.||| q needed in the second Strang lemma. It is proved in proposition 4.2.

• The second step consists in characterizing the uniform continuity of a q . For all (x, y)

∈ V 2 |a q (x, y)| = |((I -A q )x, y)|, |||x||| q y .
Using (4.17) one has w h Ch -3/2 |||w h ||| q for some constant C, so that for small h ∀(x, y) ∈ V 2 , |a q (x, y)| Ch -3/2 |||x||| q |||y||| q .

• The last step is the inequality itself. The triangular inequality yields

|||x -x h ||| q |||x -y h ||| q + |||x h -y h ||| q , ∀y h ∈ V.
On the other hand proposition 3.3 shows that

1 3 |||x h -y h ||| 2 q a q (x h -y h , x h -y h ) , a q (x -y h , x h -y h ) + a q (x -x h , x h -y h ) , Ch -3/2 |||x -y h ||| q |||x h -y h ||| q + a q (x, x h -y h ) -b q (x h -y h ) .
As

w h = x h -y h ∈ V , then 1 3 |||x h -y h ||| q Ch -3/2 |||x -y h ||| q + |a q (x, w h ) -b q (w h )| w h w h |||w h ||| q .
Using one more time w h Ch -3/2 |||w h ||| q , it yields the desired result.

We now have to estimate the error defined by

D h (x, w h ) = |a q (x, w h ) -b q (w h )|, ∀w h ∈ V.
In order to simplify the proof and to match to the physical meaning of the problem, we will assume that the right hand side (3.12) is characterized by f = 0 and g ∈ L 2 (Γ ) ⊂ V .

Lemma 4.5. There exists a constant C > 0 such that Ch q+1 x w h +Ch q+1 g w h since the second term vanishes (I -A)x = b. The third term is bounded using (3.12) like (bb q , w h ) C F -F q g w h Ch q+1 g w h . This gives exactly (4.20). It is now easy to prove the theoretical convergence of the method in dimension one.

∀w h ∈ V -{0}, D h (x, w h ) w h Ch q+1 ( x + g ) . (4.20) Proof. ∀w h ∈ V -{0}, D h (x, w h ) = |((I -A q )x, w h ) V -(b, w h ) V | , |((A -A q )x, w h ) V | + |((I -A)x, w h ) V -(b, w h ) V | + (b -b q , w h ) ,
Theorem 4.4. One has the estimation |||xx h ||| q = O(h q-1/2 ).

(4.21)

Proof. In dimension one the discrete space of approximation is equal to V whatever the method of construction of basis functions is. This is why one can choose y h = x in (4.19). So inf y h ∈V |||xy h ||| q = 0. The remaining term is bounded with (4.20).

It is useful to rewrite this inequality using a norm with the usual scaling

z = ∑ k∈[[1,N h ]] h|z k | 2 .
By construction z = h 1 2 z . Using (4.17) one gets z Ch -1 |||z||| q . Therefore a corollary of the theorem is the estimate of convergence xx h = O(h q-3/2 ).

(4.22)

Numerical examples

All the following examples are linked with Airy functions since it is the physical problem (1.3)-(1.4) we are interested in. We only consider here coefficients β (x) = x and β (x, y) = x, so that in dimension one as in dimension two that Airy functions are exact solutions of the equation. All the linear systems are assembled and solved with Matlab. The parameter γ is set to be constant equal to 1 everywhere.

One dimensional test case

The test problem considered here is the following : on an interval

Ω =]a, b[⊂ R -u ′′ (x) + x u(x) = 0, (]a, b[), (∂ ν + iγ)u(x) = (∂ ν + iγ)Ai(x), ({a, b}),
The points of the uniform mesh are denoted {x k } k∈ [[1,N h +1]] , where N h stands for the number of elements defining the mesh. For a given value of q the basis functions are designed as in paragraph 2.2.2. The solution computed corresponds to an element x h ∈ V . A simple formula to expresses the traces of u h in funcrion of x h 2iγu h = (I + Π )x h + g ({a, b}),

2iγu h = (I + Π )x h ({x k } k∈[[2,N h ]] ).
In all simulations, the accuracy is reported using a discrete l 2 norm so that the relative error is computed as

∑ k∈[[1,N h +1]] |u ex (x k ) -u h (x k )| 2 ∑ k∈[[1,N h +1]] |u ex (x k )| 2 .
Considering the domain Ω =] -5, 5[, one gets the the typical result of figure 2 where we plot the exact analytical Airy function and the numerical solution computed with our method. The rates of convergence are described in figures 1 and 3. The numerical rates of convergence are better than the theoretical estimates. 3.4e-03 -1.99 8.3e-06 -4.00 5.3e-06 -4.00 1.3e-08 -6.00 7.9e-09 -6.00 1024 8.6e-04 -2.00 5.2e-07 -4.00 3.3e-07 -4.00 2.0e-10 -6.00 1.2e-10 -6.00 2048 2.2e-04 -2.00 3.3e-08 -4.00 2.1e-08 -4.00 3.1e-12 -5.99 1.9e-12 -6.00 4096 5.4e-05 -2.00 2.0e-09 -4.00 1.3e-09 -4.00 7.3e-14 -5.43 7.5e-14 -4.69 8192 1.3e-05 -2.00 1.3e-10 -4.00 8.1e-11 -4.00 1.6e-14 -2.21 5.8e-14 -0.37 16384 3.4e-06 -2.00 7.9e-12 -4.01 5.0e-12 -4.01 5.0e-14 1.67 5.0e-14 -0.20

One can also notice that we observe better convergence rates for odd values of q compared to even values. We have no explanation for the moment.

One can see that on the finest meshes the solution is accurate to machine precision for the highest values of parameter q. FIG. 3: Convergence of the method increasing the parameter q, relative discrete L 2 error as a function of the number of elements defining the mesh.

About q convergence

In figure 3, when the number of nodes is fixed, the error decreases when the parameter q 2 increases.To obtain better understanding of this phenomenon, we plot in figure 4 for different values of q and around two points x 0 the Airy function and its approximations thanks to the two basis functions ϕ constructed in section 2.2.2.

We observe that the approximation is uniform in ]x 0 -ε, x 0 + ε[, with ε independent of q.

Two dimensional test case

A first test case in dimension two is presented here. Consider an open set Ω ⊂ R 2 and the following simple problem

-∆ u(x, y) + x u(x, y) = 0, (Ω ), (∂ ν + iγ)u(x, y) = (∂ ν + iγ)Ai(x), (∂ Ω ),
so that the exact solution is again the Airy function Ai. The domain considered here is square and meshed with regular triangles. A comparison between the exact solution and the numerical solution is diplayed in figure 5.

As explained in section 2.3.1, the design of basis functions is easy in the case of a coefficient depending on only one coordinate, performing a one dimension reduction. The basis function ϕ is defined by ϕ(x, y) = e P(x,y) with P(x, y) = p(x) + λ y where λ still has to be defined. Here we chose 5.3e-03 -2.04 1.1e-04 -2.96 2.6e-05 -3.96 7.1e-07 -6.12 1.3e-06 -3.74 12288 1.3e-03 -2.04 1.4e-05 -3.01 1.6e-06 -4.00 1.1e-07 -2.67 1.7e-07 -3.03

λ ∈ i sin 2πk 3 , k ∈ [[1, 3]] .
For each λ the corresponding functions p + and p -are constructed as in the one dimensional case since ϕ being a solution of the homogeneous equation means thate p(x) ′′ + (x -λ 2 )e p(x) = 0.

The other difference with the one dimensional case is the numerical estimation of boundary integrals. It requires numerical quadrature. The quadrature is performed with a given number of points with either Simpson or Boole method. The corresponding results are given in figures 6, 2 and 3. One can observe a clear improvement in the results obtained using Boole formulas compared to the results obtained using Simpson formulas.

Other basis functions

Figures 7 and 4 present the numerical convergence results obtained with basis functions designed with the normalization β 1,± = ± α(x k+1/2 ). Comparing to figures 1 and 3, one can see that the convergence rate is not modified by this new choice, however for a given number of mesh elements the error is smaller when the method is constructed with these new basis functions than with the basis functions described in section 2.2.2. In fact, for a given order q, the numerical results show that the constant underlying FIG. 7: Relative discrete L 2 error as a function of the number of elements defining the mesh, using the normalization β 1,± = ± α(x k+1/2 ). Different curves correspond to increasing order parameter q.

wave length of the wave. The design of new families of generalized plane waves with increased or even better optimal accuracy is of course of major interest for practical applications.

The numerical analysis of the method is also rich of new theoretical questions in dimension higher than two. We distinguish two particular problems. A first problem is to determine the best choice for ∂ P ∂ x (G j ) and ∂ P ∂ y (G j ) where G j is the center of mass of the cell and ϕ = e P is a basis function; this problem corresponds to the choice of β 1 described in section 2.2.2. A second problem, more fundamental from a theoretical perspective, is the generalization in two dimension of the inverse inequalities like those of lemma 4.3 and proposition 4.1. The main difficulty stems from the fact that V is finite dimensional in dimension one and has infinite dimension in higher dimension : as a consequence V q = V in dimension two and more. Up to this difference we think nevertheless that the functional setting that we have developed, which is based on the second Strang's lemma, is still convenient in dimension higher than two. In this context it could be worthwhile to make the connection with the Discontinuous Galerkin formalism developed for example in Huttunen et al. (16).

A. Appendix

For the sake of completeness of this work, we review some very classical results needed for the proof of the convergence of our algorithm.

A.1 On the initial problem

It concerns the solution of the system in bounded domains

-∆ u + αu = f , x ∈ Ω , (∂ ν + iγ) u = Q (-∂ ν + iγ) u + g, x ∈ Γ .
It is necessary to assume that the regularity of the boundary of Ω is sufficient so that a unique continuation principle holds. We do not want to discuss it because it is not in the scope of this work. We refer the reader to Monk (23) p.92. Moreover, for classical results on polygonal domains we refer to Grisvard (11). To simplify here Q is constant. A proof can be found in (author?) (Imbert-Gérard & Després).

Theorem A.1. Let Ω be a bounded domain in R 2 with a Lipschitz and piecewise C 2 boundary Γ . Let f ∈ L 2 (Ω ), g ∈ L 2 (Γ ) and ζ ∈ C such that ℜ(ζ ) = 0. Then there exists a unique solution u ∈ H 1 (Ω ) to the variational formulation

Ω ∇u • ∇v + Ω αuv + iζ Γ uv = Ω f v + Γ gv, ∀v ∈ H 1 (Ω ).
Using the notations of (A.1) with |Q| < 1, then ℜ( 1-Q 1+Q ) = 0 so there exists a unique solution u ∈ H 1 to (1.1), i.e. such that

Ω ∇u • ∇v + Ω αuv + i 1 -Q 1 + Q γ Γ uv = Ω f v + 1 1 + Q Γ gv, ∀v ∈ H 1 (Ω ).
Remark A.1. This result can be generalized to the case where |Q| 1 almost everywhere on Γ and |Q| < 1 on a smooth part of Γ which length is non zero.

A.2 Proof of inequality (4.2)

We need a very classical Poincaré inequality in one dimension. We obtain

   u 2 L 2 (∂ O) 1 γ f L 2 (O) u L 2 (O) + 1 γ g L 2 (∂ O) u L 2 (∂ O) , u ′ 2 L 2 (O) g L 2 (∂ O) u L 2 (∂ O) + β L ∞ (O) u 2 L 2 (O) + f L 2 (O) u L 2 (O) .
The first inequality yields

u 2 L 2 (∂ O) 2 γ f L 2 (O) u L 2 (O) + 1 γ 2 g 2 L 2 (∂ O) .
A standard inequality yields

g L 2 (∂ O) u L 2 (∂ O) 1 2γ g 2 L 2 (∂ O) + γ 2 u 2 L 2 (∂ O) 1 2γ g 2 L 2 (∂ O) + f L 2 (O) u L 2 (O) + 1 2γ g 2 L 2 (∂ O) .
Inserting in the second inequality we obtain

u ′ 2 L 2 (O) 1 γ g 2 L 2 (∂ O) + 2 f L 2 (O) u L 2 (O) + β L ∞ (O) u 2 L 2 (O) .
Then from (A.1)

u 2 L 2 (O) C h 2 γ f L 2 (O) u L 2 (O) + 1 γ 2 g 2 L 2 (∂ O) +h 2 1 2γ g 2 L 2 (∂ O) + 2 f L 2 (O) u L 2 (O) + β L ∞ (O) u 2 L 2 (O)
.

For h small enough we obtain

u 2 L 2 (O) C h γ 2 g 2 L 2 (∂ O) + h 2 γ 2 f 2 L 2 (O) . (A.3)
One can notice that the scaling of this estimate is optimal. Indeed considering that γ is the dimension of the inverse of a length which is evident from the boundary condition, all quantities have the same dimension at inspection of (A.2). Inequality (4.2) is obtained by taking f = 0 in the previous inequality.

Inequality (A.3) yields

e h L 2 (O) C h γ (β h -β )u L 2 (O) C h γ β h -β L ∞ (O) u L 2 (O) .
Using one more time (A.3) to estimate u and regarding γ which is a positive number, we get

e h L 2 (O) C h 3 2 g L 2 (∂ O) + h 2 f L 2 (O) β h -β L ∞ (O) .

  ); Cessenat & Després (3); Gittelson et al. (10); Buffa & Monk (1); Hiptmair et al. (13, 14). Definition 3.5. Let us define the norm |||v||| = (I -A)v for all v ∈ V , and the bilinear form of the formulation (3.4): a(x, y) = (x, y) -(Π x, Fy). Since I -A is injective, ||| • ||| is indeed a norm. In the rest of this paper, R(z) (resp. I(z))stands for the real (imaginary) part of z ∈ C. A fundamental property is Lemma 3.3. The bilinear form is coercive with respect to the norm ||| • ||| |||x||| 2 2R (a(x, x)) ∀x ∈ V, and is bicontinuous in the sense |a(x, y)| |||x||| × y ∀x, y ∈ V.

Theorem 4. 1 .

 1 Let O be a one-dimensional open interval with length h. Let w be the unique solution of -∆ w + αw = 0, (O), (-∂ ν + iγ) w = g, (∂ O).(4.1)

FIG. 2 :

 2 FIG. 2: Plot of the analytical Airy function, and comparison with the numerical solution. Here we used a large (≈ 200) number of cells and two high order generaalized plane waves per cell. One clearly distinguishes between the propagative medium x < 0 and the non propagative medium x > 0.

FIG. 4 :

 4 FIG.4: Approximation of Airy function by corresponding basis functions for different values of q, in the vicinity of x 0 = 0 and x 0 = -4.

Proposition A. 1 .

 1 There exists a constant C such that for all h> 0, all open interval O ⊂ R, for all u ∈ L(O) u L 2 (O) C √ h u L 2 (∂ O) + h u ′ L 2 (O) (A.1) Proof. There exists a ∈ R such that O =]a, a+h[. From u(x) = u(a)+ x a u ′ (t)dt it yields a+h a |u(x)| 2 dx 2h|u(a)| 2 + 2 a+h a ( x a |u ′ (t)|dt) 2 dx, so that u L 2 (O) √ 2h u L 2 (∂ O) + √ 2h u ′ L 2 (O). It gives the result for C = √ 2. Proof. We will show a more general inequality than (4.2) . We use u as test function in the variational formulation (A.1) corresponding to the following problem-u ′′ + β u = f , (O) (-∂ ν + iγ)u = g, (∂ O).

A. 3

 3 Proof of theorem 4.2 Proof. Suppose that u and u h are the solutions of the two following problems-u ′′ + β u = f , (O) (-∂ ν + iγ)u = g, (∂ O).and -u ′′h + β h u h = f , (O) (-∂ ν + iγ)u h = g, (∂ O). Then e h := uu h satisfies -e ′′ h + β h e h = (β h -β )u, (O) (-∂ ν + iγ)e h = 0, (∂ O).

TABLE 1 .

 1 Errors and orders of convergence for different orders of approximation q depending on the number of unknowns N.

		q=2	q=3	q=4	q=5	q=6
	l	Error Rate	Error Rate	Error Rate	Error Rate	Error Rate
	4	9.5e-01 -	9.9e-01 -	8.6e-01 -	8.6e-01 -	NaN	-
	8	9.2e-01 -0.05 9.7e-01 -0.03 9.7e-01 0.18	9.9e-01 0.20	9.9e-01 NaN
	16	7.8e-01 -0.23 9.5e-01 -0.03 9.2e-01 -0.09 9.6e-01 -0.04 9.4e-01 -0.04
	32	6.0e-01 -0.39 3.3e-01 -1.51 2.5e-01 -1.89 1.5e-01 -2.65 1.1e-01 -3.14
	64	2.0e-01 -1.59 3.2e-02 -3.4	2.0e-02 -3.61 3.2e-03 -5.6	2.0e-03 -5.75
	128	5.4e-02 -1.89 2.1e-03 -3.91 1.3e-03 -3.93 5.2e-05 -5.94 3.2e-05 -5.96
	256	1.4e-02 -1.97 1.3e-04 -3.98 8.4e-05 -3.98 8.2e-07 -5.99 5.0e-07 -5.99
	512

TABLE 2 .

 2 Errors

		q=2	q=3	q=4	q=5	q=6
	l	Error Rate	Error Rate	Error Rate	Error Rate	Error Rate
	48	3.1e-01 -	2.1e-02 -	8.4e-02 -	2.2e-02 -	3.5e-02 -
	192	8.6e-02 -1.84 5.6e-03 -1.92 5.9e-03 -3.83 1.6e-03 -3.85 8.9e-04 -5.3
	768	2.2e-02 -1.98 8.9e-04 -2.67 4.1e-04 -3.85 5.0e-05 -4.97 1.8e-05 -5.62
	3072

and orders of convergence depending on the number of unknowns N for the two dimensional case with Simpson quadrature formulas.

g L 2 (∂ O) + h 2 f L 2 (O) αα h L ∞ (O) . (4.5)

1.0e-03 -1.99 7.4e-07 -4.00 1.1e-06 -3.99 6.0e-10 -5.99 1.7e-09 -6.00 512 2.6e-04 -2.00 4.6e-08 -4.00 7.0e-08 -4.00 9.4e-12 -6.00 2.7e-11 -6.00 1024 6.5e-05 -2.00 2.9e-09 -4.00 4.4e-09 -4.00 1.6e-13 -5.92 4.3e-13 -5.98 2048 1.6e-05 -2.00 1.8e-10 -4.00 2.7e-10 -4.00 9.8e-15 -3.99 1.4e-14 -4.95 4096 4.1e-06 -2.00 1.1e-11 -4.00 1.7e-11 -4.00 2.4e-14 1.28 1.9e-14 0.42 8192 1.0e-06 -2.00 7.5e-13 -3.90 1.0e-12 -4.06 1.3e-13 2.43 1.4e-13 2.88 16384 2.6e-07 -2.00 2.1e-13 -1.84 2.0e-13 -2.32 2.1e-13 0.73 2.1e-13 0.61 in estimation (4.22) is much better : for a given number of mesh elements the numerical error can be improved by a factor ≈ 10 2 . Once again the only difference between these two different choices of basis functions relies on the fact that the leading coefficient in P ± does depend or not on the coefficient α. The theoretical tools that developed previously can be adapted without difficulty to this new family of basis functions but the vertical shift visible in figures 3 to 7 will require more research to be fully understood.

As in dimension one, one can see that with Simpson method on the finest meshes the errors increase for q = 5 and 6. It is also supposedly linked with machine precision.

Perspectives

The method proposed in this work has been designed in any dimension : its goal is to increase the accuracy of any plane wave methods in the case the coefficients of the equation are smooth. Preliminary tests in dimensions one and two for h-convergence assess effective gain in accuracy. The evolution of this method may be pursued in at least two directions, numerical and theoretical.

Our first interest is to validate the method on challenging test problems inspired by the physics of reflectometry : in particular we have in mind to use non uniform meshes to reduce the computational burden for problems such that |Ω | 1/d >> λ where |Ω | is the size of the domain and λ the characteristic