Mohamed Amara
email: mohamed.amara@univ-pau.fr

Henri Calandra
email: henri.calandra@total.com

Rabia Djellouli
email: rabia.djellouli@csun.edu

Magdalena Grigoroscuta-Strugaru
email: magdalena.grigoroscuta@inria.fr

A Modied Discontinuous Galerkin Method for Solving Helmholtz Problems

Keywords: Helmholtz equation, discontinuous Galerkin, plane waves, Lagrange multipliers, positive denite matrix, inf-sup condition, waveguide problems équation de Helmholtz, Galerkin discontinu, ondes planes, multiplicateurs de Lagrange, matrice dénie positive, condition inf-sup, guide d'ondes

A new solution methodology is proposed for solving Helmholtz problems in the mid-and high frequency regime. The proposed method falls in the category of the discontinuous Galerkin methods. The primal variable is obtained by solving in parallel a set of well posed local problems. The Lagrange multiplier is the solution of a global positive denite linear system. These two properties are the main features of the proposed method that distinguish it from the existing solution methodologies. Numerical results are presented to illustrate both the stability and the accuracy of the proposed method when applied for solving waveguide-type problems.

Introduction

The Helmholtz equation belongs to the classical equations of mathematical physics that are well understood from a mathematical view point. However, the numerical approximation of the solution is still a challenging problem in spite the tremendous progress made during the last years (see, for example, the recent monograph [START_REF] Magoulès | Computational Methods for Acoustics Problems Saxe-Coburg Publications[END_REF] and the references therein). Indeed, the standard nite element method (FEM) is not well suited for solving Helmholtz problems in the mid-and high-frequency regime since highly oscillating solutions are not accurately approximated by piecewise polynomials. This phenomenon, related to the indeniteness of the Helmholtz operator, is known as the pollution eect [START_REF] Babu²ka | Is the Pollution Eect of the FEM Avoidable for the Helmholtz Equation Considering High Wave Numbers?[END_REF]. In order to maintain a certain level of accuracy, a mesh renement is required and/or higher order FEM are used, leading to a prohibitive computational cost for high wavenumbers.

In response to this challenge, alternative techniques for alleviating the pollution eect were proposed. Numerous of these approaches use the plane waves, since they are expected to better approximate highly oscillating waves. Examples of such methods include the weak element method for Helmholtz equation [START_REF] Rose | Weak element approximations to elliptic dierential equations Numer[END_REF], the Galerkin least-squares method [START_REF] Harari | Galerkin/least-squares nite element methods for the reduced wave equation with non-reecting boundary conditions in unbounded domains[END_REF], the partition of unity method [START_REF] Babu²ka | The partition of unity method Internat[END_REF], the residual free bubbles method [START_REF] Franca | Residual-free bubbles for the Helmholtz equation Internat[END_REF], the least-squares method [START_REF] Monk | A least-squares method for the Helmholtz equation[END_REF], the ultraweak variational method [4] and recently, the discontinuous Galerkin method (DGM) designed by Farhat et al and presented in a series of papers [5,[START_REF] Farhat | A discontinuous Galerkin method with plane waves and Lagrange multipliers for the solution of short wave exterior Helmholtz problems on unstructured meshes[END_REF][START_REF] Farhat | Higher-order extensions of a discontinuous Galerkin method for mid-frequency Helmholtz problems Internat[END_REF]. In the latter method, the solution is approximated at the element mesh level using a superposition of plane waves which results in a discontinuous solution along interior boundaries of the mesh. The continuity is then restored in a weak sense through the use of Lagrange multipliers. The rectangular and quadrilateral elements constructed in [5,[START_REF] Farhat | A discontinuous Galerkin method with plane waves and Lagrange multipliers for the solution of short wave exterior Helmholtz problems on unstructured meshes[END_REF][START_REF] Farhat | Higher-order extensions of a discontinuous Galerkin method for mid-frequency Helmholtz problems Internat[END_REF] clearly outperform the standard Galerkin FEM. For example, for ka ≥ 10 and for a xed level of accuracy, the so-called R-4-1 element reduces the total number of degrees of freedom (dofs) required by the Q1-based nite element discretization for Helmholtz equation by a factor greater or equal to ve. Similar results are obtained for the R-8-2a and R-8-2b elements when compared to the Q2 element, and for Q-16-4 and Q-32-8 when compared to the Q4 element. In spite of this impressive performance, the DGM has three important drawbacks. First, the method has to satisy an inf-sup condition which is translated, in practice, as a compatibility requirement: the number of dofs of the Lagrange multiplier (corresponding to the dual variable) and of the eld (the primal variable) cannot be chosen arbitrarily. The problem here is that there is no theoretical result on how to satisfy this compatibility requirement, except for the simple case of R-4-1 element (see [START_REF] Djellouli | Convergence analysis of a discontinuous Galerkin method with plane waves and Lagrange multipliers for the solution of Helmholtz problems[END_REF]). Hence, for other elements, the existing choices are based on numerical experiments only. The second major issue with the DGM is that it becomes unstable as we rene the mesh. Such instabilities occur because of the singularity of the local problems and, to some extent, to the loss of the linear independence of the plane waves as the step size mesh discretization tends to zero. The latter aects dramatically the stability of the global system due to its ill-conditioning nature. Finally, the DGM exhibits a loss of accuracy for unstructured mesh [START_REF] Farhat | A discontinuous Galerkin method with plane waves and Lagrange multipliers for the solution of short wave exterior Helmholtz problems on unstructured meshes[END_REF].

We propose a new solution methodology for Helmholtz problems, that falls in the category of discontinuous Galerkin methods. The proposed formulation distinguishes itself from existing procedures by the well-posed character of the local problems and by the resulting global system which is associated with a positive denite Hermitian matrix. More specically, the computation domain is subdivided in quadrilateral-or triangular-shaped elements. The solution is approximated, at the element level, by a superposition of plane waves that are solution of the Helmholtz equation. The continuity of the solution at the interior interfaces of the elements is then enforced by Lagrange multipliers. Unlike the DGM, the proposed method does not require the continuity of the normal derivative. Consequently, Lagrange multipliers are introduced to restore in the weak sense the continuity of both the eld and its normal derivative across interior boundaries of the mesh. Such choice leads to solving (a) local boundary value problems that are well posed in the sense of Hadamard [START_REF] Hadamard | Lectures on Cauchy's Problem in Linear Partial Dierential Equations[END_REF] and (b) a global system, whose unknowns are the Lagrange multipliers. The Lagrange multiplier is the solution of a variational problem whose bilinear form can be written into two equivalent expressions. The approximation of both formulations leads to two linear systems corresponding to (a) a non-negative matrix and (b) a positive denite matrix. Note that the proposed technique is a two-step procedure where the local problems are rst solved and then the Lagrange multipliers are evaluated. This two-step approach allows us to consider equally structured and unstructured meshes with either triangular-or quadrilateral-shaped elements. Since the proposed solution methodology resembles in some aspect the DGM, we will refer to it as mDGM (modied Discontinuous Galerkin M ethod).

The remainder of the paper is organized as follows. In Section 2, we introduce general notations and the model problem. Section 3 is devoted to the presentation of mDGM. In Section 4, we present the algebraic framework of the formulation. We compare in Section 5 the numerical performance of both methods: DGM and mDGM. The obtained results clearly indicate that mDGM outperforms DGM in terms of stability and accuracy. Finally, Section 6 concludes this paper.

Preliminaries

In this section, we introduce the model problem and specify the nomenclature and assumptions adopted throughout this paper.

2.1

The mathematical model

We consider the following class of waveguide-type problems:

(BVP) -∆u -k 2 u = f in Ω, ∂ n u = i ku + g on ∂Ω
where Ω ⊂ R 2 is an open bounded region, with smooth boundary ∂Ω, k is a positive number representing the wavenumber, ∂ n is the normal derivative and f and g are regular complex valued functions dened respectively on Ω and ∂Ω. The second equation of BVP is a representation of a class of nonhomogeneous Robin boundary conditions, but other types of boundary condition can be considered.

INRIA

Note that BVP is considered here for its simplicity since it allows us to compute analytically the solution u for a suitable choice of Ω, f and g. Such an expression of u is used when assessing the accuracy of mDGM.

Nomenclature and assumptions

In what follows, we consider a regular triangulation τ h of Ω into quadrilateralor triangular-shaped subdomains K whose boundaries are denoted by ∂K. The step size mesh discretization is denoted by h. We introduce the space of the primal variable:

V = v ∈ L 2 (Ω) ; v| K ∈ H 1 (K) ,
that we equip with the norm:

||v|| V = K∈τ h ||v K || 2 1,K 1 2 , ∀ v ∈ V,
where || • || 1,K is the H 1 -norm on the element K. In addition, we introduce || • || 0,K and | • | 1,K to designate the L 2 -norm and the H 1 -seminorm respectively on the element K.

Note that V contains functions that are discontinuous across interior boundaries since their regularity is only L 2 (Ω). For any v ∈ V, we dene the jump across an interior edge e = ∂K ∩ ∂K by:

[v] = v K -v K .
We introduce the space of the dual variable, corresponding here to Lagrange multipliers, by:

M = µ ∈ K∈τ h L 2 (∂K) ; µ = 0 on ∂K ∩ ∂Ω
and we associate to M the norm given by:

||µ|| M = K∈τ h ||µ K || 2 0,∂K 1 2 , ∀ µ ∈ M,
where µ K designates the restriction of µ to ∂K: µ K = µ| ∂K and || • || 0,∂K is the L 2 -norm on ∂K. Moreover, for any function µ ∈ M, we dene the jump across an interior edge e = ∂K ∩ ∂K by:

[[µ]] = µ K + µ K .

The continuous approach

The basic idea of mDGM is to evaluate u, the solution of BVP, using the following splitting:

u = Φ (λ) + ϕ, (1)
where ϕ and Φ are elements of V and λ ∈ M.

To compute these three quantities, we proceed into two steps:

Step 1 : For all K ∈ τ h and µ ∈ M, we compute ϕ and Φ (µ). This is achieved by solving local Helmholtz problems. This step is called the restriction procedure.

Step 2 : We determine λ ∈ M by solving a global linear system to ensure the continuity in a weak sense of the solution u given by (1) and of the normal derivative of u. This step is called the optimization procedure.

3.1

Step 1: The restriction procedure

As stated earlier, this step is devoted to the computation of ϕ and Φ (µ), for all µ ∈ M, by solving locally Helmholtz problems. More specically, for all K ∈ τ h , we compute ϕ K by solving the following boundary value problem:

(BVP1)        Find ϕ K ∈ H 1 (K) such that: -∆ϕ K -k 2 ϕ K = f in K ∂ n ϕ K = i kϕ K + g on ∂K ∩ ∂Ω ∂ n ϕ K = i αϕ K on ∂K ∩ Ω .
Next, for all µ ∈ V and K ∈ τ h , we compute Φ µ K by solving the boundary value problem given by:

(BVP2)        Find Φ µ K ∈ H 1 (K) such that: -∆Φ µ K -k 2 Φ µ K = 0 in K ∂ n Φ µ K = i kΦ µ K on ∂K ∩ ∂Ω ∂ n Φ µ K = i αΦ µ K + µ K on ∂K ∩ Ω ,
with α ∈ R * + . Note that the presence of α ensures the uniqueness of the solution of BVP 1 and BVP 2, as it will be shown later.

It is easy to verify that the variational formulation of both problems can be expressed in a compact form as follows:

Find Ψ K ∈ H 1 (K) such that: a K Ψ K , v K = L K v K ∀v K ∈ H 1 (K) (2)
where a K (•, •) is a bilinear form given by:

a K v K , w K = K ∇v K • ∇w K dx -k 2 K v K w K dx -i α ∂K∩ Ω v K w K ds -i k ∂K∩∂Ω v K w K ds, ∀ v K , w K ∈ H 1 (K) (3)
and Ψ K is either the solution of BVP1 or BVP2, i.e.:

Ψ K = ϕ K for BVP1 Φ µ K for BVP2 , ∀ µ ∈ M.
The right-hand side L K (•) is given by:

L K v K =      K f v K dx + ∂K∩∂Ω gv K ds for BVP1 ∂K∩ Ω µ K v K ds for BVP2 , ∀ v K ∈ H 1 (K) . ( 4
)
INRIA Consequently, the solving of BVP1 and BVP2 requires solving one linear system with multiple right-hand side. Note that the bilinear form a K is neither Hermitian, nor symmetric. However, it is easy to check that a K (•, •) is continuous on H 1 (K) × H 1 (K) and satises the Gärding inequality in H 1 (K) since

a K v K , v K + k 2 v K 2 0,K = |v| 2 1,K , (5)
where designates the real part. In addition, we have:

Proposition 1. For a xed K ∈ τ h , the variational problem (2) admits a unique solution.

Proof of Proposition 1. Let K be a xed element of τ h . From (5), it follows that the bilinear form a K (•, •) satises the Fredholm alternative. Hence, the uniqueness ensures the existence of Ψ K ∈ H 1 (K), solution of [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graph[END_REF].

To prove the uniqueness, we consider the homogeneous problem associated to the bilinear form a K (•, •) and let w K be its solution. We therefore have:

Find w K ∈ H 1 (K) such that: a K w K , v K = 0 ∀ v K ∈ H 1 (K) (6)
In particular, for v K = w K , we have:

α w K 2 0,∂K∩ Ω + k w K 2 0,∂K∩∂Ω = 0.
Since α > 0, we must have:

w K = 0 on ∂K and ∂ n w K = 0 on ∂K.
Therefore, using the continuation theorem [START_REF] Hörmander | The Analysis of Linear Partial Dierential Operator[END_REF][START_REF] Taylor | Partial Dierential Equations I: Basic Theory[END_REF], we deduce that w K = 0 in K and the problem (2) has a unique solution. 2

Remark 1. The presence of the Robin condition on ∂K ∩ Ω with α > 0 is crucial to ensure that 0 is the only solution of the variational problem [START_REF] Farhat | A discontinuous Galerkin method with plane waves and Lagrange multipliers for the solution of short wave exterior Helmholtz problems on unstructured meshes[END_REF]. Indeed, if α = 0 and K satises ∂K ∩∂Ω = ∅, then the resulting homogeneous Neumann boundary condition on ∂K is not sucient to guarantee the uniqueness of the solution of ( 6), since k 2 may become an interior eigenvalue.

Next, we dene ϕ such that, for all element K in the mesh, the restriction of ϕ to K is ϕ K , the solution of BVP1, i.e. ϕ| K = ϕ K . Similarly, for all element K and for all µ in M, we dene Φ (µ) such that we have Φ (µ) | K = Φ µ K , where Φ µ K is the solution of BVP2. Using the denition of ϕ and Φ (µ) , ∀µ ∈ M we have:

ϕ ∈ V and Φ (µ) ∈ V, ∀µ ∈ M (7)
In summary, Step 1 allows us to compute, for all µ in M:

ϕ + Φ (µ) ∈ V (8)
by solving one variational problem given by ( 2) with dierent right-hand side given by (4).

Step 1 can be viewed, to some extent, as a prediction step.

3.2

Step 2: The optimization procedure

The objective here is to determine λ ∈ M for which the function given by ( 8) is in H 1 (Ω). This requirement can be viewed as a correction stage since we select the best-t Lagrange multiplier λ. The determination of λ is accomplished by solving the following global variational problem:

(VF) Find λ ∈ M such that A (λ, µ) = F (µ), ∀ µ ∈ M, ( 9
)
where the bilinear form A(•, •) is given by:

A(η, µ) = e-interior edge β e e [Φ (η)] [Φ (µ)]ds + e-interior edge e=∂K∩∂K γ e e η K + η K + i α Φ η K + Φ η K µ K + µ K + i α (Φ (µ K ) + Φ (µ K )) ds ( 10
)
and the linear form F (•) is given by:

F (µ) = - e-interior edge β e e [ϕ] [Φ (µ)]ds -i α e-interior edge e=∂K∩∂K γ e e ϕ K + ϕ K µ K + µ K + i α (Φ (µ K ) + Φ (µ K )) ds. ( 11
)
The parameters β e and γ e are two positive numbers that can be viewed as weight parameters. This problem expresses the continuity in the weak sense of the solution and its normal derivative. Note that the bilinear form A is Hermitian. Consequently, only half of the corresponding matrix will be stored.

Remark 2. Unlike the DGM, where only the primal variable u is discontinuous, the mDGM leads to the discontinuity of both variables: the primal variable u and the Lagrange multiplier λ, the dual variable. Consequently, the normal derivative of u is discontinuous. Alternatively, for numerical approximation purpose, we rewrite Eq. ( 10) as follows:

A(η, µ) = e-interior edge β e e [Φ (η)] [Φ (µ)]ds + e-interior edge γ e e [[∂ n Φ (η)]] [[∂ n Φ (µ)]]ds + e⊂∂Ω ω e e (∂ n Φ (η) -i kΦ (η)) (∂ n Φ (µ) -i kΦ (µ))ds, (12)
where the weight parameter ω e is a positive number. Note that the second integral in Eq. ( 12) is equal to the second integral in Eq. ( 10), whereas the third INRIA integral in Eq. ( 12) is in fact equal to 0. Consequently, the right-hand side given by ( 11) is modied depending on the use of Eq. ( 10) or [START_REF] Hörmander | The Analysis of Linear Partial Dierential Operator[END_REF].

The next result states the equivalence between solving BVP and solving the problem arising in the proposed two-step procedure.

Theorem 1.

(i) Let u = Φ(λ) + ϕ, where for all K, ϕ K is solution of BVP1 and Φ λ K is solution of BVP2 with λ solution of VF. Then u is the unique solution of BVP.

(ii) Conversely, let u be the solution of BVP. For each K ∈ τ h , we dene λ by:

λ K = 0 on e ⊂ ∂K ∩ ∂Ω ∂ n u K -i αu K on e ⊂ ∂K ∩ Ω (13)
Let ϕ K be the solution of BVP1 and Φ(λ K ) the solution of BVP2. Then λ is solution of VF and u = Φ(λ) + ϕ.

Remark 3. Eq. ( 13) indicates a clear distinction between mDGM and DGM, in which λ K = ∂ n u K and is continuous along the interior edges.

The algebraic approach

The implementation of mDGM requires rst to introduce two nite-dimensional spaces V h and M h such that V h ⊂ V and M h ⊂ M. Similarly to the DGM formulation, we have considered in this paper spaces of plane waves functions. However, other shape functions satisfying the Helmholtz equation can also be considered. Moreover, unlike the DGM, mDGM allows -in principle -to choose the spaces V h and M h independently.

For any element K ∈ τ h , we denote by

V h (K) (resp. M h (K)) the set of functions of V h (resp. M h ) restricted to K (resp. ∂K). Furthermore, n K (resp. n λ K ) denotes the dimension of V h (K) (resp. M h (K)).
Last, the dimension of M h , which corresponds to the total number of dofs, is denoted by n λ .

We show that when formulated in nite dimensional spaces, the proposed twostep procedure consists in solving linear algebraic systems in each step. Note that in Step 2, the resulting linear system to be solved depends on the approximation of the continuous formulations ( 10) and ( 12) respectively. As stated earlier, equations ( 10) and ( 12) are equivalent only at the continuous level. At the discontinuous level, the second and the third equation of BVP2 are satised in the weak sense.

4.1

Step 1: the restriction procedure

For an element K ∈ τ h and for any

µ K h ∈ M h (K), we denote by ϕ K h ∈ V h (K) and Φ h µ K h ∈ V h (K) the approximation of ϕ K and Φ µ K h respectively.
Similarly to the continuous formulation, ϕ h , Φ h (µ h ) and µ h are given by:

ϕ h | K = ϕ K h , Φ h (µ h ) | K = Φ h µ K h and µ h | K = µ K h , for any element K in the mesh.
To compute ϕ h and Φ h (µ h ), for all K ∈ τ h , we set the variational problem (2) in the nite dimensional space V h (K), that is:

Find Ψ K h ∈ V h (K) such that: a K Ψ K h , v K h = L K v K h , ∀ v K h ∈ V h (K) (14)
where the forms a K (•, •) and L K (•) are given by ( 3) and ( 4) respectively, and

Ψ K h = ϕ K h for BVP1 Φ h µ K h for BVP2, ∀ µ h ∈ M h . ( 15
)
Consequently, the variational problem ( 14)-( 15) can be written in the following matrix form:

K K -k 2 M K -i αS ∂K∩ Ω -i kS ∂K∩∂Ω X K = rhs, (16)
where K K (resp. M K ) is the stiness (resp. mass) matrix at the element level K. S ∂K∩ Ω and S ∂K∩∂Ω are mass-like matrices dened on ∂K ∩ Ω and ∂K ∩ ∂Ω respectively. X K is the vector in C n K whose components are the values of Ψ K h in the basis of V h (K).

The linear system ( 16) possesses the following properties:

All the entries of the corresponding matrix can be evaluated analytically for plane waves shape functions.

The linear system admits a unique solution, even when ∂K ∩ ∂Ω = ∅. Thanks to the positive number α since the presence of the matrix S ∂K∩ Ω guarantees the invertibility of the system. Note that this is not the case for the DGM, for which S ∂K∩ Ω does not appear, leading to possibly a (weakly) singular system when ∂K ∩ ∂Ω = ∅.

The corresponding matrix is neither Hermitian, nor symmetric. This cannot be viewed as a deciency of the approach since the size of the system is small and thus can be solved easily using LU factorization. More specifically, the size of the system is n K × n K , where n K (the number of shape functions at the element level) does not exceed few hundreds.

For an element K ∈ τ h , the number of rhs is n λ K + 1. We must point out that the obtained problems can be solved in parallel since they are independent from an element K to another.

4.2

Step 2: The optimization procedure

In this step, we set the global problem VF in nite dimension. We have:

Find λ h ∈ M h such that: A h (λ h , µ h ) = F h (µ h ), ∀ µ h ∈ M h ( 17
)
where the forms A h (•, •) and F h (•) are obtained from A (•, •) and F (•) respectively by replacing ϕ with ϕ h and Φ (µ h ) with Φ h (µ h ), for µ h ∈ M h . Hence, solving the variational problem [START_REF] Rose | Weak element approximations to elliptic dierential equations Numer[END_REF] comes to solve the following linear algebraic system:

AΛ = b, (18)
INRIA where the entries of the matrix A and of the vector b depend on the expression of the continuous bilinear form A. More specically, when using the bilinear form given by [START_REF] Harari | Galerkin/least-squares nite element methods for the reduced wave equation with non-reecting boundary conditions in unbounded domains[END_REF], the entries of the matrix A and of the vector b are given by:

A lm = e-interior edge β e e [Φ h (µ m )] Φ h (µ l ) ds + e-interior edge e=∂K∩∂K γ e e µ K m + µ K m + i α Φ h µ K m + Φ h µ K m µ K l + µ K l + i α(Φ h µ K l + Φ h µ K l ds ( 19
)
and

b l = - e-interior edge β e e [ϕ h ] Φ h (µ l ) ds -i α e-interior edge e=∂K∩∂K γ e e ϕ K h + ϕ K h µ K l + µ K l + i α(Φ h µ K l + Φ h µ K l ) ds (20) for 1 ≤ l, m ≤ n λ .
On the other hand, when the bilinear form A is given by ( 12), the entries of the corresponding matrix, denoted by A, are dened by:

A lm = e-interior edge β e e [Φ h (µ m )] [Φ h (µ l )]ds + e-interior edge γ e e [[∂ n Φ h (µ m )]] [[∂ n Φ h (µ l )]]ds + e⊂∂Ω ω e e (∂ n Φ h (µ m ) -i kΦ h (µ m )) (∂ n Φ h (µ l ) -i kΦ h (µ l ))ds, (21)
and consequently, the right-hand side, denoted by b, is given by:

b l = - e-interior edge β e e [ϕ h ] [Φ h (µ l )]ds - e-interior edge γ e e [[∂ n ϕ h ]][[∂ n Φ h (µ l )]]ds - e⊂∂Ω ω e e (∂ n ϕ h -i kϕ h -g) (∂ n Φ h (µ l ) -i kΦ h (µ l ))ds (22) for 1 ≤ l, m ≤ n λ .
The unknown Λ is a vector in C n λ whose components are the values of λ h in the basis of M h . Hence, from a numerical point of view, two approaches are possible at Step 2 for determining the Lagrange multiplier:

Approach 1: solving the linear system ( 19)- [START_REF] Taylor | Partial Dierential Equations I: Basic Theory[END_REF].

Approach 2: solving the linear system (21)-( 22).

Note that the matrices A and A are both Hermitian. Next, we show that the system in Approach 1 is positive semi-denite, whereas the one in Approach 2 is positive denite. Therefore, Approach 2 is expected to be more stable.

Proposition 3.

(i) The matrix A is positive semi-denite. Moreover, for Φ h µ K h satisfying, in the strong sense, the following equations:

∀K ∈ τ h , ∀µ h ∈ M h , ∂ n Φ h µ K h = i αΦ h µ K h + µ K h on ∂K ∩ Ω, (23)
and

∀K ∈ τ h , ∀µ h ∈ M h , ∂ n Φ h µ K h = i kΦ h µ K h on ∂K ∩ ∂Ω, (24)
A is a positive-denite matrix.

(ii) The matrix A is positive denite.

Proof of Proposition 3. Let {µ j } 1≤j≤n λ be a basis of M h and let y = t [y 1 , y 2 , ..., y n λ ] ∈ C n λ be an ordinary vector. We set:

η h = 1≤l≤n λ y l µ l . ( 25
)
Therefore, it is easy to verify that:

Φ h (η h ) = 1≤l≤n λ y l Φ h (µ l ). (26)
Consequently, in each K ∈ τ h , Φ h (η h ) satises:

a K Φ h η K h , v K h = ∂K∩ Ω η K h v K h ds ∀v K h ∈ V h (K) . ( 27
)
(i) Using the denition of η h , we have:

y * Ay = e-interior edge β e [Φ h (η h )] 2 L 2 (e) + γ e [[η h + i αΦ h (η h )]] 2 L 2 (e) . ( 28
)
Hence, y * Ay ≥ 0, ∀y ∈ C n λ that is A is a positive semi-denite matrix.

Next, we assume that condition (23) is satised. Let y be a vector in C n λ such that y * Ay = 0. Then, it follows from Eq. (28) that:

e-interior edge

β e [Φ h (η h )] 2 L 2 (e) + γ e [[η h + i αΦ h (η h )]] 2 L 2 (e) = 0. (29)
INRIA Since for any interior edge e, β e > 0 and γ e > 0, we have:

[Φ h (η h )] L 2 (e) = 0 and [[η h + i αΦ h (η h )]] L 2 (e)
= 0 on all interior edges. (30)

Therefore, [Φ h (η h )] = 0 on all interior edges and, using Eq. ( 23), we have also:

[[∂ n Φ h (η h )]] = 0 on all interior edges.
Consequently, Φ h (η h ) ∈ H 1 (Ω) and using (24), we have:

-∆Φ h (η h ) -k 2 Φ h (η h ) = 0 in Ω ∂ n Φ h (η h ) = i kΦ h (η h ) on ∂Ω. ( 31
)
Hence, Φ h (η h ) = 0 because the boundary value problem (31) admits a unique solution. It follows from Eq. ( 23) that η K h = 0, ∀K ∈ τ h . Consequently, η h = 0 and therefore, y l = 0 for all 1 ≤ l ≤ n λ . Thus, y * Ay > 0, ∀y ∈ C n λ \ {0}, that is A is a positive denite matrix.

(ii) Observe that the matrix A satises:

y * Ay = e-interior edge β e [Φ h (η h )] 2 L 2 (e) + γ e [[∂ n Φ h (η h )]] 2 L 2 (e) + e⊂∂Ω ω e ∂ n Φ h (η h ) -i kΦ h (η h ) 2 L 2 (e) ≥ 0,
that is A is a positive semi-denite matrix.

Let y be a vector in C n λ such that y * Ay = 0. Following the same reasoning as for (i), we have:

[Φ h (η h )] = 0 on all interior edges e. Hence, Φ h (η h ) ∈ H 1 (Ω). [[∂ n Φ h (η h )]] = 0 on all interior edges. Using the fact that ∆Φ h η K h ∈ L 2 (K), we deduce that ∆Φ h (η h ) ∈ L 2 (Ω). ∂ n Φ h (η h ) = i Φ h (η h )
on all the boundary edges of the domain.

Using the same argument as for (i), we conclude that y l = 0 for all 1 ≤ l ≤ n λ , and then A is a positive denite matrix. 2

Numerical investigation

In order to illustrate and assess the potential of mDGM for solving eciently Helmholtz problems, we have performed numerical experiments using discrete spaces in which the shape functions are plane waves, as done in DGM [5,[START_REF] Farhat | A discontinuous Galerkin method with plane waves and Lagrange multipliers for the solution of short wave exterior Helmholtz problems on unstructured meshes[END_REF][START_REF] Farhat | Higher-order extensions of a discontinuous Galerkin method for mid-frequency Helmholtz problems Internat[END_REF]. More specically, V h are the spaces introduced in [5]. Once the local space of shape functions, V h (K) is chosen, the Lagrange multiplier is approximated on each edge using a subset or all set of shape functions that occur when evaluating

∂ n v K h -i αv K h , for v K h ∈ V h (K).
From now on, we suppose that Ω is an a × a square domain. We use a uniform partition of Ω in rectangular-shaped elements K. The functions f and g are such that the exact solution u of BVP is a plane wave propagating in a direction d = (cos θ, sin θ). We vary the propagation angle θ in the interval [0, 2π). In order to compare the results obtained with mDGM to those delivered by DGM, we measure, for each propagating angle θ, the relative error using the following modied H 1 norm [5]:

v H 1 =   K v 2 H 1 (K) + e-interior edge [v] 2 L 2 (e)   1 2 , ∀v ∈ V. ( 32
)
Note that (32) is a modied H 1 norm since it takes into account the H 1 norm at the element level and the jump of the numerical solution along the interior interfaces of the mesh. We also use the total relative error, that is the mean value of the relative error obtained when θ ∈ [0, 2π). For all numerical experiments, we have set α = k, and β e = 1 and γ e = h for all the interior edges e. The choice of these parameters results from our numerical investigation, given the lack of theoretical guidelines. We present the results of two classes of numerical experiments: experiments using four plane waves per element, and experiments using eight plane waves per element. All the results are compared to the ones obtained with DGM.

In all numerical experiments, the linear systems ( 19)-( 20) and ( 21)-( 22) are solved using an LU decomposition method for sparse matrices designed by Pardiso [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graph[END_REF][START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with PARDISO[END_REF][START_REF] Schenk | On fast factorization pivoting methods for symmetric indenite systems Elec[END_REF].

Four plane waves per element

We equip each rectangular element with four plane waves positioned as indicated in Figure 1. More specically, for each K ∈ τ h , we consider:

V h (K) = v K h = 1≤p≤4 e i kθp•x u p , θ p = t [cos θ p , sin θ p ] θ p = π/4 + (p -1)π/2, 1 ≤ p ≤ 4, u p ∈ C .
.

As stated earlier, the choice of the basis of M h is related to the computation

∂ n v K h -i kv K h on the edges of the mesh, for v K h ∈ V h (K). Observe that, ∀v K h ∈ V h (K)
, on each interior edge e we have:

∂ n v h -i kv h = µ 1 e i k √ 2 2 s + µ 2 e -i k √ 2 2 s , ( 33
)
where s represents the curvilinear abscissa and µ 1 , µ 2 ∈ C. Consequently, the Lagrange multiplier is approximated in the following discrete space:

M h = µ h ∈ M; ∀K ∈ τ h , µ K h | e = µ K 1 e i k √ 2 2 x + µ K 2 e -i k √ 2 2 x if e -→ x , µ K h | e = µ K 1 e i k √ 2 2 y + µ K 2 e -i k √ 2 2 y if e -→ y , µ 1 , µ 2 ∈ C .
The spaces V h and M h dened above correspond to the so-called R-4-2 element in the nomenclature of DGM (see [5]). Note that for the DGM, considering two dofs per edge leads also to a complete approximation of the Lagrange multiplier.

We must point out that, unlike the DGM, the Lagrange multipliers in mDGM are not continuous across interior boundaries. This is why on an edge shared by two elements, the Lagrange multipliers on each side of the edge are dierent. Consequently, the stencil of the matrix given by ( 19) is equal to 40 (see Figure 2). Note that in the DGM the Lagrange multipliers are equal on both sides of the edge, which leads to a stencil equal to 14.

The rst experiments consist in comparing the error delivered by both numerical methods (DGM and mDGM) for dierent values of ka, while maintaining kh constant. More specically, we consider ka = 10, 20, 30 and we choose the step size of the mesh discretization h/a such that kh = 1 5 , which is about 30 elements per wavelength. The results are depicted in Figure 3 and Figure 4. These results indicate the following:

The two methods deliver results with the same level of accuracy, as indicated in Figure 3: both curves are superposed.

As expected, the relative error is π/2 periodic (see Figure 3). On each period [(l-1)π/2, lπ/2] (with l = 1, 2, 3, 4), the error is symmetric with respect to the propagation angle (2l -1)π/4. Moreover, the error is minimal for π/4, 3π/4, 5π/4, 7π/4 (less than 1%) and maximal for 0, π/2, π, 3π/2 (about 4%). This is due to the chosen basis, which includes the exact solution when the propagation angle is (2l -1)π/4, with l = 1, 2, 3, 4 and to the fact that the Lagrange multiplier eld contains all the functions needed to have a complete approximation.

Figure 4 indicates that the R-4-2 element (for both methods) exhibits little pollution: increasing ka, while maintaining kh constant, leads to an increase in the relative error which is less than 0.5% at most (see Figure 4 at angles θ = lπ/2, with l = 0, 1, ..., 8).

Next, we compare the sensitivity of the total relative error (the mean value over the propagation angles) to the mesh size. The result depicted in Figure 5 is obtained for ka = 1. One can observe the following:

For h/a > 1 100 , the errors delivered by the two methods are comparable. The two curves are on top of each other.

For h/a < 1 100 mDGM outperforms DGM. As we rene the mesh (h/a < 1 100 ), DGM becomes unstable. Indeed, there is a dramatic loss in the accuracy of more than one order of magnitude. The error jumps from 0.09% (for h/a = 1 100 )to 1.5% (for h/a = 1 190 ). The instability observed in DGM seems to be related to the severe ill conditioning of the local matrices. Observe that mDGM remains stable as we rene the mesh. The last point of the curve was obtained for h/a = 1 450 , the limit of our computing platform. The total relative error for this mesh size is 0.04%.

We must point out that the performance of mDGM in this case is not sensitive to the choice of the approach for solving the linear system in Step 2. Both approaches deliver results with the same level of accuracy.

Eight plane waves per element

We approximate the primal variable using eight plane waves, positioned at: This choice corresponds to the following discrete space for the primal variable:

θ p = (p -1)π/4, ∀1 ≤ p ≤ 8.
V h (K) = v h | K = 1≤p≤8 e i kθp•x u p , θ p = t [cos θ p , sin θ p ] , θ p = (p -1)π/4, 1 ≤ p ≤ 8, u p ∈ C .
For an element v K h ∈ V h (K), the full approximation of ∂ n v K h -i kv K h leads to ve dofs per edge. More specically in mDGM, the discrete space M h corresponding to the full approximation of the Lagrange multiplier is given by:

INRIA M h = µ h ∈ M; ∀K ∈ τ h , µ K h | e = µ K 1 + µ K 2 e i kx + µ K 3 e -i kx + µ K 4 e i k √ 2 2 x +µ K 5 e -i k √ 2 2 x if e -→ x , µ K h | e = µ K 1 + µ K 2 e i ky + µ K 3 e -i ky + µ K 4 e i k √ 2 2 y + µ K 5 e -i k √ 2 2 y if e -→ y , µ 1 , µ 2 , µ 3 , µ 4 , µ 5 ∈ C .
Following the nomenclature introduced in [5], such an approximation is called the R-8-5 element. Note that mDGM can be implemented using less dofs per edge for the Lagrange multiplier. We recall that in DGM, the maximum number of dofs considered on each edge is three. Indeed, the computation of the normal derivative of the numerical solution leads to the following complete approximation:

λ h = µ 1 + µ 2 e i k √ 2 2 s + µ 3 e -i k √ 2 2 s ,
where s represents the curvilinear abscissa. This choice of approximation corresponds to the so-called R-8-3 element. We rst present the results obtained in the case of Approach 1, than the ones obtained in the case of Approach 2.

Performance assessment in the case of Approach 1

Since the full approximation in DGM requires three dofs per edge, rst compare the performance of mDGM and DGM when using the R-8-3 element. The result depicted in Figure 6 compares the relative error delivered by both methods, as a function of the propagation angle. This result is obtained for ka = 10 and h/a = 1/20, that is kh = 1 2 , corresponding to about 12 elements per wavelength. It shows a clear superiority of mDGM over DGM. In addition, we have: DGM delivers the exact solution for lπ/4, with l = 0, 1, ..., 8. Note that in each of these cases, the exact solution is represented by one of the basis functions of the considered element and all the functions obtained when computing the normal derivative are in the Lagrange multiplier eld. On the other hand, mDGM computes exactly the solution at angles π/4, 3π/4, 5π/4, 7π/4 only, as it will be shown later (see Figure 8). This is not surprising since two dofs are removed from the full approximation of the Lagrange multiplier and therefore, for the plane waves propagating in the directions parallel to the axis, the approximation is not complete.

The total relative error is about 0.091% for mDGM and is about 5% for DGM. This means that mDGM improves the accuracy by about one and a half order of magnitude.

Observe that DGM R-8-3 is an unstable element. Indeed, the error obtained for θ = (2l -1)π/8 should be the same for all l = 1, 2, ..., 8 and symmetric with respect to 0, π/4, 2π/4, ..., 8π/4. Figure 6 shows that these values are not equal. The instability observed in the DGM approach for R-8-3 element is not surprising since this element does not satisfy the numerical inf-sup condition required by DGM [11]. Almost everywhere in the mesh, for an element K there are twelve dofs for the Lagrange multiplier and only eight dofs for the primal variable. A dof must be removed from each edge. For this reason, two discrete spaces were suggested in [5] for the discrete dual variable, leading respectively to the so-called R-8-2a and R-8-2b elements. Since in the cited paper, the R-8-2b element was shown to deliver more accurate results than the R-8-2a element, we have compared the two methods when employing the R-8-2b element. We

λ h = µ 1 e i k √ 2 4 s + µ 2 e -i k √ 2 4 s , INRIA
where s is the curvilinear abscissa. Similarly to the previous numerical experiment, we have set ka = 10 and h/a = 1/20, which corresponds to kh = 1 2 . The result depicted in Figure 7 suggests the following:

As expected, both methods preserve the symmetry of the error with respect to the propagation angles θ = π/4, 3π/4, 5π/4, 7π/4.

The total relative error obtained with the mDGM is about 5% while the one obtained with DGM is about 0.099%. This superiority of DGM over mDGM is most likely due to the poor approximation of the Lagrange multiplier in the mDGM (three out of ve dofs are neglected), compared to the DGM, where only one dof out of three is neglected.

Next, we enrich the approximation of the Lagrange multiplier in mDGM. We use the elements R-8-3 and R-8-5 and compare mDGM to DGM, equipped with the best element R-8-2b. The results reported in Figure 8 are obtained for ka = 10 and h/a = 1/20. One can observe the following:

There is a little improvement in the accuracy of the results delivered by R-8-3 and R-8-5 mDGM elements over DGM R-8-2b element. The total errors obtained with these elements are 0.091% and 0.048% respectively, whereas the one delivered by DGM equipped with the best element, R-8-2b, is 0.099%. The R-8-5 element seems to be unstable, as illustrated by the loss of the symmetry at the propagation angles θ = lπ/4, with l = 0, 1, 2, ..., 7, corresponding to the basis functions. This numerical instability is also noticeable when we compare R-8-5 to R-8-3: for some propagation angles, R-8-3 is more accurate than R-8-5, which is contrary to what is expected since using ve dofs per edge leads to the full approximation of the Lagrange multiplier. We believe that the instability of the R-8-5 element is due to the loss of the linear independence of the shape functions, but also to the fact that the matrix of the system is positive semi-denite.

Next, we investigate the sensitivity of the total relative error with respect to h, the step size of the mesh discretization. We consider the case of R-8-2b element and we set ka = 1. We evaluate the total relative error, as well as the smallest eigenvalue of the local system. The results are reported in Table I. The results reveal that in the DGM approach, the error decreases as long as h/a > 1 6 . Then, the error jumps from 0.01% to about 33% for kh = 1 23 . This is not surprising. Indeed, the local systems in the DGM are nearly singular and therefore extremely ill-conditioned when h becomes small, as indicated by the values corresponding to the smallest eigenvalues.

The smallest error delivered by mDGM, which is about 0.01%, is obtained for kh = 1 8 . Then the error jumps to 37% for kh = 1 23 . This instability is, a priori, unexpected and very surprising since we have introduced the Robin-type boundary condition to address specically this issue, as demonstrated in the case of the R-4-2 element. A quick look at Table I indicates that the local system corresponding to mDGM becomes nearly singular too. Hence, contrary to our goal, the presence of α seems to be not sucient to avoid the singularity of the local systems. We believe that the singularity of the local system in the mDGM formulation is due to the loss of the linear independence of the shape functions (eight plane waves) as h becomes small, as well as to the non-negativeness nature of the matrix A of Approach 1. The loss of the linear independence can be demonstrated as follows: let K be an element of the mesh. For simplicity, we assume K to be the square

[0, h] × [0, h]. Then, for a function w K h ∈ V h (K), there exist c 1 , c 2 , ..., c 8 ∈ C such that: w K h = 8 m=1 c m e i kθm•x .
Assume that:

a K w K h , w K h = 0. (34)
INRIA Consequently, w K h = 0 on ∂K, which means:

8 m=1
c m e i kθm•x = 0 on ∂K.

We write this equality on each of the four edges of K. For all x ∈ [0, h] and y ∈ [0, h], we have:

                         c 1 + c 5 + (c 2 + c 4 )e i k √ 2 2 y + c 3 e i ky + (c 6 + c 8 )e -i k √ 2 2 y + c 7 e -i ky = 0, c 3 + c 7 + (c 2 + c 8 )e i k √ 2 2 x + c 1 e i kx + (c 4 + c 6 )e -i k √ 2 2 x + c 5 e -i kx = 0, c 1 e i kh + c 5 e -i kh + (c 2 e i k √ 2 2 h + c 4 e -i k √ 2 2 h )e i k √ 2 2 y + c 3 e i ky + (c 6 e -i k √ 2 2 h + c 8 e i k √ 2 2 h )e -i k √ 2 2 y + c 7 e -i ky = 0, c 3 e i kh + c 7 e -i kh + (c 2 e i k √ 2 2 h + c 8 e -i k √ 2 2 h )e i k √ 2 2 x + c 1 e i kx + (c 6 e -i k √ 2 2 h + c 4 e i k √ 2 2 h )e -i k √ 2 2 x + c 5 e -i kx = 0.
Therefore, we deduce that:

             c 1 = c 3 = c 5 = c 7 = 0 c 4 = -c 2 c 6 = c 2 c 8 = -c 2 c 2 e i k √ 2 2 h -e -i k √ 2 2 h = 0. (35)
The problem here is that when h → 0, we have:

e i k √ 2 2 h → 1 and e -i k √ 2
2 h → 1 and hence, numerically speaking, it is not necessary to have c 2 = 0 in order to have w K h = 0 on ∂K. Consequently, w K h may not be equal to 0 when (34) is satised. This computation shows that when h tends to 0, the eight plane waves become linearly dependent, which leads to the singularity in the local matrix. Remark 4. We have performed additional numerical experiments for higher frequencies and observed that both methods become unstable as we rene the mesh. DGM is stable as long as kh > 1 6 , while mDGM seems to remain stable longer (kh > 1 9 ). We must point out that the source for the instability is however dierent. DGM is unstable not only because of the singularity of the local systems, but also due to the loss of the linear independence of the plane waves as h becomes small. We believe that the resulting local system in DGM is also very sensitive to this loss. On the other hand, mDGM restores the stability of the local problems and leads to a more stable formulation when the shape functions remain linearly independent, as it has been demonstrated in the case of the R-4-2 element (see Figure 5) and in the next experiments. However, the previous numerical results suggest that the loss of the linear independence aects dramatically the stability of mDGM when using Approach 1, that is when the linear system is non-negative only.

The following experiment reveals the behavior of mDGM and DGM when the used shape functions remain linearly independent as we rene the mesh. This experiment consists in approximating the solution, at the element level, using seven plane waves, positioned at:

θ p = 2(p -1)π/7, 1 ≤ p ≤ 7.
(36)

We have maintained the same two dofs per edge as in R-8-2b (see Section 5.2). Following the nomenclature introduced in [5], we will refer to this element as R-7-2. For ka = 1, we have compared the sensitivity of the total relative error to the mesh size for the mDGM R-8-2b and R-7-2 elements. The result depicted in Figure 9 illustrates the following: The error delivered by the R-7-2 element decreases as long as kh > 1 18 , unlike the R-8-2b element, which delivers the smallest error for h/a = 1 9 . This means that there is a reduction of factor 2 on the mesh size, while maintaining the stability.

The R-7-2 element is shown to be more accurate. For each mesh size, the error delivered by this element is smaller than the one obtained with the R-8-2b element. Moreover, the most accurate approximation (about 0.001% for h/a = 1 18 ) obtained with R-7-2 is one order of magnitude lower than R-8-2b element (about 0.01% for h/a = 1 9 ). Last, note that the R-7-2 element remains stable while rening the mesh. For h/a = 1 100 , the total relative error is about 1%, unlike the R-8-2b element in which the error jumps from 0.01% (for h/a = 1 9 ) to 59% (for h/a = 1 24 ). We have also compared the errors delivered by the mDGM R-7-2 element to the ones obtained with the DGM R-7-2 element. The result is reported in Figure 10. The following observations are noteworthy:

The accuracy of the two methods is comparable for h/a > 1 8 . In this region the two curves are superposed. The DGM R-7-2 element delivers the most accurate approximation (which is about 0.005%) for h/a = 1 12 . Observe that mDGM becomes unstable slightly later: the smallest error (about 0.001%) is obtained for h/a = 1 18 ). Although for both methods we observe numerical instabilities as soon as h/a < 1 18 , mDGM is more accurate than DGM. For any mesh size the error delivered mDGM is smaller than the one obtained with DGM. Moreover, for some mesh sizes, mDGM outperforms DGM by one order of magnitude. We believe that this is due to the local problems which are nearly singular in DGM.

The results depicted in Figure 9 seem to be surprising since one may expect that approximating the primal variable with eight plane waves leads to more accurate results than when using seven plane waves. Here, it seems that the linear independence when using seven plane waves is less sensitive (remains longer) to the mesh size h/a. Consequently, the linear system corresponding to Approach 1 is more stable in this case than when using eight plane waves. In Table II we report the total relative error, as well as the smallest eigenvalue of the local system for the R-7-2 element for both DGM and mDGM methods, when ka = 1.

As in the case of the R-8-2b element, the eigenvalues of the local matrices corresponding to each mesh size have the same real part in DGM and mDGM. The Robin-type condition used in mDGM leads to more important imaginary parts in mDGM.

The values of the smallest eigenvalues and the comparison to the ones reported in Table I show that the linearly independence of the seven plane waves is less sensitive to the mesh renement. Indeed, for example for h/a = 1 9 the real part of the smallest eigenvalue obtained when using eight plane waves at the element level is 8.64 • 10 -13 . This is four orders of magnitude larger than the real part of the smallest eigenvalue of the matrix obtained with seven plane waves (4.00 • 10 -9 ). A quick comparaison of the errors reported in Tables I and II (see also Figure 9 for mDGM) shows that in both formulations a better conditioning of the local matrices leads to more accurate approximations. Moreover, in both methods there is a reduction of factor 2 on the mesh size, while maintaining the stability. These two important points are related to the fact that the seven shape functions remain linearly independent as we rene the mesh.

As it was shown in Figure 10 and reported in Table II, the numerical instabilities appear earlier in DGM than in mDGM. We believe that this is due to the local systems which become nearly singular with the mesh rening. In mDGM the observed numerical instabilities are due to the non-negativeness nature of the global matrix corresponding to Approach 1.

Last, we must point out the fact that the seven shape functions are becoming linearly dependent with the mesh renement, but more slowly than the eight plane waves corresponding to the R-8-2b element. This behavior of the shape functions is predictible when observing the dependence of the smallest eigenvalue of the local matrix with respect to the mesh size. Consequently, this element will be ultimately unstable. Note that the numerical stability occurs only in the case of mDGM when equipped with Approach 2. Indeed, one can observe that the errors (see dashed curves in Figure 13) are oscillating as a/h is increasing, but their magnitude remains steadily about 1%, whereas the magnitude of the errors delivered by the two other methods increase to over 100%.

Summary and conclusions

We have designed a new solution methodology, called mDGM, for Helmholtz problems which is easy to understand and implement. At the element level, we approximate the solution by a superposition of plane waves. Consequently, the obtained solution is discontinuous and Lagrange multipliers are introduced to ensure the continuity in a weak sense. Unlike the DGM, the Lagrange multiplier is also discontinuous, which allows us to consider well-posed local problems. The algebraic approach requires solving local linear systems with multiple right-hand side: the system's size is given by the number of plane waves considered in the local basis. These problems are independent from one element to another and therefore can be solved in paraller. The global system, whose size is the number of total dofs used for approximating the Lagrange multiplier, can be either INRIA positive semi-denite or denite depending on the approach adopted at the continuous level. The numerical results we have presented show that the proposed method is more stable than the DGM. When using Approach 2, mDGM is not only more stable than DGM, but also exhibits a better level of accuracy. More specically, as indicated by the reported numerical results, mDGM reduces the level of errors by one to two orders of magnitude depending on the mesh size and on the element.

Figure 1 : 4 INRIAA

 14 Figure 1: Four plane waves positioned at θ p = π/4 + (p -1)π/2, ∀ 1 ≤ p ≤ 4

Figure 2 :

 2 Figure 2: A 40 dofs stencil corresponding to 2 dofs per edge of∂ n v K h -i kv K h on the edges of the mesh, for v K h ∈ V h (K). Observe that, ∀v K h ∈ V h (K), on each interior edge e we have:

Figure 3 :

 3 Figure 3: Performance of the two methods for kh=1/5

Figure 4 :Figure 5 :

 45 Figure 4: Pollution eect for the R-4-2 element

Figure 6 :

 6 Figure 6: Performance of DGM and mDGM equipped with the R-8-3 element for ka=10, h/a=1/20

Figure 7 :

 7 Figure 7: Performance of DGM and mDGM equipped with the R-8-2b element for ka=10, h/a=1/20

Figure 8 :

 8 Figure 8: Performance of DGM equipped with the R-8-2b and of mDGM equipped with the R-8-3 and R-8-5 for ka=10, h/a=1/20

Figure 9 :

 9 Figure 9: Convergence rates for mDGM R-8-2b and R-7-2 elements

Figure 10 :

 10 Figure 10: Convergence rates for mDGM R-7-2 and DGM R-7-2 elements

5. 2 . 2

 22 Performance assessment in the case of Approach 2 We have performed several numerical experiments to assess the performace of INRIA mDGM when at Step 2 the linear system is given by Approach 2 (see Eqs.

Figure 11 :Figure 12 :Figure 13 :

 111213 Figure 11: Performance of the three methods, kh=1/2

Table 1 :

 1 Performance • 10 -10 + 2.9 • 10 -18 i 0.054% 5.7 • 10 -10 -3.1 • 10 -11 i 1/5 0.011% 9.5 • 10 -11 + 8.0 • 10 -19 i 0.036% 9.5 • 10 -11 -4.1 • 10 -12 i 1/6 0.019% 2.2 • 10 -11 -2.4 • 10 -18 i 0.025% 2.2 • 10 -11 -7.9 • 10 -13 i 1/7 0.092% 6.5 • 10 -12 + 5.9 • 10 -19 i 0.019% 6.4 • 10 -12 -2.0 • 10 -13 i 1/8 0.394% 2.2 • 10 -12 + 1.8 • 10 -18 i 0.015% 2.2 • 10 -12 -5.9 • 10 -14 i 1/9 1.206% 8.6 • 10 -13 -1.6 • 10 -18 i 0.017% 8.6 • 10 -13 -2.1 • 10 -14 i 1/20 4.734% 1.5 • 10 -15 -6.0 • 10 -21 i 3.590% 1.5 • 10 -15 + 5.9 • 10 -17 i 1/21 11.664% 9.8 • 10 -16 + 7.0 • 10 -18 i 5.889% 9.6 • 10 -16 + 2.4 • 10 -17 i 1/23 33.435% 4.7 • 10 -16 -3.9 • 10 -18 i 37.378% 4.7 • 10 -16 -6.7 • 10 -17 i

			of the R-8-2b element for DGM
		and mDGM when ka=1
			DGM		mDGM
	h/a	Total	The	Total	The
		relative	smallest	relative	smallest
		error	eigenvalue	error	eigenvalue
	1/4	0.016% 5.7

Table 2 :

 2 Performance • 10 -09 + 2.5 • 10 -17 i 0.005% 4.0 • 10 -09 -1.1 • 10 -10 i 1/12 0.005% 7.1 • 10 -10 + 6.0 • 10 -18 i 0.003% 7.1 • 10 -10 -1.5 • 10 -11 i 1/18 0.035% 6.3 • 10 -11 + 7.7 • 10 -18 i 0.001% 6.3 • 10 -11 -8.5 • 10 -13 i 1/20 0.222% 3.3 • 10 -11 + 3.1 • 10 -19 i 0.010% 3.3 • 10 -11 -4.1 • 10 -13 i 1/25 1.075% 8.7 • 10 -12 -2.0 • 10 -18 i 0.052% 8.7 • 10 -12 + 2.4 • 10 -17 i 1/34 0.589% 1.4 • 10 • 10 -17 i 1/100 9.266% 2.1 • 10 -15 + 4.5 • 10 -18 i 1.390% 2.1 • 10 -15 + 2.5 • 10 -18 i 1/150 243.2% 1.9 • 10 -16 + 1.5 • 10 -19 i 67.75% 1.8 • 10 -16 -8.3 • 10 -18 i

			of the R-7-2 element for DGM
		and mDGM when ka=1
			DGM		mDGM
	h/a	Total	The	Total	The
		relative	smallest	relative	smallest
		error	eigenvalue	error	eigenvalue
	1/4	0.020% 5.2 • 10 -07 -3.0 • 10 -17 i 0.030% 5.2 • 10 -07 -3.2 • 10 -08 i
	1/5	0.014% 1.4 • 10 -07 -4.8 • 10 -17 i 0.019% 1.4 • 10 -07 -6.6 • 10 -09 i
	1/6	0.010% 4.6 • 10 -08 -8.0 • 10 -18 i 0.012% 4.6 • 10 -08 -1.9 • 10 -09 i
	1/7	0.008% 1.8 • 10 -08 + 2.0 • 10 -17 i 0.009% 1.8 • 10 -08 -6.3 • 10 -10 i
	1/8	0.007% 8.1 • 10 -09 + 3.9 • 10 -17 i 0.007% 8.1 • 10 -09 -2.5 • 10 -10 i
	1/9	0.006% 4.0

-12 + 4.4 • 10 -18 i 0.099% 1.4 • 10 -12 -9.9 • 10 -15 i 1/50 0.493% 1.4 • 10 -13 -1.1 • 10 -19 i 0.285% 1.4 • 10 -13 -6.8 • 10 -16 i 1/75 2.046% 1.2 • 10 -14 + 4.0 • 10 -18 i 0.373% 1.2 • 10 -14 -5.1

RR n°7050

The continuous approach 5

http://www.inria.fr ISSN 0249-6399

Acknowledgements

The authors acknowledge the support by TOTAL and INRIA/CSUN Associate Team Magic, INRIA Bordeaux Sud-Ouest Center. Any opinions, ndings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reect the views of TOTAL, INRIA or CSUN.


A Reineix

O Maurice

Optimizations of systems using the tensorial analysis of networks for electromagnetic compatibility

Introduction

The tensorial analysis of network (TAN) developed by G.Kron in 1939 has today proved its capability to give answers to the difficult problem of the electromagnetic compatibility (EMC) prediction of complex systems [START_REF] Leman | Use of the circuit approach to solve large EMC problems[END_REF]. The purpose of this article is to show how this technique can allow to optimize the choice of parameters in system functions. We recall briefly the principles of the TAN, before to speak of the hybrid approach. Using this last technique, we detail a very simple example to illustrate the mechanisms involved. After what we can speak about the optimization of the system studied through a classical experience plan. The techniques can be then used for many more complex systems and efficient mathematical techniques of optimizations.

Recall on the tensorial analysis of network

The tensorial analysis of networks (TAN) developed by G.Kron benefits today of many publications describing its application and capabilities [START_REF] Leman | Use of the circuit approach to solve large EMC problems[END_REF][START_REF] Maurice | Electromagnetic compatibility of complex systems[END_REF][START_REF] Kaufmann | Equations and electrical network solving in the time domain using the tensorial method[END_REF]. G.Kron has understood all the powerfull mechanism of the tensorial algebra to resolve complex problems (at his time about electrical machines, quantum mechanics, computing, etc.). The reader can see the web description of the Kron's work in "http:\\www.quantum-chemistryhistory.com\Kron Dat\KronGabriel1.htm". Until today the Kron's method was used by many various people in many different jobs. More than to be a method to solve the problems, the Kron's formalism is first of all a theoretical technique to establish its equations even in multiphysic [START_REF] Hubert | Combined effects induced by electromagnetic pulses and single event effect on electronics[END_REF]. The Kron's idea was to reuse as far as possible all the work made previously in modelling by an engineer in any new problem he had to solve. To do that, he shows that the TAN gives all the mathematical tools needed to connect old networks in a new one, or to make the inverse, i.e. to cut a complex problem in littleler ones, without missing any information from the global network. Kron has called this second possibility "Diakoptic" [START_REF] Stenbakken | Diakoptic and large change sensitivity[END_REF]. The Kron's method is based on the use of two dual spaces: the voltage space and the current one, each of them defined in R n . Whatever is the representation of a network (using edges, vertices, nodes, meshes, ...), Kron had remarked that its total energy should not change. So, the global power of any network can be seen as an invariant.

Once two dual spaces and an invariant are defined, next step is to define a metric [START_REF] Hoffmann | Kron's Non-Riemannian Electrodynamics[END_REF]. In case of electrical circuit, this metric is not a Riemannian one. But its stills possible to define it through operators, and the result is a tensorial algebra to theoretically study any complex networks [START_REF] Kaufmann | Equations and electrical network solving in the time domain using the tensorial method[END_REF]. The technique used is so a real topological one as deeply shown by Branin [START_REF] Branin | The algebraic topological basis for network analogies and the vector calculus[END_REF]. In all the article we use the index notation defined by Einstein [START_REF] Einstein | Four Princeton conferences[END_REF]. On the principle, once defined the covector of the effort (voltage in an electrical circuit) e µ and the flux (current) vector f ν . The invariant W results from their contracted product W = e µ f µ . The metric can be deduced from the same invariant writing: W = z µν f µ f ν . Classically, the electrical circuit solvers work in the nodal description. They often use the Modified Nodal Analysis which includes the voltage sources, but it doesn't change the fundamental concepts 1 . It means they translate the impedance matrix in the nodal space using the incidence matrix [START_REF] Peikari | Fundamentals of network analysis and synthesis[END_REF]. This matrix can be seen as a connectivity between the voltage differences V defined in the edge space, and the electrical potentials ψ defined in the nodes space. Starting from the gradient Maxwell's law we can simply write: V a = B p a ψ p . Kron has shown that this kind of translation could be used between the currents f defined in the edge space and the meshes currents m. To do that, he define a connectivity C writing: f a = C a µ m µ . This transformation defined, it began possible to compute all the problem in the mesh description, changing the metric through:

z µν = C a µ C b ν z ab .
What is remarkable, is the fact that in this space, the edge voltage vanishes. Effectively, as the Maxwell's equation tell us, Rot(GradV) = 0. The whole problem is simply compacted in a unique equation: e µ = z µν m ν . But for various 1 More, our hybrid formalism get out this problem as we will show it in the example. reasons, it can be interesting to keep a nodal description, at least to reuse cases already computed under this space for example. To allow this mixted approach, the authors have proposed to create mix tensors [START_REF] Maurice | Hybrid tensorial analysis as theoretical tool for EMC and electromagnetism[END_REF]. This mechanism give all the agility needed to manage the complexity and the diversity in systems.

A simple example as illustration

We consider the system described fig. 1. All the next equations used are extracted from Paul's formalism on lines [START_REF] Paul | Sinusoidal excitation of a line[END_REF].

The current in the input of the line can be defined depending on the incident voltage wave and the coefficient reflexion Γ L :

i(0) = V + Z c 1 -Γ L e -2jβL (1)
But i(0) is a proportion of the current source i 0 :

i(0) = i 0 R 0 R 0 + Z in = αi 0 ( 2
)
V + is defined as a function of the voltage at the input of the line:

V + = V (0) 1 + Γ L e -2jβL (3)
Combining the previous equations we obtain:

V (0) = Z c αi 0 1 + Γ L e -2jβL 1 -Γ L e -2jβL (4)
On the other side of the line we can define the relation between the transmitted voltage and the input voltage and as a consequence, the relation between the transmitted current and the input voltage:

V (L) = V (0)e -jβL 1+Γ L 1+Γ L e -2jβL . . . . . . ⇒ i(L) = 1 R L V (L) (5)
Writing:

Y = 1 αZ c 1 -Γ L e -2jβL 1 + Γ L e -2jβL (6)
The nodal relation is synthetized in i 0 = Y V (0). On the second network, the mesh relation is:

e 0 -jγωi(L) = (Z 1 + Z 2 ) J = ZJ (7)
Where γ is the mutual inductance coefficient.

It is there assumed that the back interaction from the second network to the first one is negligeable. Using the relation between i(L) and V (0) we can write: i(L) = gV (O). Finally the whole system can be written:

i 0 e 0 = Y 0 jγgω Z V (0) J (8)
This is the basic principle of hybrid tensor applied to connected networks. This can be applied to large networks where the previous elements become vectors and the metric is hybrid with both impedance, admittance metrics and matrix too making the connection between the element of same physical dimension. The previous equation becomes:

i a e µ = Y ab M a ν M b µ Z µν V b J ν (9)
4 System optimization Writing ζ = jγgω, λ = Y -1 and S = Z -1 , the solution of the previous system is given by:

V (0) J = λ 0 -ζλ S i 0 e 0 (10)
which can be compacted in the form: f ν = Ω ν µ h µ . Now we wonder how J depends on γi 0 or e 0 in frequency, because we hope to keep the coupling influence (γi 0 ) inferior to the self generator e 0 . We can theoretically calculate for the first request:

∂J ∂γ = ∂f 2 ∂γ = ∂ ∂γ (-ζλi 0 + Se 0 ) (11)
We found easily: ∂J ∂γ = |gωλ|. The second request is directly given by |S|. The first term gives:

∂J ∂γ = ZcR 0 (1+Γ L ) (Zc+R 0 )R L ω . . . . . . [1 + Γ 0 Γ L cos(2βL)] 2 + [Γ 0 Γ L sin(2βL)] 2 -1 (12)
With Γ 0 the coefficient of reflexion: Zc-R 0 Zc+R 0 . So, to keep the dependencies equals in a matched line case, it is necessary to have the admittance S of the form: R 1 + ω 0 ω 2 -1

.

Until the cutoff frequency ω 0 , the generated current will keep the same distance to the one induced by i 0 .

Conclusion

Through this simple example we try to illustrate the method and how it can comply with the various needs of any engineers. The optimization of the system can be studied numerically of course, through experience plans [START_REF] Goupy | Introduction of experience plans[END_REF]. Once an equation is available to represent the system, any analytical mathematical study can be imagine to realize a theoretical analysis of the system behavior. Several studies was already conducted using this technique in conjonction with uncertainties computations. Very good results were obtained showing that the Kron's method really gives a support for these kinds of complex problems, allowing all the formulations to take part of its formalism.

Figure 1 :

 1 Figure 1: System considered as illustration


Kamal Hammani

Benjamin Wetzel

Bertrand Kibler
email: bertrand.kibler@u-bourgogne.fr

Julien Fatome

Christophe Finot

Guy Millot

Nail Akhmediev

John M Dudley

Spectral dynamics of modulation instability described using Akhmediev breather theory

Keywords: Nonlinear optics, fibers), 060.5530 (Pulse propagation

A central process of nonlinear physics is modulation instability (MI), where periodic perturbations on a plane wave are amplified to generate a parametric cascade of spectral sidebands [1,2]. During the initial evolution, the sidebands experience exponential growth at the expense of the pump, but the subsequent dynamics is more complex and different scenarios of energy exchange between the spectral modes are possible. This dynamics is closely related to the universal Fermi-Pasta-Ulam (FPU) recurrence phenomenon [3].

In nonlinear fiber optics, MI has been the subject of extensive study using both numerical simulations as well as varying analytical approaches based on truncated sideband solutions to the nonlinear Schrödinger equation (NLSE) [4,5]. However, it has recently been shown that MI can be accurately described by a family of exact analytic solutions known as Akhmediev breathers (ABs) [6], and this has led to the design of experiments exciting the rational Peregrine soliton [7,8]. An additional important feature of the AB theory is that it has provided an analytical description of the spectral decay of the sideband modes at the peak of the MI gain, and this description has provided new insights into the development of the noise-driven supercontinuum generation [9,10].

In this Letter, we complete the analytic description of the AB spectral dynamics for the case of arbitrary gain, allowing the development of sideband generation in induced MI to be compared with theory for general experimental conditions. Experiments measuring pump and multiple sideband generation over a growth-return cycle of MI are used to quantitatively test the theoretical results. In fact, with the excellent signal to noise ratio in our experiments, we can compare experiment and theory out to more than 10 spectral sidebands and over 30 dB dynamic range. To our knowledge this is the most complete characterization of MI dynamics in any NLSE-governed system, and significantly extends previous studies in either optics [3,11] or hydrodynamics [12,13] where only a greatly reduced number of harmonics (~ 3) could be measured.

Our analysis is based on the normalized NLSE:

ψ ξ + 1/2 β 2 ψ ττ + |ψ| 2 ψ = 0 where ψ (ξ,τ) is the field envelope.
The AB solution describes growth and decay of a harmonically-perturbed plane wave, and its temporal description in terms of localized compressed pulses is now well-known [6,7]. In addition to the time-domain description, expansion in a Fourier series and appropriate integration yields exact solutions for the pump and spectral harmonic amplitudes. Labelling the amplitudes of the pump and the n th sideband A 0 and A n respectively (n = ±1, ±2, etc), and ignoring factors of constant amplitude and phase, we derive the following results describing pump and sideband evolution with distance:

2 0 2 2 2 2 | | sinh cosh ( ) 1 (1a) cosh 2 cosh cosh 2 sinh cosh ( ) (1b) 2 cosh 2 n n ib b p b A b a b b a ib b p b A a b a ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ + = - -   - - + =   -  
The parameters a and b are the normalized perturbation frequency and instability growth rate, and in dimensional units (for dispersion β 2 and nonlinearity γ)

we have 2a = [1 -(ω mod / ω c ) 2 ] and b = [8a (1-2a)] 1/2 .
Here ω mod is the dimensional modulation frequency in the range 0< ω mod < ω c where ω c 2 = 4γ P 0 / |β 2 | for an initial plane wave power P 0 . The parameter p = 2 1/2 ω mod /ω c. The normalized distance is related to dimensional distance z via ξ = z / L NL with the nonlinear length L NL = (γP 0 ) -1 . We stress that this general result exactly predicts the dynamics of an arbitrary number of sidebands without any assumption of an undepleted pump [3][4][5].

The predictions of Eq (1) have been compared with simulations and experiment by considering induced MI in standard fiber around 1550 nm. Fig. 1(a) shows our experimental setup. A laser diode at 1550.08 nm passes through a 16 GHz intensity modulator (IM) to generate two equal amplitude sidebands at -15.6 dB relative to the pump. This weak modulation is that where the AB theory is expected to well-describe the propagation dynamics and a phase modulator (PM) is used for linewidth broadening to suppress Brillouin scattering. The power after modulation was adjusted using a low noise erbium doped fiber amplifier (EDFA). The use of standard SMF-28 fiber in our experiments allows fusion splicing to be performed with negligible loss, so that we could readily study the spectral evolution dynamics using cutback experiments, and so that we could also repeat experiments for a range of different conditions. The fiber parameters were β 2 = -21.4 ps 2 km -1 , β 3 = 0.012 ps 3 km -1 , γ = 1.2 W -1 km -1 and the linear loss was 0.19 dB km -1 . Here P 0 = 0.46 W so that the normalised modulation frequency a = 0.45 such that we are far from the MI gain peak and we perform a valid test of the generalised AB spectral result in Eq. (1). Fig. 1(b) shows typical results using a fiber length of 5 km and with the input power adjusted to generate a maximallybroadened spectrum. With the low noise levels in our experiment, we could resolve more than 14 sidebands on either side of the pump (solid line). The experimental results are compared with the predictions of Eq. (1) (the crosses) and we see very good agreement in predicting spectral decay over more than four orders of magnitude. The figure also shows results of numerical simulations which are also in good agreement. The simulations used an extended NLSE including third-order dispersion and Raman scattering, but we found that these additional effects have little effect at the point of maximum spectral broadening.

For the same parameters as above (16 GHz and a = 0.45) we then performed experiments to examine the evolution dynamics more fully. Specifically, the fiber length was changed over the range 0-8.5 km spectral measurements were made at ~250 m intervals. Fig. 2(a) plots these measurements. For visual clarity we plot only 9 sidebands and use a plotting floor at -30 dB, although the dynamic range of our measurements extended further as discussed above and shown in Fig. 1(b). These results are significant in clearly showing the full instability dynamics of multiple sideband generation, as well as the recovery back towards the initial state with propagation beyond 5 km. The experimental results in Fig. 2(a) were compared with both numerical simulations and the theoretical predictions from the spectral AB theory in Eq. ( 1). The simulation results (again using an extended NLSE) explicitly show that although propagation beyond 5 km sees the expected spectral compression as the sideband energy begins to be converted back to the pump, ideal FPU recurrence is not observed. Rather, after around 7 km, we see the onset of more complex dynamics. In this context, we note that we have been able to perform simulations with and without the higher-order terms in the NLSE, and this has allowed us to clarify that loss, third order dispersion and the Raman effect in fact have no influence on this deviation from ideal propagation behaviour. Rather, it is the non-ideal initial conditions (that do not correspond to the mathematical ideal case for a single AB cycle) that induce the deviation from recovery to the initial state.

A full quantitative comparison between experiment and the predictions of Eq. ( 1) is shown in Fig. 3. Here, we extract the experimental pump and sideband intensity for various sideband order n (circles) and compare the evolution as a function of distance with corresponding results from simulations (blue short dashes) and the AB spectral theory (black solid line). The spectral normalisation in all cases is relative to the total initial energy.

In addition for completeness, for the pump and n=1 sideband, the red long dashes show the predictions of the well-known truncated three wave (TW) model [4,5]. It is important to note that no free parameters were used in the comparison between experiment and simulation or the three wave theory. Also, only the distance offset for the AB theoretical curve was adjusted to match the position of maximum pump depletion; once this was fixed, no relative normalisation of pump or sideband intensity to experimental or simulation data was made.

The agreement between simulation and experiment is very good, but this of course might be expected given the fact that the simulation initial conditions were chosen to match exactly those used in experiment. We also note that the TW theory is accurate in predicting the distance to maximum compression and the extent of pump depletion. On the other hand, because it assumes that all the depleted energy is transferred only to the first sideband pair, the TW theory significantly overestimates the sideband amplitude.

What is particularly impressive, however, is the degree to which the AB theory describes both the dynamical evolution and the quantitative extent of pump depletion and sideband growth, even out to the n = 11 sideband. In fact, our ability to study dynamics to such a high sideband order reveal a previously uncommented feature of the dynamics where we see that the temporal compression is associated with an increased spatial localisation of the energy transfer to higher sideband orders.

Although the AB theory has been previously shown to describe field characteristics at the maximal spectral broadening, these results now confirm its validity in describing the dynamical evolution from growth to decay over an MI cycle. Realistic simulations using an extended NLSE identify that the primary physical factors that preclude perfect recurrence are due to imperfect initial conditions, specifically a non-zero relative phase between injected pump and sidebands. With suitable additional control of the injected field, we anticipate that creation of ideal nonlinear breather states should be possible. These results show how analytical breather solutions of the NLSE provide new insights into seemingly well-known processes such as fiber MI.

Fig. 1 .

 1 Fig. 1. (a) Experimental setup. Diagnostics recorded power and spectra. (b) Solid line shows measured spectral intensity at the point of maximum sideband generation. Circles and crosses show simulation and analytic results.

Fig. 2 .

 2 Fig. 2. Dynamics of multiple sideband generation. (a) Experimental results. (b) Results from numerical simulation using realistic NLSE model. The generation of multiple sidebands around 5 km and recovery back towards the initial state with further propagation is clearly seen

Fig. 3

 3 Fig. 3 Spectral dynamics for the pump and sideband orders as shown. Experiment (circles) is compared with simulation (blue dashed line), the exact AB result Eq. (1) (black line) and the truncated three wave (TW) model (red dot-dashed line).

We thank the French Agence Nationale de la Recherche projects MANUREVA ANR-08-SYSC-019 and IMFINI ANR-09-BLAN-0065, the Conseil Régional de Bourgogne, and the Australian Research Council Discovery Project DP110102068 for support.