Lucas Chesnel
email: lucas.chesnel@cmap.polytechnique.fr

Xavier Claeys
email: claeys@ann.jussieu.fr

A numerical approach for the Poisson equation in a planar domain with a small inclusion

Keywords: Small hole, asymptotic analysis, singular perturbation, finite element method

We consider the Poisson equation in a domain with a small hole of size δ. We present a simple numerical method, based on an asymptotic analysis, which allows to approximate robustly the far field of the solution as δ goes to zero without meshing the small hole. We prove the stability of the scheme and provide error estimates. We end the paper with numerical experiments illustrating the efficiency of the technique.

Introduction

In the present article, we consider ω, Ω ⊂ R 2 two bounded Lipschitz domains such that ω ⊂ Ω. Define ω δ := {x ∈ R 2 , x/δ ∈ ω} and Ω δ := Ω \ ω δ . Given a data f ∈ L 2 (Ω), we are interested in devising a robust and accurate numerical method for approximating, for small values of δ, the far field of the function satisfying

u δ ∈ H 1 0 (Ω δ ) and -∆u δ = f in Ω δ . (1)
In ( 1), H 1 0 (Ω δ ) denotes the subspace of the elements of the Sobolev space H 1 (Ω δ ) vanishing on ∂Ω δ . On the other hand, we call far field of u δ the restriction of u δ to Ω \ D r , where D r := D(0, r) is the disk with fixed arbitrary radius r > 0. Problem [START_REF] Antoine | Wide frequency band numerical approaches for multiple scattering problems by disks[END_REF], or variants of it, arises as a simple but relevant model in many applications ranging from electrical engineering [START_REF] Bérenger | A multiwire formalism for the FDTD method[END_REF][START_REF] Railton | The treatment of thin wires in the FDTD method using a weighted residuals approach[END_REF] to flow transport around wells [START_REF] Chen | Numerical homogenization of well singularities in the flow transport through heterogeneous porous media[END_REF][START_REF] Peaceman | Interpretation of well-block pressures in numerical reservoir simulations[END_REF]. This kind of problem also appears when considering wave scattering by small impenetrable inclusions [START_REF] Cassier | Multiple scattering of acoustic waves by small sound-soft obstacles in two dimensions: mathematical justification of the Foldy-Lax model[END_REF].

In order to solve numerically Problem (1), a crude but rather natural idea would consist in neglecting the influence of the small inclusion on the total field u δ . Indeed (see for example [START_REF] Maz'ya | Asymptotic theory of elliptic boundary value problems in singularly perturbed domains[END_REF]), as δ → 0 the function u δ converges toward u 0 the solution to the limit problem where the inclusion has disappeared u 0 ∈ H 1 0 (Ω) and -∆u 0 = f in Ω.

(2)

However, in the general case (more precisely, when u 0 (0) = 0), the convergence turns out to be very slow: for any arbitrary radius r > 0 such that D r ⊂ Ω, we have u δ -u 0 H 1 (Ω\Dr) ≥ C | ln δ| -1 f L 2 (Ω) , for some constant C > 0 independent of δ. To give an idea | ln δ| -1 ≈ 0.0434 for δ = 10 -10 . Thus, neglecting the presence of the small inclusion is not satisfactory from a computational point of view, and a reasonable numerical approach for (1) should reproduce accurately the perturbation induced by the presence of the small inclusion ω δ .

Most of the numerical approaches that could be considered for dealing with this problem suffer from a numerical locking effect [START_REF] Babuška | On locking and robustness in the finite element method[END_REF]: performances of standard strategies deteriorate as δ → 0. Admittedly, robust strategies already exist in the literature, like the multi-scale finite element method coupled with some mesh refinement strategy [START_REF] Efendiev | of Surveys and Tutorials in the Applied Mathematical Sciences[END_REF], or the boundary element method (see e.g. [START_REF] Antoine | Wide frequency band numerical approaches for multiple scattering problems by disks[END_REF]). These techniques provide satisfying results in many cases, but they require careful and thorough implementation efforts, and/or rely on strong assumptions such as homogeneity of the coefficients of the equation under consideration. Other numerical strategies are based on an approximation of u δ of the form "u 0 + corrector", where both terms of this sum are computed separately (see for example [START_REF] Dambrine | Influence of a boundary perforation on the Dirichlet energy[END_REF][START_REF] Bonnaillie-Noël | On moderately close inclusions for the Laplace equation[END_REF]). These approaches may induce substantial additional computational cost in a real life simulation. In many practical situations, small inclusions are not the main subject of concern, and it would be desirable to devise a simple, general purpose and implementation friendly method that would not rely on any kind of mesh refinement technique, while remaining robust as δ → 0. This is the purpose of the present article to describe and analyse a method matching these requirements, while relying on only one numerical resolution.

The outline of this article is the following. In Section 2, we summarize the main results concerning the asymptotic expansion of u δ with respect to the size of the small hole. Section 3 is dedicated to the construction of a model problem, based on the matched expansion of u δ , whose solution has the same far field asymptotics as u δ , up to a remainder in O(δ 1-), ∀ > 0. We prove this in Section 4 (see Proposition 4.2) and also show that consistency of any Galerkin discretization of this model problem is quasi-optimal and uniform with respect to δ. Finally, in Section 6, we present and comment numerical results that confirm and illustrate our theoretical conclusions.

Asymptotic expansion of the solution

Asymptotic analysis for problems involving small inclusions can be found in many works. We refer the reader to [START_REF] Campbell | Une justification de la méthode de raccordement des développements asymptotiques appliquée à un problème de plaque en flexion. Estimation de la matrice d'impédance[END_REF][START_REF] Il'in | Study of the asymptotic behavior of the solution of an elliptic boundary value problem in a domain with a small hole[END_REF][START_REF] Il'in | Matching of asymptotic expansions of solutions of boundary value problems[END_REF][START_REF] Maz'ya | Asymptotic expansions of eigenvalues of boundary value problems for the Laplace operator in domains with small openings[END_REF][START_REF] Maz | Asymptotic behavior of energy integrals under small perturbations of the boundary near corner and conic points[END_REF][START_REF] Maz'ya | Asymptotic theory of elliptic boundary value problems in singularly perturbed domains[END_REF][START_REF] Nazarov | Asymptotic conditions at a point, selfadjoint extensions of operators, and the method of matched asymptotic expansions[END_REF][START_REF] Nazarov | Asymptotic analysis of shape functionals[END_REF][START_REF] Nazarov | Self-adjoint extensions for the Neumann Laplacian and applications[END_REF]. The asymptotic expansion for the particular problem we are considering in this paper is described in detail in [START_REF] Maz'ya | Asymptotic theory of elliptic boundary value problems in singularly perturbed domains[END_REF] and here, we just wish to remind the main results provided by the method of matched expansions at order one. To proceed, we need first to introduce two particular functions: the Green function G and the logarithmic capacity potential P . These two functions are defined in normalized geometries by the following equations

         -∆G = 0 in Ω \ {O} G = 0 on ∂Ω G(x) = 1 2π ln(1/|x|) + O |x|→0 (1)
         -∆P = 0 in Ξ := R 2 \ ω P = 0 on ∂Ξ P (ξ) = 1 2π ln(1/|ξ|) + O |ξ|→∞ (1).
(

Classical techniques of separation of variables (see e.g. [START_REF] Lebedev | Special functions and their applications[END_REF]) show that there exist constants G 0 , P 0 that depend only on the domains Ω,

ω such that G(x) -(2π) -1 ln |x| -1 -G 0 = O(|x|) for |x| → 0, and
P (ξ) -(2π) -1 ln |ξ| -1 -P 0 = O(|ξ| -1 ) for |ξ| → ∞.
The asymptotic analysis of Problem (1) also involves the gauge function (see [START_REF] Il'in | Study of the asymptotic behavior of the solution of an elliptic boundary value problem in a domain with a small hole[END_REF])

λ(δ) := 2π ln δ + 2π(P 0 -G 0 ) . ( 4
)
Finally, the global approximation of u δ is defined as an interpolation between a far field and a near field contribution as follows:

ûδ (x) := ψ(x/δ) v δ (x) + χ(x) V δ (x/δ) -χ(x)ψ(x/δ) m δ (x) (5)
where

v δ (x) := u 0 (x) + u 0 (0) λ(δ) G(x) V δ (ξ) := u 0 (0) λ(δ) P (ξ) m δ (x) := u 0 (0) λ(δ) 1 2π ln(δ/|x|) + P 0 . ( 6
)
In the expression above, the cut-off functions χ, ψ are two elements of

C ∞ (Ω) := {v| Ω | v ∈ C ∞ (R 2
)} such that χ(x) = 1 for |x| ≤ r 0 /2, χ(x) = 0 for |x| ≥ r 0 and ψ := 1 -χ. Here, r 0 > 0 is a given parameter such that D r 0 ⊂ Ω. The following well-known result provides an error estimate for u δ -ûδ H1 (Ω δ ) . For the proof, we refer the reader, for example, to Section 2.4.1 of [START_REF] Maz'ya | Asymptotic theory of elliptic boundary value problems in singularly perturbed domains[END_REF].

Proposition 2.1.

Considering u δ defined by [START_REF] Antoine | Wide frequency band numerical approaches for multiple scattering problems by disks[END_REF] and ûδ defined by [START_REF] Bonnaillie-Noël | Interactions between moderately close circular inclusions: the Dirichlet-Laplace equation in the plane[END_REF], there exist constants C, δ 0 > 0 independent of δ such that

u δ -ûδ H 1 (Ω δ ) ≤ C δ | ln δ| f L 2 (Ω) ∀δ ∈ (0, δ 0 ].
Note that, the constant C involved in the estimate above a priori depends on χ, ψ. Note also that, in the definition of χ, ψ, the parameter r 0 could be any positive number such that D r 0 ⊂ Ω. In particular, it can be chosen arbitrarily small. Looking at the explicit definition of ûδ given by ( 6), this implies the following result.

Proposition 2.2.

Consider u δ , v δ defined by ( 1), [START_REF] Bonnaillie-Noël | Interactions between moderately close circular inclusions: the Dirichlet-Laplace equation in the plane[END_REF]. For any disk D r ⊂ Ω with 0 < r ≤ r 0 , there exist constants C r , δ 0 > 0 independent of δ such that

u δ -v δ H 1 (Ω\Dr) ≤ C r δ | ln δ| f L 2 (Ω) ∀δ ∈ (0, δ 0 ]. (7)
This last result shows that v δ = u 0 +u 0 (0)λ(δ)G provides a reasonable approximation (for example δ | ln δ| ≈ 2.3 10 -9 for δ = 10 -10 ) of u δ at any fixed distance from the small hole. Thus, the far field of u δ appears as the superposition of the limit field u 0 and a field "radiated" by a point source located at the center of the hole. Numerically, u 0 , G can be approximated by functions u h 0 , G h using a standard finite element method (here, h refers to some mesh size) and define v h δ = u h 0 + u h 0 (0)λ(δ)G h . Then, [START_REF] Bonnaillie-Noël | On moderately close inclusions for the Laplace equation[END_REF] ensures that v h δ is a good approximation of the far field of u δ . This procedure is rather simple to implement and it has been proven in [START_REF] Bonnaillie-Noël | Interactions between moderately close circular inclusions: the Dirichlet-Laplace equation in the plane[END_REF] (see also [START_REF] Bonnaillie-Noël | Interactions between moderately close inclusions for the Laplace equation[END_REF][START_REF] Dambrine | Influence of a boundary perforation on the Dirichlet energy[END_REF][START_REF] Dambrine | A multiscale correction method for local singular perturbations of the boundary[END_REF][START_REF] Bonnaillie-Noël | On moderately close inclusions for the Laplace equation[END_REF][START_REF] Ben Hassen | Asymptotic formulas for the voltage potential in a composite medium containing close or touching disks of small diameter[END_REF][START_REF] Bonnetier | An elliptic regularity result for a composite medium with "touching" fibers of circular cross-section[END_REF] for slightly different problems 1 ) that it gives good results. However, it requires to solve two problems which we would like to avoid because it may be time consuming. Adapting this approach to the case of N inclusions would lead to N + 1 numerical solves. Similarly, looking for an approximation of u δ as sharp as the first M terms of its asymptotic expansion would lead to M numerical solves. From this perspective, for practical computations, a method involving only one numerical solve would be much more interesting. In the next section, we propose a model problem that can be discretized by means of any standard Galerkin method (with classical finite elements for example) with quasi-optimal approximation properties with respect to both h and δ. In addition, the numerical schemes obtained in this manner do not deteriorate as δ → 0.

Construction of a model problem

The model problem we wish to propose is formulated in [START_REF] Campbell | Une justification de la méthode de raccordement des développements asymptotiques appliquée à un problème de plaque en flexion. Estimation de la matrice d'impédance[END_REF]. The goal of the present section is to explain how we obtain this problem. To avoid having to compute both u 0 and G in the decomposition v δ = u 0 + u 0 (0)λ(δ)G, we will use the fact that the regular part of G belongs to H 1 (Ω). Let us decompose G under the form

G = s log + G with s log (x) := 1 2π χ(x) ln(1/|x|) and G ∈ H 1 0 (Ω) ∩ C 0 (Ω). ( 8
)
This allows us to write v δ as v δ = w δ + u 0 (0)λ(δ) s log with w δ := u 0 + u 0 (0)λ(δ) G. Let us express the coefficient u 0 (0)λ(δ) by means of w δ . According to the definition of u 0 and G it is clear that

w δ belongs to H 1 0 (Ω) ∩ C 0 (Ω)
, where C 0 (Ω) refers to the space of continuous functions on Ω. Moreover, observing that G = G 0 + Ĝ for some function Ĝ ∈ H 1 0 (Ω) ∩ C 0 (Ω) vanishing at 0, we find w δ (0) = u 0 (0)(1 + λ(δ) G 0 ) and so u 0 (0)λ(δ) = w δ (0)λ(δ)/(1 + λ(δ) G 0 ). Using Definition (4) of λ(δ), we deduce that

v δ = w δ + b δ (w δ ) s log with w δ := u 0 + u 0 (0)λ(δ) G and b δ (w δ ) := 2π w δ (0) ln δ + 2πP 0 . (9)
Let us emphasize that this expression for v δ is interesting because it involves only one unknown function which belongs to the variational space H 1 0 (Ω). Now, we need to derive a problem characterizing w δ . In the sense of distributions in Ω, there holds

-∆w δ = -∆(u 0 + u 0 (0)λ(δ) G) = f -b δ (w δ ) ∆ G.
Multiplying by w ∈ H 1 0 (Ω) and using Green's formula, we find that w δ verifies

a(w δ , w ) + b δ (w δ ) b log (w ) = Ω f w dx ( 10
)
where a(w δ , w

) := Ω ∇w δ • ∇w dx b log (w ) := Ω ∆ G w dx ∆ G(x) := (2π) -1 (∆χ)(x) ln |x| + 2∇χ(x) • ∇(ln |x|) . ( 11
)
Note that χ is equal to one in a neighbourhood of 0 so that ∆ G indeed belongs to C ∞ (Ω). As a remark, let us observe that for test functions w such that 0 / ∈ supp(w ), we have b log (w ) =

Ω ∇s log • ∇w dx.
Of course [START_REF] Bourlard | Coefficients of the singularities for elliptic boundary value problems on domains with conical points. III: Finite element methods on polygonal domains[END_REF] is not a valid variational formulation in H 1 (Ω) because the functional w δ → w δ (0) is not defined on this space. So it cannot be exploited directly for discretization and then numerical computation. This is the motivation for considering a regularized counterpart of [START_REF] Bourlard | Coefficients of the singularities for elliptic boundary value problems on domains with conical points. III: Finite element methods on polygonal domains[END_REF] where in b δ (w δ ), we replace w δ (0) by (2πδ) -1 ∂D δ w δ dσ, ∂D δ denoting the circle centered at 0 and of radius δ. Finally, this leads us to examine the following model problem,

Find wδ ∈ H 1 0 (Ω) such that a( wδ , w ) + bδ ( wδ ) b log (w ) = Ω f w dx ∀w ∈ H 1 0 (Ω), (12)
where a(•, •), b log (•) are defined in [START_REF] Brenner | The mathematical theory of finite element methods[END_REF] and where bδ ( wδ

) := 2π ln δ + 2πP 0 1 2πδ ∂D δ wδ dσ. ( 13
)
The variational formulation [START_REF] Campbell | Une justification de la méthode de raccordement des développements asymptotiques appliquée à un problème de plaque en flexion. Estimation de la matrice d'impédance[END_REF] perfectly makes sense for δ small enough and, in the next section, we show that it admits a unique solution so that wδ is well defined. Since Problem (12) differs from [START_REF] Bourlard | Coefficients of the singularities for elliptic boundary value problems on domains with conical points. III: Finite element methods on polygonal domains[END_REF], its solution wδ is a priori different from w δ . However we are going to show that wδ and w δ (defined by ( 9)) are close to each other, and that wδ + bδ ( wδ ) s log is a good approximation of the far field of u δ .

Remark 3.1. We could have proposed a formulation where, in [START_REF] Bourlard | Coefficients of the singularities for elliptic boundary value problems on domains with conical points. III: Finite element methods on polygonal domains[END_REF], the term w δ (0) is replaced by (πδ 2 ) -1 D δ w δ dx. The analysis we will develop and the results we will obtain would have been the same with this alternative choice. [START_REF] Campbell | Une justification de la méthode de raccordement des développements asymptotiques appliquée à un problème de plaque en flexion. Estimation de la matrice d'impédance[END_REF], a simple perturbation of a usual formulation allows to take into account the small hole. Therefore, with this approach, we can adapt classical codes at little cost. In this respect, this technique shares similarities with the extended finite element method (XFEM) [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF][START_REF] Dolbow | A finite element method for crack growth without remeshing[END_REF] and the generalized finite element method (GFEM) [START_REF] Duarte | An h-p adaptive method using clouds[END_REF][START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF]. Remark 3.3. In this paper, we have chosen to investigate the problem of the small hole with Dirichlet boundary condition in 2D only because in this case, the logarithmic term which appears in the asymptotic expansion of u δ makes the zero order approximation clearly unsatisfactory (see the discussion in the introduction). However, the present approach could allow to consider other problems of singular perturbation. It could also be adapted to obtain higher orders of approximation. In this case, new perturbation terms would have to be considered in the left hand side of [START_REF] Campbell | Une justification de la méthode de raccordement des développements asymptotiques appliquée à un problème de plaque en flexion. Estimation de la matrice d'impédance[END_REF].

Remark 3.2. It is worth noting that in

Analysis and discretization of the model problem

We first prove that ãδ (•, •) := a(•, •) + bδ (•) b log (•), the bilinear form appearing in the left hand side of ( 12), differs from a(•, •) by a small perturbation. This will allow to show that wδ is a relevant approximation of w δ .

Proposition 4.1.

There exists a constant C > 0 independent of δ such that

sup ϕ∈H 1 0 (Ω)\{0} | bδ (ϕ)| ϕ H 1 (Ω) ≤ C | ln δ| ∀δ ∈ (0, 1). ( 14
)
As a consequence, for δ small enough, ãδ (•, •) is coercive and Problem ( 12) has a unique solution wδ . Moreover, for any ε > 0, there exist constants C ε , δ 0 > 0 independent of δ such that

wδ -w δ H 1 (Ω) ≤ C ε δ 1-ε f L 2 (Ω) ∀δ ∈ (0, δ 0 ], (15)
where w δ is the function defined in [START_REF] Bonnetier | An elliptic regularity result for a composite medium with "touching" fibers of circular cross-section[END_REF].

Proof: First, we prove [START_REF] Chen | Numerical homogenization of well singularities in the flow transport through heterogeneous porous media[END_REF]. Consider the disk D r 0 introduced in the definition of χ that satisfies D r 0 ⊂ Ω. We have in particular supp(χ) ⊂ D r 0 . Take an arbitrary ζ ∈ C ∞ (Ω) such that supp(ζ) ⊂ D r 0 . Integration by parts and Cauchy-Buniakowski-Schwarz inequality show that

1 2πδ ∂D δ ζ dσ = 1 2π Dr 0 \D δ ∇(ln |x|) • ∇ζ dx ≤ | ln(r 0 /δ)| 2π ζ H 1 (Ω) . ( 16
)
As a consequence, for any ϕ ∈ C ∞ 0 (Ω), considering χϕ instead of ζ in (16) and using (13), we see that there exist constants C, C > 0 (whose values may change from one occurrence to another) independent of δ such that | bδ (ϕ

)| = | bδ (χϕ)| ≤ C| ln δ| -1/2 χϕ H 1 (Ω) ≤ C | ln δ| -1/2 ϕ H 1 (Ω) . Since C ∞ 0 (Ω) is dense into H 1 0 (Ω)
, this shows [START_REF] Chen | Numerical homogenization of well singularities in the flow transport through heterogeneous porous media[END_REF]. We deduce that for all ϕ ∈ H 1 0 (Ω), we have

|ã δ (ϕ, ϕ)| = |a(ϕ, ϕ) + bδ (ϕ) b log (ϕ)| ≥ C (1 -C | ln δ| -1/2 ) ϕ 2 H 1 (Ω) . ( 17
)
This guarantees that for δ small enough, Problem (12) has a unique solution wδ . To establish the second part of the statement, we use [START_REF] Ciarlet | Numerical solution of maxwell's equations in axisymmetric domains with the fourier singular complement method[END_REF] and write, for δ small enough,

w δ -wδ 2 H 1 (Ω) ≤ C |ã δ (w δ -wδ , w δ -wδ )| ≤ C |b δ (w δ ) -bδ (w δ )| |b log (w δ -wδ )|. ( 18
)
Let us focus on

|b δ (w δ ) -bδ (w δ )| = 2π ln δ + 2π P 0 w δ (0) - 1 2πδ ∂D δ w δ dσ . ( 19
)
From ( 9), we know that there holds

w δ = u 0 + u 0 (0)λ(δ) G with G = G 0 + Ĝ, Ĝ ∈ C ∞ (Ω).
For β ∈ R, we define the weighted norm

ϕ V 1 β (Ω) := ( |x| β ∇ϕ 2 L 2 (Ω) + |x| β-1 ϕ 2 L 2 (Ω) ) 1/2 , ( 20
)
and let V 1 β (Ω) refer to the completion of C ∞ (Ω \ {O}) := {v| Ω | v ∈ C ∞ (R 2
), v = 0 in a neighbourhood of 0 } with respect to this norm. We refer the reader to [START_REF] Kondratiev | Boundary-value problems for elliptic equations in domains with conical or angular points[END_REF] for more details on weighted Sobolev spaces. On the other hand, classical Kondratiev analysis (see [START_REF] Kozlov | Elliptic Boundary Value Problems in Domains with Point Singularities[END_REF]Chap.6]) allows to prove the decomposition u 0 = u 0 (0) + ũ0 , where ũ0 ∈ H 1 (Ω) ∩ V 1 -1+ε (Ω) for all ε > 0, with the estimate

|u 0 (0)| + ũ0 V 1 -1+ε (Ω) ≤ C ε f L 2 (Ω) . ( 21
)
This implies w δ = w δ (0) + (ũ 0 + u 0 (0)λ(δ) Ĝ). Conducting a calculus analogue to [START_REF] Ciarlet | The Fourier Singular Complement Method for the Poisson problem. Part II: axisymmetric domains[END_REF], replacing formally ζ(x) by |x| -1+ε ũ0 (x), we find that

| ∂D δ ũ0 dσ| ≤ C δ 2-ε f L 2 (Ω)
. Writing a Taylor expansion of Ĝ at x = 0 and using [START_REF] Duarte | An h-p adaptive method using clouds[END_REF], we obtain |u 0 (0)λ(δ

) ∂D δ Ĝ dσ| ≤ C δ 2 f L 2 (Ω) . We deduce w δ (0) - 1 2πδ ∂D δ w δ dσ ≤ C δ 1-ε f L 2 (Ω) . ( 22
)
Plugging this estimate in [START_REF] Dambrine | A multiscale correction method for local singular perturbations of the boundary[END_REF]

leads to |b δ (w δ ) -bδ (w δ )| ≤ C δ 1-ε f L 2 (Ω) .
Combining this inequality with [START_REF] Dambrine | Influence of a boundary perforation on the Dirichlet energy[END_REF], we obtain (15) as a direct consequence.

We have just proved that wδ is close to w δ . From the relation linking w δ to v δ , we deduce that wδ + bδ ( wδ ) s log is a good approximation of the far field of u δ .

Proposition 4.2.

For any disk D r ⊂ Ω with 0 < r ≤ r 0 and for any ε > 0, there exist constants C, δ 0 > 0 depending on r, ε but not on δ such that

u δ -( wδ + bδ ( wδ ) s log ) H 1 (Ω\Dr) ≤ C δ 1-ε f L 2 (Ω) ∀δ ∈ (0, δ 0 ]. (23)
Proof: Remembering that v δ = w δ + b δ (w δ ) (see ( 9)), where v δ , w δ are defined in ( 6), ( 9), and using the triangular inequality, we can write

u δ -( wδ + bδ ( wδ ) s log ) H 1 (Ω\Dr) ≤ u δ -v δ H 1 (Ω\Dr) + w δ -wδ H 1 (Ω\Dr) + |b δ (w δ ) -bδ ( wδ )| s log H 1 (Ω\Dr) . ( 24
)
In the previous proof, we have established that |b δ (w δ ) -bδ (w δ )| ≤ C δ 1-ε f L 2 (Ω) for some constant C > 0 independent of δ. Combining this with ( 14), we find

|b δ (w δ ) -bδ ( wδ )| ≤ |b δ (w δ ) -bδ (w δ )| + | bδ (w δ ) -bδ ( wδ )| ≤ C (δ 1-ε f L 2 (Ω) + w δ -wδ H 1 (Ω) ), (25)
for some constant C > 0 independent of δ. Plugging ( 25) in [START_REF] Il'in | Study of the asymptotic behavior of the solution of an elliptic boundary value problem in a domain with a small hole[END_REF] and using ( 7), ( 15), we finally obtain [START_REF] Hazard | A singular field method for maxwell's equations: Numerical aspects for 2d magnetostatics[END_REF].

Remark 4.1.

Working as in the previous proof, one can obtain a slightly more general result where the norm • H 1 (Ω\Dr) in the right hand side of ( 23) is replaced by the norm [START_REF] Dolbow | A finite element method for crack growth without remeshing[END_REF]) with β > 0. Remark 4.2. Making the additional assumption that the source term f verifies |x| -β f L 2 (Ω) < +∞ for some β > 0, and revisiting Estimate [START_REF] Efendiev | of Surveys and Tutorials in the Applied Mathematical Sciences[END_REF], we find that [START_REF] Ciarlet | The Fourier singular complement method for the Poisson problem. Part I: Prismatic domains[END_REF]

• V 1 β (Ω) (see
can be improved in wδ -w δ H 1 (Ω) ≤ C δ |x| -β f L 2 (Ω) , ∀δ ∈ (0, δ 0 ]. In this case, (23) becomes u δ -( wδ + bδ ( wδ ) s log ) H 1 (Ω\Dr) ≤ C δ | ln δ| |x| -β f L 2 (Ω) , ∀δ ∈ (0, δ 0 ].
To conclude, assume that we want to solve Formulation [START_REF] Campbell | Une justification de la méthode de raccordement des développements asymptotiques appliquée à un problème de plaque en flexion. Estimation de la matrice d'impédance[END_REF] by means of a Galerkin approach associated with a family of discrete subspaces (V h ) h>0 (in the numerical experiments, h will refer to the mesh size). We assume that there holds V h ⊂ H 1 0 (Ω) for all h > 0 . The natural discrete variational formulation associated with [START_REF] Campbell | Une justification de la méthode de raccordement des développements asymptotiques appliquée à un problème de plaque en flexion. Estimation de la matrice d'impédance[END_REF] writes

Find wh δ ∈ V h such that a( wh δ , ϕ h ) + bδ ( wh δ ) b log (ϕ h ) = Ω f ϕ h dx ∀ϕ h ∈ V h . ( 26
)
The coercivity of ãδ (•, •) = a(•, •) + bδ (•) b log (•) proven in Proposition 4.1 shows straightforwardly, by Cea's lemma, the result of quasi-optimal convergence wδwh

δ H 1 (Ω) ≤ C inf ϕ h ∈V h wδ -ϕ h H 1 (Ω) . ( 27
)
Combining this with Estimate [START_REF] Hazard | A singular field method for maxwell's equations: Numerical aspects for 2d magnetostatics[END_REF] proves that wh δ + bδ ( wh δ ) s log is a reasonable approximation of the far field expansion of u δ . The following proposition is one of the two main results (with Proposition 5.1 hereafter) of the present article. It establishes quasi-optimal convergence of the numerical method ( 26) both in δ and h.

Proposition 4.3.

Consider a finite dimensional space V h ⊂ H 1 0 (Ω). For any disk D r ⊂ Ω with 0 < r ≤ r 0 and for any ε > 0, there exists a constant C > 0 depending on r, ε but not on δ and h such that, for δ small enough,

u δ -( wh δ + bδ ( wh δ ) s log ) H 1 (Ω\Dr) ≤ C (δ 1-ε + | ln δ| -1 inf ϕ h ∈V h G -ϕ h H 1 (Ω) ) f L 2 (Ω) + C inf ϕ h ∈V h u 0 -ϕ h H 1 (Ω) . ( 28
)
Proof: The continuity estimate of bδ (see ( 14)) implies that there exists a constant C > 0 independent of δ such that u δ -( wh δ + bδ ( wh δ )s log ) H 1 (Ω\Dr) ≤ C u δ -( wδ + bδ ( wδ )s log ) H 1 (Ω\Dr) +C wδwh δ H 1 (Ω) . Proposition 4.2 already yields that u δ -( wδ + bδ ( wδ ) s log ) H 1 (Ω\Dr) ≤ C δ 1-ε f L 2 (Ω) for any ε > 0, so we only need to focus on the second term of the previous inequality. Since we have w δ = u 0 + u 0 (0)λ(δ) G, [START_REF] Kozlov | Elliptic Boundary Value Problems in Domains with Point Singularities[END_REF] allows us to write wδwh

δ H 1 (Ω) ≤ C w δ -wδ H 1 (Ω) + C inf ϕ h ∈V h w δ -ϕ h H 1 (Ω) ≤ C w δ -wδ H 1 (Ω) + C inf ϕ h ∈V h u 0 -ϕ h H 1 (Ω) + C |u 0 (0)λ(δ)| inf ϕ h ∈V h G -ϕ h H 1 (Ω) ≤ C (δ 1-ε + | ln δ| -1 inf ϕ h ∈V h G -ϕ h H 1 (Ω) ) f L 2 (Ω) + C inf ϕ h ∈V h u 0 -ϕ h H 1 (Ω) .
This finishes the proof.

To illustrate what kind of result the above proposition implies, assume for example that V h is the space of P 1 -Lagrange finite element functions constructed on a quasi-uniform regular triangulation of the domain Ω. In this situation, according to [START_REF] Lebedev | Special functions and their applications[END_REF], for any ε > 0 there exists a constant C > 0 independent of δ and h, such that u δ -( wh δ + bδ ( wh

δ ) s log ) H 1 (Ω\Dr) ≤ C (δ 1-ε + h) f L 2 (Ω)
. We also emphasize that the result of quasi-optimal convergence for (26) with a constant independent of δ discards any numerical locking effect. In other words, this assert the robustness of (26) as δ → 0.

Practical implementation of the perturbation

From the point of view of practical implementation, a natural idea consists in computing the perturbation term bδ (•) by means of the crude quadrature formula ∂D δ ϕ h dσ 2πδϕ h (0) for any ϕ h ∈ V h , which boils down to actually considering b δ (•) instead of bδ (•). In this section we examine the validity of such a substitution. We introduce the discrete formulation

Find w h δ ∈ V h such that a(w h δ , ϕ h ) + b δ (w h δ ) b log (ϕ h ) = Ω f ϕ h dx ∀ϕ h ∈ V h . ( 29
)
assuming that V h ⊂ C 0 (Ω) (this implies in particular that (29) has indeed a sense) is a Lagrange finite element space constructed on a quasi-uniform regular triangulation of the domain Ω. Let us prove that Problem (29) yields to a good approximation of the far field of u δ .

Proposition 5.1.

For any given h > 0, for δ > 0 small enough, Problem [START_REF] Maz | Asymptotic behavior of energy integrals under small perturbations of the boundary near corner and conic points[END_REF] has a unique solution w h δ . Moreover, if f ∈ H 2 (Ω) and if Ω is smooth, then for any disk D r ⊂ Ω with 0 < r ≤ r 0 and for any ε > 0, there exists a constant C > 0 depending on r, ε but not on δ and h such that, for δ small enough,

u δ -(w h δ + b δ (w h δ ) s log ) H 1 (Ω\Dr) ≤ C (δ| ln δ| + γ(δ, h) + | ln δ| -1 inf ϕ h ∈V h G -ϕ h H 1 (Ω) ) f H 2 (Ω) + C inf ϕ h ∈V h u 0 -ϕ h H 1 (Ω) . ( 30
)
In [START_REF] Maz'ya | Asymptotic expansions of eigenvalues of boundary value problems for the Laplace operator in domains with small openings[END_REF], the constant γ(δ, h) can be chosen such that γ(δ, h) = (δ+h 2 | ln h|)/(1-(1+| ln h|) 1/2 /| ln δ|).

Remark 5.1. Observe that for any given h > 0, there holds |γ(δ, h)| ≤ C (δ + h 2 | ln h|) for δ small enough. Actually, in the proof, we will see that the condition | ln h| 1/2 /| ln δ| = O( 1) is sufficient to guarantee well-posedness for Problem [START_REF] Maz | Asymptotic behavior of energy integrals under small perturbations of the boundary near corner and conic points[END_REF]. Note that this assumption is in accordance with the situation we want to consider, namely an obstacle small compare to the mesh size (δ << h).

Remark 5.2. The additional smoothness assumption on the source term is needed for technical reasons (see the proof of Lemma 5.1). The authors do not know if it can be weakened.

Proof: We first recall the discrete Sobolev inequality (see [11, Lemma 4.9.2])

ϕ h L ∞ (Ω) ≤ C (1 + | ln h|) 1/2 ) ϕ h H 1 (Ω) ∀ϕ h ∈ V h . ( 31
)
Here and in the sequel of this proof, C > 0 denotes a constant independent of δ, h which may change from one occurrence to another. Since b δ (ϕ h ) = 2π ϕ h (0)/(ln δ + 2πP 0 ), we deduce from (31) that, for δ small enough, for all ϕ h ∈ V h , we have

|a(ϕ h , ϕ h ) + b δ (ϕ h ) b log (ϕ h )| ≥ C α(δ, h) ϕ h 2 H 1 (Ω) , (32)
where α(δ, h)

:= 1 -β(δ, h) and β(δ, h) := (1 + | ln h|) 1/2 /| ln δ|.
It is clear that for a given h, β(δ, h) tends to zero as δ goes to zero. Therefore, Estimate [START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF] shows that a(•,

•) + b δ (•) b log (•) is
coercive for δ small enough. In this case, from the Lax-Milgram theorem, we infer that Problem (29) has a unique solution. Now, we wish to establish [START_REF] Maz'ya | Asymptotic expansions of eigenvalues of boundary value problems for the Laplace operator in domains with small openings[END_REF]. Thanks to the triangular inequality, we can write

u δ -(w h δ + b δ (w h δ ) s log ) H 1 (Ω\Dr) ≤ u δ -( wh δ + bδ ( wh δ ) s log ) H 1 (Ω\Dr) + w h δ -wh δ H 1 (Ω\Dr) + |b δ (w h δ ) -bδ ( wh δ )| s log H 1 (Ω\Dr) . ( 33
)
The first term of the right hand side of (33) has already been studied in Proposition 4.3. To handle the last term, we use [START_REF] Maz'ya | Asymptotic theory of elliptic boundary value problems in singularly perturbed domains[END_REF] to obtain

|b δ (w h δ ) -bδ ( wh δ )| ≤ |b δ (w h δ ) -b δ ( wh δ )| + |b δ ( wh δ ) -bδ ( wh δ )| ≤ C β(δ, h) w h δ -wh δ H 1 (Ω) + |b δ ( wh δ ) -bδ ( wh δ )|. ( 34
)
Let us estimate the quantity w h δ -wh δ H 1 (Ω) which appears both in [START_REF] Nazarov | Asymptotic conditions at a point, selfadjoint extensions of operators, and the method of matched asymptotic expansions[END_REF] and [START_REF] Nazarov | Asymptotic analysis of shape functionals[END_REF]. The coercivity inequality [START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF] and the definition of Problems ( 26), [START_REF] Maz | Asymptotic behavior of energy integrals under small perturbations of the boundary near corner and conic points[END_REF] provide

C α(δ, h) w h δ -wh δ 2 H 1 (Ω) ≤ |a(w h δ -wh δ , w h δ -wh δ ) + b δ (w h δ -wh δ ) b log (w h δ -wh δ )| ≤ |b δ ( wh δ ) -bδ ( wh δ )| |b log ( wh δ -wh δ )|.
Observing that b log is bounded on H 1 (Ω), we deduce that

w h δ -wh δ H 1 (Ω) ≤ C α(δ, h) -1 |b δ ( wh δ ) -bδ ( wh δ )|. ( 35
)
Plugging [START_REF] Nazarov | Self-adjoint extensions for the Neumann Laplacian and applications[END_REF] in [START_REF] Nazarov | Asymptotic conditions at a point, selfadjoint extensions of operators, and the method of matched asymptotic expansions[END_REF] and [START_REF] Nazarov | Asymptotic analysis of shape functionals[END_REF], we conclude that it is sufficient to control |b δ ( wh δ ) -bδ ( wh δ )| to prove [START_REF] Maz'ya | Asymptotic expansions of eigenvalues of boundary value problems for the Laplace operator in domains with small openings[END_REF]. We have

|b δ ( wh δ ) -bδ ( wh δ )| ≤ |b δ ( wδ -wh δ ) -bδ ( wδ -wh δ )| + |b δ ( wδ ) -bδ ( wδ )|. ( 36
)
Then, by definition of b δ , bδ , we find, for δ small enough,

|b δ ( wδ -wh δ ) -bδ ( wδ -wh δ )| = 2π ln δ + 2πP 0 1 2πδ ∂D δ ( wδ -wh δ )(0) -( wδ -wh δ ) dσ ≤ C wδ -wh δ L ∞ (Ω) (37)
Lemma 5.1 hereafter guarantees that if wδ

∈ W 2,∞ (Ω) := {v ∈ L ∞ (Ω) | ∂ α v ∈ L ∞ (Ω), |α| ≤ 2} 2 ,
then ( wh δ ) uniformly converges to wδ as h tends to zero, with the estimate

wδ -wh δ L ∞ (Ω) ≤ C h 2 | ln h| wδ W 2,∞ (Ω) . ( 38
)
To ensure such a regularity for wδ , let us assume that the source term f verifies f ∈ H2 (Ω). Using ( 12) and ( 13), we see that in the sense of distributions in Ω, there holds

-∆ wδ = f δ with f δ := f -bδ ( wδ )∆ G.
Thus, if f ∈ H 2 (Ω), then the theory of elliptic regularity asserts that wδ ∈ H 4 (Ω). Besides, Proposition 4.1 and Estimate [START_REF] Ciarlet | Numerical solution of maxwell's equations in axisymmetric domains with the fourier singular complement method[END_REF] imply

wδ H 1 (Ω) ≤ C f L 2 (Ω) . ( 39
)
We emphasize that in [START_REF] Scott | Optimal L ∞ estimates for the finite element method on irregular meshes[END_REF], the constant C > 0 is independent of δ. As a consequence, from Proposition 4.1 and (39), we get f δ H 2 (Ω) ≤ C f H 2 (Ω) . In this case, from the Sobolev imbedding theorem, we deduce that wδ ∈ C 2 (Ω) with

wδ W 2,∞ (Ω) ≤ C wδ H 4 (Ω) ≤ C f δ H 2 (Ω) ≤ C f H 2 (Ω) . ( 40
)
Collecting ( 37), ( 38) and (40), we find

|b δ ( wδ -wh δ ) -bδ ( wδ -wh δ )| ≤ C h 2 | ln h| f H 2 (Ω) . ( 41
)
On the other hand, concerning the second term of the right hand side of (36), writing the Taylor expansion of wδ at x = 0 and coming back to the definition of b δ , bδ yields

|b δ ( wδ ) -bδ ( wδ )| ≤ C δ f H 2 (Ω) . ( 42
)
Plugging (41) and (42) into (36) leads to

|b δ ( wh δ ) -bδ ( wh δ )| ≤ C (δ + h 2 | ln h|) f H 2 (Ω) . ( 43
)
Finally, combining (33), ( 34), ( 35), (43) and using Remark 4.2 allows to obtain [START_REF] Maz'ya | Asymptotic expansions of eigenvalues of boundary value problems for the Laplace operator in domains with small openings[END_REF].

In order to complete the previous analysis, now we state a result of uniform approximation of wδ by wh δ whose proof can be obtained working exactly as in [START_REF] Scott | Optimal L ∞ estimates for the finite element method on irregular meshes[END_REF].

On Figures 5,6, we display the behaviour of u

δ -w h δ1 -b δ (w h δ1 )s log H 1 (Ω\Dρ) , u δ -w h δ2 - b δ (w h δ2 )s log H 1 (Ω\Dρ) , u δ -w h δ3 -b δ (w h δ3 )s log H 1 (Ω\Dρ)
with respect to the mesh size in logarithmic scale. For the experiments of Figure 5, the cut-off function χ appearing in the definition of b log (•) (see [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]) is chosen equal to χ exp , an element of C ∞ (Ω) built with the exponential function. For the simulations of Figure 6, we take χ equal to χ pol ∈ C 3 (Ω) \ C 4 (Ω), a piecewise polynomial function of degree 7. We take g = 1 and δ = 10 -10 . We notice that with χ = χ exp , we obtain optimal rates of convergence. This is not the case for P3 approximation when we choose χ = χ pol ∈ C 3 (Ω) \ C 4 (Ω). However, we also remark that for the mesh sizes h considered here, the error is smaller when χ is a polynomial function (χ = χ pol ) than when χ is built with the exponential function (χ = χ exp ).

On Figure 7, we observe the behaviour of u δ -w h δ2 -b δ (w h δ2 )s log H 1 (Ω\Dρ) with respect to the mesh size in logarithmic scale and for different values of δ. Here, w h δ2 is the solution of ( 29) with V h = V h 2 . We take g = 1 + sin(θ). We see some thresholds in the convergence with respect to h: according to the value of δ, the error stops decreasing at some h 0 . This corresponds again to the error of the model. Estimates [START_REF] Bonnaillie-Noël | On moderately close inclusions for the Laplace equation[END_REF] indicates that this error behaves like δ | ln δ|. This is better than the error of the 0 order model (in | ln δ| -1 ), but when δ is not so small, it is natural that it appears. These thresholds are absent in the curves of Figure 1 because of the value of the source term. A natural approach to decrease the error consists in considering a model of higher order. Then, working as in Section 3, one can derive a model problem which does not suffer from numerical locking effect and whose solution yields a better approximation of u δ . We emphasize that at any order, this technique requires only one numerical resolution, and remains robust as δ → 0.

Figure 1 :

 1 Figure 1: Convergence w.r.t. the mesh sizeδ = 10 -1 , g = 1.

  approx. -H 1 norm -a = 0.00042099 0 order approx. -L 2 norm -a = 6.5877e-05 Model problem -H 1 norm -a = 1.0066 Model problem -L 2 norm -a = 2.0019

Figure 2 :

 2 Figure 2: Convergence w.r.t. the mesh sizeδ = 10 -2 , g = 1.

  approx. -H 1 norm -a = 0.00042099 0 order approx. -L 2 norm -a = 6.5877e-05 Model problem -H 1 norm -a = 1.0066 Model problem -L 2 norm -a = 2.0049

Figure 3 :

 3 Figure 3: Convergence w.r.t. the mesh sizeδ = 10 -4 , g = 1.

  approx. -H 1 norm -a = 0.00042099 0 order approx. -L 2 norm -a = 6.687e-05 Model problem -H 1 norm -a = 1.0066 Model problem -L 2 norm -a = 2.0069

Figure 4 :

 4 Figure 4: Convergence w.r.t. the mesh sizeδ = 10 -10 , g = 1.

Figure 5 :

 5 Figure 5: Convergence w.r.t. the mesh size for several orders of approximationδ = 10 -10 , g = 1.The cut-off function χ appearing in the definition of b log (•) (see[START_REF] Brenner | The mathematical theory of finite element methods[END_REF]) is in C ∞ (Ω).

  1 norm -a = 1.0237 P2 -H 1 norm -a = 2.1533 P3 -H 1 norm -a = 2.3395

Figure 6 :

 6 Figure 6: Convergence w.r.t. the mesh size for several orders of approximationδ = 10 -10 , g = 1.The cut-off function χ appearing in the definition of b log (•) (see[START_REF] Brenner | The mathematical theory of finite element methods[END_REF]) is chosen in C 3 (Ω) \ C 4 (Ω).

  1 norm -δ = 0.1 P2 -H 1 norm -δ = 0.031623 P2 -H 1 norm -δ = 0.01 P2 -H 1 norm -δ = 0.0031623 P2 -H 1 norm -δ = 0.001

Figure 7 :

 7 Figure 7: Convergence w.r.t. the mesh size for several values of δg = 1 + sin(θ).

This technique is also very close to singular complement methods (or singular function methods) which are used to compute efficiently the solution of elliptic partial differential equations in non smooth domains (see for example[10, 17, 15, 16,

23])

In this definition, we use the classical multi-index notation.

FreeFem++, http://www.freefem.org/ff++/.

Matlab, http://www.mathworks.se/.

Acknowledgments

The authors would like to thank Sergey A. Nazarov, of the Faculty of Mathematics and Mechanics of St. Petersburg State University, for useful discussions and remarks. Besides, the work of the first author was supported by the Academy of Finland (decision 140998).

Lemma 5.1. Assume that the solution of Problem [START_REF] Campbell | Une justification de la méthode de raccordement des développements asymptotiques appliquée à un problème de plaque en flexion. Estimation de la matrice d'impédance[END_REF] verifies wδ ∈ W 2,∞ (Ω). Then, for δ small enough, we have the estimate

where C > 0 is independent of δ, h.

Numerical experiments

Now, let us present the results of the numerical tests that we conducted in order to validate our theoretical conclusions. First, we detail the parameters used for the experiments. Let Ω (resp. ω δ ) be the disk centered at 0 of radius 1 (resp. δ). Remember that we denote Ω δ = Ω \ ω δ . We consider the problem of finding u δ ∈ H 1 0 (Ω δ ) such that

Admittedly ( 44) is not exactly of the same form as [START_REF] Antoine | Wide frequency band numerical approaches for multiple scattering problems by disks[END_REF]. However the analysis developed in the previous sections can be adapted in a straightforward manner to deal with (44) and the results are the same. For such a configuration, the exact solution u δ is given by

Note that, with ω δ = D δ , we have ω = ω 1 = D 1 so that the logarithmic capacity potential P defined by (3) verifies P (ξ) = (2π) -1 ln |ξ| -1 . As a consequence, the parameter P 0 appearing in the definition of b δ (•) (see ( 9)) satisfies P 0 = 0. For the computation of this parameter in other geometries, we refer the reader to [START_REF] Ransford | Computation of logarithmic capacity[END_REF]. Let us consider Ω h a polygonal approximation of the domain Ω. Introduce (T h ) h a shape regular family of triangulations of Ω h . Here, h denotes the average mesh size. Define the family of finite element spaces

) such that ϕ| τ ∈ P κ (τ ) for all τ ∈ T h , where P κ (τ ) is the space of polynomials of degree at most κ ∈ {1, 2, 3} on the triangle τ . We will denote w h δ1 , w h δ2 and w h δ3 the numerical solutions of (29) obtained respectively with V h 1 , V h 2 and V h 3 . The cut-off function χ appearing in the definition of b log (•) (see [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]) is chosen in C ∞ (Ω) (except for the simulation of Figure 6) with χ(|x|) = 1 for |x| ≤ 0.25 and χ(|x|) = 0 for |x| ≥ 0.5. The errors are expressed in the norms • L 2 (Ω\Dρ) and • H 1 (Ω\Dρ) with ρ = 0.15. For the computations, we use the FreeFem++ 3 software while we display the results with Matlab 4 .

On Figures 1, 2, 3 and 4, we represent the behaviour of u

with respect to the mesh size in logarithmic scale. Figures 1,2, 3 and 4 correspond respectively to δ = 10 -1 , δ = 10 -2 , δ = 10 -4 and δ = 10 -10 . Here, u h 0 is the standard P1 approximation of u 0 , the 0 order approximation of u δ defined by [START_REF] Babuška | On locking and robustness in the finite element method[END_REF]. Moreover, w h δ1 is the solution of (29) with V h = V h 1 (again P1 approximation). We take g = 1. As predicted at the end of Section 2, we observe that the approximation of u δ by u h 0 does not provide satisfactory results (even for δ = 10 -10 ). This is due to the error in the model, of order | ln δ| -1 , which decays very slowly as δ tends to zero. Conversely, w h δ1 + b δ (w h δ1 )s log appears as a good approximation of u δ and the rates of convergence are as expected. Moreover, the curves for δ = 10 -10 confirm the absence of any locking phenomenon for this numerical scheme.


Eduardo F Costa
email: efcosta@icmc.usp.br.

Benoîte De Saporta
email: benoite.de-saporta@umontpellier.fr

Linear minimum mean square filters for Markov jump linear systems

New linear minimum mean square estimators are introduced in this paper by considering a cluster information structure in the filter design. The set of filters constructed in this way can be ordered in a lattice according to the refines of clusters of the Markov chain, including the linear Markovian estimator at one end (with only one cluster) and the Kalman filter at the other hand (with as many clusters as Markov states). The higher is the number of clusters, the heavier are pre-compuations and smaller is the estimation error, so that the cluster cardinality allows for a trade-off between performance and computational burden. In this paper we propose the estimator, give the formulas for pre-computation of gains, present some properties, and give an illustrative numerical example.

Introduction

There is a vast number of applications benefiting from the nice properties of the Kalman filter (KF). Among these properties, the possibility of pre-computation of gains [1,7] is of much relevance for applications. However, in some cases pre-computation is not possible or viable due to missing a priori relevant information. This is the case when using the KF to estimate the state of Markov jump linear systems (MJLS), since the parameters are not known prior to the current time instant k, in fact they depend on the Markov chain current state θ(k). Then, to use KF for MJLS one needs to do either online computation of the gains or offline pre-computation of a number of sample path dependent gains, a figure that grows exponentially with time.

This drawback of KF for MJLS is one of the main motivations behind the emergence of other filters for this class of systems, see e.g. [2,4,5,6]. Among them, the one that is closer to the KF in terms of structure and performance is the linear minimum mean square estimator (LMMSE) first introduced for MJLS in [3]. Indeed, as we shall see later, the LMMSE computation relies on coupled Riccati equations that are quite similar to the ones arising in Kalman filtering. Another similarity is that both are optimal in the mean square error sense, although under different constraints. The main dissimilarity with KF lies in the fact that, instead of having path dependent gains, the LMMSE has sets of N gains, N being the cardinality of the Markov state space (which we assume finite in this paper), that are precomputed based on the system matrices and the initial distribution of θ. During application, after observing θ(k) one picks the corresponding gain from the precomputed sets of gains. In this way, one obtains the best estimate for the system state x k among all estimators that are linear and Markovian, these being precisely the constraints we mentioned before.

However, there is currently no intermediary solution between the KF and the LMMSE in literature. In this paper we provide a "lattice" of filters bridging the KF to the LMMSE by relaxing the markovianity constraint and allowing clustered information of the Markov chain to be considered when designing the filter gains. By clustered information, we mean that we have a partition of the state space of the Markov chain into several classes called clusters, and we observe the trajectory of classes the chain belongs to along time. The clusters may be chosen as they are considered as a design parameter, establishing a trade-off between complexity and performance that can be explored in practical systems aiming at the best feasible performance. At one extreme when only one cluster is taken into account our filter is equivalent to the LMMSE, and at another extreme with N clusters we retrieve the Kalman filter; intermediary number of clusters leads to filters with variable performances and computational burden. Reasonably enough, the higher is the number of clusters the smaller is the attained estimation error and higher is the number of gains to compute.

We start with a simple, precise formulation of the optimal estimation problem in Section 2, with the estimator in the classical form of Luenberger observers. We then proceed in Section 3 to a constructive proof that evaluates the estimation error and uses the completion of squares method to obtain the optimal gains. Some remarks on how the proposed class of estimators includes both the KF and LMMSE, and on the number of gains and Riccati-like equations to be precomputed, are presented. Some variants of the studied optimization problem and how to extend optimality to general estimators are briefly discussed in Section 4. We have also included a numerical example in Section 5 comparing the computational burden in terms of CPU time and the estimation error computed both via Monte Carlo simulation and via the proposed formula. The example makes clear that the performance is strongly dependent on the number of clusters and how the Markov states are distributed in the clusters.

Problem formulation

Consider the MJLS

x k+1 = A θ(k) x k + G θ(k) w k y k = L θ(k) x k + H θ(k) w k , k ≥ 0, (1)
with initial condition x 0 ∼ N (x, Ψ) where N (x, Ψ) is the normal distribution with mean x and covariance matrix Ψ. The variable θ(k) denotes the state of a Markov chain with finite state space {1, 2, . . . , N } and initial distribution π 0 = [Pr(θ(0) = 1)

• • • Pr(θ(0) = N )].
The noise sequence w is independent from x 0 and the Markov chain θ, E[w k ] = 0 and E[w k w k ] is the identity matrix for all k. We assume that G i H i = 0 and H i H i > 0, 1 ≤ i ≤ N . We shall consider Luenberger observers in the form

1 x k+1 = A θ(k) x k + M k (y k -L θ(k) x k ), (2)
where matrix M k is referred to as the filter gain, and the initial estimate is given by x 0 = x. This produces an estimation error x = xx satisfying

x k+1 = (A θ(k) -M k L θ(k) ) x k + (G θ(k) -M k H θ(k) )w k , (3)
and

x 0 = x 0 -x ∼ N (0, Ψ).
As for the Markov chain, we consider a partition S 1 , . . . , S N C for its state space, and employ the variable ρ(k) to indicate the partition being visited at time k, that is,

ρ(k) = N C m=1 m × 1 {θ(k)∈Sm} .
We assume that the observations of ρ up to time k are available to calculate the filter gains. We also assume that the jump variable θ at time k is available, however in the clustered information filter we do not take into account its past values, that is, θ(0), . . . , θ(k -1) are not taken into account when calculating the gain (to avoid an excessive number of branches, as explained earlier). Moreover, the gain should not depend on future information, as it has to be implemented at every time instant k. Therefore, we impose that the filter matrices at time instant k are in the form

M k = h k (ρ(0), . . . , ρ(k -1), θ(k)), k ≥ 0 (4)
for measurable functions h k . Gains satisfying this constraint are referred to as feasible gains. We are interested in obtaining the minimum mean square state estimation at a given time s ≥ 0, leading to the optimization problem min M0,...,Ms

E{ x s -x s 2 |R s }, s.t. (4), (5)
where we write R s = {ρ(0), . . . , ρ(s), y(0), . . . , y(s), θ(s)} for ease of notation. We refer to the filter (3) satisfying ( 5) as the clustered information LMMSE, or CLMMSE for short.

Clustered information LMMSE computation

Consider the MJLS in ( 1) and the filter in ( 3) with an arbitrary sequence of feasible gains

M = {M k , k ≥ 0}. For each k ≥ 0, 0 ≤ i ≤ N and 0 ≤ m ≤ M , 0 ≤ m ≤ k -1, we define X 0,..., k-1 ,i,k (M ) = E( x k x k 1 {ρ(0)= 0 ,...,ρ(k-1)= k-1 ,θ(k)=i} ). (6)
The variable X plays an important role in the derivation of the formula for the optimal filter because the optimal gains

M = {M k , 0 ≤ k ≤ s} are such that X 0,..., k-1 ,i,k (M ) = Y 0,..., k-1 ,i,k
where Y is the solution of a Riccati-like equation, as we shall see in the next theorem. The physical interpretation for X is that, when divided by Pr(ρ(0 Consider trivial gains M k = 0, so that x k = x k and E( x k+1 x k+1 ) = E( x k x k ) + 1. Moreover, the estimation error and the Markov state are independent because all modes are identical, then

) = 0 , . . . , ρ(k -1) = k-1 , θ(k) = i) it
A i = G i = 1, 1 ≤ i ≤ 3, Ψ = 1, x = 0, the clusters S 1 = {1, 2}, S 2 = {3}, initial distribution π 0 = [0.5 0.
X 0 ,..., k-1 ,i,k (M ) = E( x k x k )Pr(ρ(0) = 0 , . . . , ρ(k -1) = k-1 , θ(k) = i). For instance, for 0 = 1, i = 2 and k = 1 we have X 1,2,1 = E( x 1 x 1 )Pr(ρ(0) = 1, θ(1) = 2) = 2 Pr(θ(1) = 2|ρ(0) = 1)Pr(ρ(0) = 1) = 2 Pr(ρ(0) = 1) j∈S 0 Pr(θ(1) = 2|ρ(0) = 1, θ(0) = j)Pr(θ(0) = j|ρ(0) = 1) = 2 × 0.8 × (0.4 × (0.5/0.8) + 0 × (0.3/0.8)) = 0.4. Similarly, X 1,1,1 = 1.1, X 1,3,1 = 0.1, X 2,1,1 = 0.2, X 2,2,1 = 0 and X 2,3,1 = 0.2.
Note from (6) that summing X in the indexes corresponding to and i we obtain E( x k x k ), e.g. from the above we have E( x 1 x 1 ) = 2.

Coupled Riccatis of the CLMMSE

the optimal gain sequence can be (pre-)computed based on the following sets of matrices. Let

Y i,0 = π i (0)Ψ for each 1 ≤ i ≤ N . For each k ≥ 1, let 0 ≤ m ≤ k -1 and compute for each 0 ≤ i ≤ N and 0 ≤ m ≤ M , Y 0,..., k-1 ,i,k =                    0, p 0,..., k-1 ,i,k = 0, j∈ S p ji A j Y 0 ,..., k-2 ,j,k-1 A j + p 0,..., k-2 ,j G j G j + A j Y 0,..., k-2 ,j,k-1 L j L j Y 0,..., k-2 ,j,k-1 L j + p 0,..., k-2 ,j,k-1 H j H j -1 L j Y 0,..., k-2 ,j,k-1 A j , otherwise, (7)
where we denote

p 0,..., k-1 ,i,k = Pr(ρ(0) = 0 , . . . , ρ(k -1) = k-1 , θ(k) = i), S = {j ∈ S k-1 : p 0,..., k-2 ,j,k-1 = 0}.
Theorem 3.2 Given the realization of θ(k), k ≤ s, and the corresponding cluster observations ρ(0), . . . , ρ(s -1), the gains M = {M 0 , . . . , M s } of the LMMSE can be computed for each k ≤ s as

M k =      0, Pr(ρ(0), . . . , ρ(k -1), θ(k)) = 0, A θ(k) Y ρ(0),...,ρ(k-1),θ(k),k L θ(k) L θ(k) Y ρ(0),...,ρ(k-1),θ(k),k L θ(k) + Pr(ρ(0), . . . , ρ(k -1), θ(k))H θ(k) H θ(k) -1 , otherwise, (8)
where Y is given in (7). Moreover, the conditional second moment of the estimation error is given, for each k ≤ s, by

X ρ(0),...,ρ(k-1),θ(k),k (M ) = Y ρ(0),...,ρ(k-1),θ(k),k , (9)
and it is optimal in the sense that, for any gain sequence M = {M 0 , . . . , M s },

X ρ(0),...,ρ(k-1),θ(k),k (M ) ≤ X ρ(0),...,ρ(k-1),θ(k),k (M ), 0 ≤ k ≤ s. ( 10
)
Proof We start showing that (9) and ( 10) are true for the gains prescribed in [START_REF] Zhang | Mode-dependent h∞ filtering for discrete-time Markovian jump linear systems with partly unknown transition probabilities[END_REF]. We proceed by induction in k. For the time instant k = 0 we have that the initial estimate is given by x 0 = x, yielding x 0 = x 0x ∼ N (0, Ψ) (irrespectively of the filter gains), hence

X i,0 (M ) = X i,0 (M ) = E( x 0 x 0 1 {θ(0)=i} ) = E( x 0 x 0 )E(1 {θ(0)=i} ) = Ψπ i (0) = Y i,0 . (11)
By the induction hypothesis we assume that (9) and ( 10) are valid for some 0 ≤ k < s. In order to complete the induction, we now consider the time instant k + 1. For a filter with an arbitrary sequence of easible gains, denoted by M = {M 0 , . . . , M s }, and a given realization ρ(0) = 0 , . . . , ρ(k) = k , θ(k + 1) = i we have X 0,..., k ,i,k+1 (M ) = E( x k+1 x k+1 1 {ρ(0)= 0,...,ρ(k)= k ,θ(k+1)=i} ).

Note that the above quantity turns out to be zero whenever p 0,..., k ,i,k+1 = 0, irrespectively of M , which makes (9) and (10) trivially true for k + 1 in this case. Now, in case p 0,..., k ,i,k+1 = 0 we write

X 0,..., k ,i,k+1 (M ) = E( x k+1 x k+1 1 {ρ(0)= 0 ,...,ρ(k)= k ,θ(k+1)=i} ) = E (A θ(k) -M k L θ(k) ) x k + (G θ(k) -M k H θ(k) )w k × (A θ(k) -M k L θ(k) ) x k + (G θ(k) -M k H θ(k) )w k 1 {ρ(0)= 0 ,...,ρ(k)= k ,θ(k+1)=i} = E [(A θ(k) -M k L θ(k) ) x k x k (A θ(k) -M k L θ(k) ) + (G θ(k) -M k H θ(k) )w k w k (G θ(k) -M k H θ(k) ) ]1 {ρ(0)= 0,...,ρ(k)= k ,θ(k+1)=i} = E N j=1 [(A θ(k) -M k L θ(k) ) x k x k (A θ(k) -M k L θ(k) ) + (G θ(k) -M k H θ(k) )w k w k (G θ(k) -M k H θ(k) ) ]1 {ρ(0)= 0,...,ρ(k)= k ,θ(k+1)=i,θ(k)=j} = E j∈S k [(A j -M k L j ) x k x k (A j -M k L j ) + (G j -M k H j )w k w k (G j -M k H j ) ]1 {ρ(0)= 0,...,ρ(k-1)= k-1 ,θ(k+1)=i,θ(k)=j} (12)
where the last equality comes from the fact that Pr(ρ(k) = k , θ(k) = j) = 0 whenever j is not in the cluster S k , and Pr(ρ(k) = k , θ(k) = j)) = Pr(θ(k) = j) otherwise. Resuming the above calculation and writing R = {ρ(0) = 0 , . . . , ρ(k -1) = k-1 , θ(k + 1) = i, θ(k) = j} for ease of notation, we have:

X 0 ,..., k ,i,k+1 (M ) = j∈S k Pr(R) (A j -M k L j )E x k x k R)(A j -M k L j ) + (G j -M k H j )E w k w k R)(G j -M k H j ) . ( 13
)
From basic properties of the Markov chain we have that ρ( ) and θ(k + 1) are conditionally independent given θ(k), for any 0 ≤ ≤ k -1. Moreover, from (1), ( 3) and ( 4) it can be shown that x k and θ(k + 1) are conditionally independent given θ(k), hence we may eliminate θ(k + 1) = i from the first conditional expectation in (13), leading to

X 0,..., k ,i,k+1 (M ) = j∈S k Pr(θ(k + 1) = i|ρ(0) = 0 , . . . , ρ(k -1) = k-1 , θ(k) = j)p 0,..., k-1 ,j,k × (A j -M k L j )E x k x k ρ(0) = 0 , . . . , ρ(k -1) = k-1 , θ(k) = j)(A j -M k L j ) + (G j -M k H j )I(G j -M k H j ) = j∈ S p ji (A j -M k L j )X 0,..., k-1 ,j,k (M )(A j -M k L j ) + p 0,..., k-1 ,j,k (G j -M k H j )(G j -M k H j ) , (14)
where we denote S = {j ∈ S k : p 0,..., k-1 ,j,k = 0}. We now turn our attention to the optimality of M . Consider a feasible gain sequence in the form

M = {M 0 , . . . , M k-1 , M k },
where M k is the variable to be determined; since X 0,..., k-1 ,j,k ( M ) is a function of M 0 , . . . , M k-1 only, we can use the induction hypothesis to write X 0,..., k-1 ,j,k ( M ) = Y 0,..., k-1 ,j,k , (15) X 0 ,..., k-1 ,j,k ( M ) ≤ X 0,..., k ,j,k (M ).

(16)

Eq. ( 16) allows to write (A j -M k L j )(X 0,..., k-1 ,j,k ( M ) -X 0,..., k-1 ,j,k (M ))(A j -M k L j ) ≤ 0, irrespectively of M k , and using (14) we evaluate

X 0,..., k ,i,k+1 ( M ) -X 0,..., k ,i,k+1 (M ) = j∈ S p ji (A j -M k L j )(X 0,..., k-1 ,j,k ( M ) -X 0,..., k-1 ,j,k (M ))(A j -M k L j ) ≤ 0 (17)
Also, by plugging (15) into (14),

X 0 ,..., k ,i,k+1 ( M ) = j∈ S p ji (A j -M k L j )Y 0 ,..., k-1 ,j,k ( M )(A j -M k L j ) + p 0 ,..., k-1 ,j,k (G j -M k H j )I(G j -M k H j ) , (18)
and, by completing squares and denoting Φ = L j Y 0,..., k-1 ,j,k L j + Pr(ρ(0) = 0 , . . . , ρ(k -1) = k-1 , θ(k) = j)H j H j for brevity, we obtain

X 0,..., k ,i,k+1 ( M ) = j∈ S p ji A j Y 0,..., k-1 ,j,k A j + p 0 ,..., k-1 ,j,k G j G j + M k -A j Y 0,..., k-1 ,j,k L j Φ -1 Φ M k -A j Y 0,..., k-1 ,i,k L j Φ -1 -A j Y 0,..., k-1 ,j,k L j Φ -1 L j Y 0,..., k-1 ,j,k A j , (19)
thus making clear that the minimal X is attained by setting

M k = M k = g( 0 , . . . , k-1 , j) = A j Y 0 ,..., k-1 ,i,k L j Φ -1
= A j Y 0,..., k-1 ,i,k L j L j Y 0,..., k-1 ,j,k L i + p 0,..., k-1 ,j,k H j H j whenever p 0,..., k-1 ,j,k = 0, confirming the second equation in [START_REF] Zhang | Mode-dependent h∞ filtering for discrete-time Markovian jump linear systems with partly unknown transition probabilities[END_REF]; the inverse always exists because we have assumed H i H i > 0. If j is such that p 0,..., k-1 ,j,k = 0 then the gain M k is immaterial for the error covariance, indeed we see from (18) that such gain is not accounted for, so that one can pick M k = 0, confirming the first equation in [START_REF] Zhang | Mode-dependent h∞ filtering for discrete-time Markovian jump linear systems with partly unknown transition probabilities[END_REF]. Chosing the gain as above we get the gain sequence M = {M 0 , . . . , M k } and X 0,..., k ,j,k+1 (M ) ≤ X 0 ,..., k ,j,k+1 ( M ), so that (16) produces X 0,..., k ,j,k+1 (M ) ≤ X 0 ,..., k ,j,k+1 (M ), which confirms (10) for the time instant k + 1. Substituting M k = M k in (19) we get after some manipulations

X 0,..., k ,i,k+1 (M ) = j∈ S p ji A j Y 0,..., k-1 ,j,k A j + p 0,..., k-1 ,j,k G j G j + A j Y 0,..., k-1 ,j,k L j L j Y 0,..., k-1 ,i,k L j + p 0 ,..., k-1 ,j,k H j H j L j Y 0,..., k-1 ,j,k A j = Y 0 ,..., k ,i,k+1 ,
which confirms (9) for k + 1, thus completing the induction. It remains only to show the optimality of M in terms of ( 5). This follows directly from (10), in fact,

E{( x k ) x k } = E(( x k ) x k 1 {ρ(0)= 0 ,...,ρ(k-1)= k-1 ,θ(k)=i} ) = tr(X 0 ,..., k-1 ,i,k (M )) ≤ tr(X 0 ,..., k-1 ,i,k (M )) = E( x k x k 1 {ρ(0)= 0,...,ρ(k-1)= k-1 ,θ(k)=i} ) = E{ x k x k },
where all sums are in the indexes 0 ≤ i ≤ N, 0 ≤ m ≤ M, 0 ≤ m ≤ k -1 and we denote by x k the estimation error associated with the gain M .

Properties of the CLMMSE

Number of matrices to be computed and stored

For each 0 ≤ k ≤ s -1, we compute N N k C matrices on the left hand side of ( 7), hence we have (up to) this number of recursive Riccati equations to solve. We also have the computation and storage of an equal number of gains. Then, to obtain the state estimate at time s, we have to store a total of N (N s C -1)(N C -1) -1 gains when N C = 1, and sN gains otherwise. Regarding the number of matrix inverses, one may invert each Y given by (7) and store it at time step k for the forthcoming iterates, hence we have a total of (up to) N (M

(s-1) C -1)(M C -1) -1 inverses.

Filtering in the entire interval 0 ≤ k ≤ s

Note from (8) that, given a realization of the Markov chain θ(k), k ≥ 0, the time instant s involved in the problem formulation (5) affects only the cardinality of the optimal gain sequence M . More precisely, if {M k , 0 ≤ k ≤ s} is the gain sequence attaining (5), and, if we replace s with ≤ s in (5) and obtain the new optimal gain sequence { M k , 0 ≤ k ≤ } (considering the same Markov chain realization), then we have that M k = M k , 0 ≤ k ≤ min( , s). This is consistent with the sense of optimality in (10), and is in perfect harmony with the theory of both Kalman filter and the standard LMMSE. As a consequence, the provided clustered information LMMSE is also a solution for the multiobjective problem min M0,...,Ms

{E{ x 0 -x 0 2 |R 0 }, . . . , E{ x s -x s 2 |R s },
or for any linear combination of mean square errors writen in the form

min M0,...,Ms   0≤k≤s α k E{ x k -x k 2 |R k }   .

Linking the Kalman filter and the standard LMMSE

It is simple to see that we retrieve the standard LMMSE when we consider only one partition S 1 = {1, . . . , N }. In fact, in this setup we have p 0,..., k-1 ,i,k = P (θ(k) = i) and one can check by inspection that (7) and the LMMSE Riccati equation [3, Eq. XX] are identical. As for the Kalman filter, if we set

S i = {i}, 1 ≤ i ≤ N , then p 0 ,..., k-1 ,i,k = Pr(ρ(0) = 0 , . . . , ρ(k -1) = k-1 , θ(k) = i) = Pr(θ(0) = 0 , . . . , θ(k -1) = k-1 , θ(k) = i) = π 0 (0)p 0, 1 • • • p k-1 ,i , S = {j ∈ S k-1 : p 0,..., k-2 ,j,k-1 = 0} = { k-1 }, if p 0,..., k-2 , k-1 ,k-1 = 0 and (7) is reduced to Y 0,..., k-1 ,i,k = A j Y 0,..., k-2 ,j,k-1 A j + p 0,..., k-2 ,j G j G j + A j Y 0,..., k-2 ,j,k-1 L j L j Y 0,..., k-2 ,j,k-1 L j + p 0,..., k-2 ,j,k-1 H j H j -1 L j Y 0,..., k-2 ,j,k-1 A j .
From Theorem 3.2 we have

Y 0,..., k-1 ,i,k = X 0,..., k-1 ,i,k (M ) =E( x k x k |θ(0) = 0 , . . . , θ(k -1) = k-1 , θ(k) = i)π 0 (0)p 0, 1 • • • p k-1 ,i so that, writing Z 0,..., k-1 ,i = E( x k x k |θ(0) = 0 , . . . , θ(k -1) = k-1 , θ(k) = i)
, substituting in the above equation for Y and manipulating (cancelling the ps and πs) yields

Z 0,..., k-1 ,i,k = A j Z 0,..., k-2 ,j,k-1 A j + G j G j + A j Z 0,..., k-2 ,j,k-1 L j L j Z 0,..., k-2 ,j,k-1 L j + H j H j -1 L j Z 0,..., k-2 ,j,k-1 A j , (20)
which is the usual Riccati difference equation appearing in Kalman filters. This means that Y 0,..., k-1 ,i,k is equal to the Kalman covariance matrix multiplied by the probability that the Markov chain visits 0 , . . . , k-1 , i. In the gain formula [START_REF] Zhang | Mode-dependent h∞ filtering for discrete-time Markovian jump linear systems with partly unknown transition probabilities[END_REF], this probability is cancelled, yielding that the Kalman gain coincide with M k . Concluding, we have the Kalman filter and the markovian LMMSE in opposite "extremes" of the CLMMSE, and a lattice of estimators between them, depending on how the Markov states are arranged in clusters.

General LMMSE

Consider linear estimators of the general form

z k+1 = F k z k + Ḡk y k , (21)
where matrices F k and Bk , k ≥ 0, are the optimization variables replacing M k in the problem (5); consider also that F k = f k (ρ(0), . . . , ρ(k -1), θ(k)) and Ḡk = g k (ρ(0), . . . , ρ(k -1), θ(k)), where f k , g k are measurable functions. It can be demonstrated that the optimal estimate satisfies

z k = x k , 0 ≤ k ≤ s, a.s.,
which is produced by setting Ḡk = M k and

F k = A θ(k) -M k L θ(k)
where M k is as in [START_REF] Zhang | Mode-dependent h∞ filtering for discrete-time Markovian jump linear systems with partly unknown transition probabilities[END_REF], thus retrieving the Luenberger observer form (3) and the solution given in Theorem 3.2. This is not surprising since the "innovation form"

F k = A θ(k) -M k L θ(k)
for some M k is necessary for some basic properties of an observer to be fulfilled, e.g. for x k = z kx k to remain zero a.s. for k ≥ 1 in cases when z 0 = x 0 a.s. and there is no additive noise in the state (G i = 0, i = 1, 2, . . . , N ).

Illustrative example

We have applied the CLMMSE to the system given in [?]. The system data is reproduced below for ease of reference.

A 1 = 0 -0.
Every possible cluster configuration has been taken into account, and the aim was to estimate the state at time instant s = 10. The mean square error was calculated by means of ( 6) and ( 9), which lead to E( x10 2 ) = 0,..., 9,i trace(Y 0,..., 9,i,10 ); results have been confirmed by Monte Carlo simulation. Figure 1 shows the obtained results in groups according to the number of clusters N C . As expected, the standard LMMSE (with N C = 1) presented the largest estimation error and the Kalman filter (N C = 4) features the smallest one. The performance of other filters with "intermediary configurations" (N C = 2, 3) is similar to the LMMSE regarding the error, hence they are not much appealing in view of their higher complexity. Note that in this particular example, the modes are similar to each other (only two parameters of A change).

N C . As expected, the standard LMMSE (with N C = 1) presented the largest estimation error and the Kalman filter (N C = 4) features the smallest one. The performance of other filters with "intermediary configurations" (N C = 2, 3) is similar to the LMMSE regarding the error, hence they are not much appealing in view of their higher complexity. Note that in this particular example the modes are similar to each other (only two parameters of A change). Now, let us

respectively. The results are displayed in Figure 2. We now can clearly distinguish two groups of filters, one with average errors around 10 6 and a second one around 10 4 . In this setup, there is a tendency for better performance when θ = 4 is isolated from other states, e.g. with S 1 = {1, 2, 3} and S 2 = {4} we have E(∥ x 10 ∥ 2 ) ≈ 7, 566.69, while S 1 = {1, 2} and S 2 = {3, 4} lead to E(∥ x 10 ∥ 2 ) ≈ 1, 286, 306.18. Moreover, considering that the (hard to implement) Kalman filter yields E(∥ x 10 ∥ 2 ) ≈ 5, 152.95 we see that the filter with configuration S 1 = {1, 2, 3} and S 2 = {4} is quite competitive in this scenario.

respectively. The results are displayed in Figure 2. We now can clearly distinguish two groups of filters, one with average errors around 106 and a second one around 104. In this setup, there is a tendency for better performance when θ = 4 is isolated from other states, e.g. with S 1 = {1, 2, 3} and S 2 = {4} we have E( x10 2 ) 7, 566.69, while S 1 = {1, 2} and S 2 = {3, 4} lead to E( x102 ) 1, 286, 306.18. Moreover, considering that the (hard to implement) Kalman filter yields E( x102 ) 5, 152.95 we see that the filter with configuration S 1 = {1, 2, 3} and S 2 = {4} is quite competitive in this scenario.

Concluding remarks

We have explored the Markov state information structure in MJLS leading to new filters whose estimation error and complexity lie in the between the standard Kalman filter and the LMMSE, and establish a trade-off between performance and computational burden. This allows one to explore the computational resources at hands more deeply, seeking for the best possible estimates. As illustrated in the example, the new filters can provide competitive alternatives to the existing ones. We note that, for a given plant, the computational burden depends only on N C , see the formulas given in Section 4.1. Then, a relevant question is how to select the Markov states to form each cluster so that the estimation error is minimized, which will be considered in future research.

Example 3 . 1

 31 gives the conditional error covariance matrix, e.g. it coincides with the filtering Riccati equation of the Kalman filter when N = N C (as many clusters as Markov states). We illustrate the notation introduced in (6) on a simple example. Consider the MJLS (1) with

Fig. 1 .

 1 Fig. 1. Mean square error E(∥ x10∥ 2 ) versus number of clusters NC for every cluster configuration for the system with data given in (22).

Figure 1 :

 1 Figure 1: Mean square error E( x102 ) versus number of clusters N C for every cluster configuration for the system with data given in (22).

6 Fig. 2 .

 62 Fig. 2. Mean square error E(∥ x10∥ 2 ) versus number of clusters NC for every cluster configuration for the system with the modified data (23).

Figure 2 :

 2 Figure 2: Mean square error E( x102 ) versus number of clusters N C for every cluster configuration for the system with data given in (23).

General recursive linear estimators are briefly addressed in Section 4.4.

This work was supported by Inria associate team CDSS and ANR grant Piece, FAPESP and CNPq.