Gabriele Ciaramella

Martin J Gander

Laurence Halpern

Julien Salomon

METHODS OF REFLECTIONS: RELATIONS WITH SCHWARZ METHODS AND CLASSICAL STATIONARY ITERATIONS, SCALABILITY AND PRECONDITIONING

Keywords: 1991 Mathematics Subject Classification. 65N55, 65F10, 65N38, 35J05, 35J57, 35J25 Methods of reflections, Domain decomposition methods, Schwarz methods, coarse correction, two-level methods, substructured methods, elliptic PDE, Laplace equation

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

We start by briefly tracing the invention of the method of reflections and its historical development, for a more detailed treatment, see [START_REF] Ciaramella | Review of the methods of reflections[END_REF]. The fundamental idea for the method of reflections can already be found in the book of Murphy [39, page 93] from 1833 under the name "principle of successive influences", which clearly indicates that the method of reflections is a method of successive approximations. Lamb used in 1906 a similar approach for Laplace's equation [33, page 122]. A related method in the work of Lorentz [37, page 29] was later called "Spiegelungsmethode" (method of reflections) in [25, page 928]. The method of reflections itself was then presented in concrete mathematical notation by Smoluchowski in 1911 with the goal to understand how the motion of a sphere in a viscous fluid is influenced by the presence or motion of another sphere [START_REF] Smoluchowski | Über die Wechselwirkung von Kugeln, die sich in einer zähen Flüssigkeit bewegen[END_REF]. In the case of the Stokes equation, the analytical solution for one sphere is already a series expansion, and thus Smoluchowski assumes that the radii of the two spheres are small compared to their distance, and then uses the method of reflections to compute a series expansion of the coupled solution up to some order in the inverse distance of the spheres. The method is thus a direct method to obtain a series solution in the inverse of the distance, up to some order. Smoluchowski then generalizes the method of reflections to the case of more than two spheres, leading to similar series approximations. In 1934, Golusin introduced a parallel method of reflections for Laplace's equation for J objects [START_REF] Golusin | Auflösung einiger ebenen Grundaufaben der mathematischen Physik im Fall der Laplaceschen Gleichung und mehrfachzusammenhängender Gebiete, die durch Kreise begrenzt sind[END_REF][START_REF] Golusin | Auflösung des dreidimensionalen Dirichleteschen Problems für die Laplacesche Gleichung und Gebiete, die durch endlich viele Sphären ohne gemeinsame Punkte begrenzt sind[END_REF], and derived a condition for its convergence, which indicates that the parallel method of reflections in the case of more than two objects converges only under certain additional conditions. Golusin conjectured that these conditions depend on the distance between the objects and their radius. In 1942, Burgers [START_REF] Burgers | On the influence of the concentration of a suspension upon the sedimentation velocity (in particular for a suspension of spherical particles)[END_REF][START_REF] Burgers | Hydrodynamics. -On the influence of the concentration of a suspension upon the sedimentation velocity[END_REF] investigated the influence of the concentration of spherical particles on the sedimentation velocity for the Stokes equation, mentioning the work of Smoluchowski, but without describing precisely an algorithm, and using to a large extend physical intuition. In 1959, Kynch presented for the Stokes equation the idea of simply summing two solutions corresponding to two particles alone in [32, page 197], under the assumption that the distance between their centers is again large. This could be interpreted as a parallel method of reflections, where the separate contributions are also summed, but again, no general algorithm is given. Happel and Brenner explained in 1983 a different parallel version of the method of reflections which alternates between one fixed object and the group of all the others treated in parallel, see [START_REF] Happel | Low Reynolds number hydrodynamics with special applications to particulate media, volume 1 of Mechanics of fluids and transport processes[END_REF], with the goal to increase the order of approximation of the expansion of the solution in a neighborhood of a given object. Their method has to be applied (independently) for each particle.

Luke gave then in 1989 a first convergence analysis for the alternating method of reflections applied to the Stokes equation [START_REF] Luke | Convergence of a multiple reflection method for calculating Stokes flow in a suspension[END_REF], using a variational characterization of the method based on projections, similar in style of one of the classical convergence analyses of the Schwarz method given by Lions in the first of his three seminal papers [START_REF] Lions | On the Schwarz alternating method[END_REF]. Kim and Karrila dedicated an entire chapter in their book from 1991 to the parallel method of reflections for the Stokes problem [START_REF] Kim | Microhydrodynamics: principles and selected applications[END_REF]. The method is first motivated like in [START_REF] Kynch | The slow motion of two or more spheres through a viscous fluid[END_REF] by just summing two one-particle solutions, and only the first terms in the series expansion are obtained. Dhont also dedicated a special section to the alternating method of reflections for the Stokes equation in 1996 [START_REF] Dhont | An introduction to dynamics of colloids[END_REF]Section 5.12]. The case of two objects is first treated, and then an extension to three objects is given, where Dhont goes cyclically through the three object in the algorithm. Balabane proved in 1997 convergence of the alternating method of reflections for the Helmholtz equation in unbounded domains in [START_REF] Balabane | Décomposition de domaine pour un calcul hybride de l'équation de Helmholtz[END_REF], and generalized his results to the parallel method of reflections in [START_REF] Balabane | Boundary decomposition for Helmholtz and Maxwell equations 1: disjoint subscatterers[END_REF]. These convergence results are valid however only in low frequency regimes. In 2001, Ichiki and Brady presented the parallel method of reflections [START_REF] Ichiki | Many-body effects and matrix inversion in low-Reynolds-number hydrodynamics[END_REF] for Stokes type problems. They started with the two particle case, and then gave a generalization just summing all contributions that were computed in parallel. They presented this iterative approach also in matrix form, relating it to a stationary iteration based on a matrix splitting. By numerical experiments, they showed that the method does not converge for three particles, if the separation distance of the particles is not large enough. They thus proposed to use the method as a preconditioner for a Krylov method. Traytak posed in 2006 in a short note directly the parallel method of reflections for N objects, written in PDE form for Laplace's equation [44, Section 2], and then used a theorem proved by Golusin [START_REF] Golusin | Auflösung des dreidimensionalen Dirichleteschen Problems für die Laplacesche Gleichung und Gebiete, die durch endlich viele Sphären ohne gemeinsame Punkte begrenzt sind[END_REF] to derive sufficient conditions for the convergence based on the distances between the objects. More recently, Höfer and Velázquez used the parallel method of reflections also as an analytic tool to prove homogenization results [START_REF] Höfer | The method of reflections, homogenization and screening for Poisson and Stokes equations in perforated domains[END_REF] (see also [START_REF] Jabin | Identification of the dilute regime in particle sedimentation[END_REF]), and they modified the usual parallel method by adding different weighting coefficients. Since they were interested in the theoretical case of an infinite number of objects, they could not use an alternating method. Laurent, Legendre and Salomon studied the alternating and parallel methods of reflections in [START_REF] Laurent | On the method of reflections[END_REF] for various types of boundary conditions, introducing also an averaged version of the parallel method. They proved convergence based on the alternating projection method in Hilbert spaces, see for example [START_REF] Reich | The optimal error bound for the method of simultaneous projections[END_REF], and also using techniques like in [START_REF] Balabane | Décomposition de domaine pour un calcul hybride de l'équation de Helmholtz[END_REF][START_REF] Balabane | Boundary decomposition for Helmholtz and Maxwell equations 1: disjoint subscatterers[END_REF].

So there are two main variants of the method of reflections: the alternating one and the parallel one. There are also two different approaches to analyze the convergence of the method of reflections: first, people worked on direct estimates performed on the single/double layer formulation of the boundary value problems involved in the iterations, see [START_REF] Golusin | Auflösung des dreidimensionalen Dirichleteschen Problems für die Laplacesche Gleichung und Gebiete, die durch endlich viele Sphären ohne gemeinsame Punkte begrenzt sind[END_REF][START_REF] Traytak | Convergence of a reflection method for diffusion-controlled reactions on static sinks[END_REF][START_REF] Balabane | Décomposition de domaine pour un calcul hybride de l'équation de Helmholtz[END_REF][START_REF] Balabane | Boundary decomposition for Helmholtz and Maxwell equations 1: disjoint subscatterers[END_REF]. There is however also the interpretation of the method as alternating projections in Hilbert spaces, see [START_REF] Luke | Convergence of a multiple reflection method for calculating Stokes flow in a suspension[END_REF][START_REF] Laurent | On the method of reflections[END_REF]. In the case of orthogonal projections this interpretation leads to convergence estimates.

When we started studying methods of reflections more than three years ago, we thought immediately that the methods must be intimately related to the Schwarz domain decomposition methods. This intuition was confirmed when we studied the literature and found that analysis techniques based on projections were already used to study methods of reflections like for Schwarz methods, and such a possible relation was even mentioned in the literature ("This paper considers a reflection method in the spirit of Schwarz's alternating procedure" [START_REF] Luke | Convergence of a multiple reflection method for calculating Stokes flow in a suspension[END_REF]). We will show here that the methods of reflections can be indeed identified with Schwarz domain decomposition methods in the case of two particles/subdomains, which leads to a new convergence proof with sharp convergence estimates using maximum principle techniques. In the case of many particles/subdomains however the methods of reflections are fundamentally different from Schwarz domain decomposition methods. The main difference is that in domain decomposition, the interface data is to be determined by the subdomain iteration, whereas in the method of reflections, the interface data is given on the boundary of the particles by the problem. Substructuring techniques from domain decomposition allow us however to reformulate the methods of reflections in the form of standard block Gauss-Seidel and block Jacobi iterations in the traces for integral formulations. This reformulation leads to a new understanding of the method of reflections also for many particles, and reveals shortcomings, especially in the parallel variant. It also allows us to develop new coarse space corrections, again using domain decomposition techniques, in order to obtain scalable methods of reflections when the number of particles grows.

Our paper is structured as follows: we start by presenting in Section 2 the class of Laplace problems on perforated domains we will use to understand the methods of reflections. We then present the different forms of the methods of reflections: in Section 3 the alternating method of reflections, in Section 4 the parallel method of reflections, in Section 5 some other variants of the method of reflections. We give in each case reformulations of the methods using domain decomposition techniques, which allows us to study their convergence properties and give new convergence proofs and convergence estimates. In Section 6, we then give a scalability analysis, and also propose for the first time for the methods of reflections a coarse space

Ω O 1 O 2 O 3 0 a1 a2 a3 b1 b2 b3 1 ϕ a 1 ϕ b 1 ϕ a 2 ϕ b 2 ϕ a 3 ϕ b 3 Figure 1.
Left: example of a domain Ω (a disk) with three objects O j (holes). Right: example of the functions ϕ a j and ϕ b j for J = 3.

correction to make them scalable in the case when the number of particles becomes large. We show numerical experiments in Section 7, followed by a brief conclusion in Section 8.

The Laplace problem in perforated domains

We consider a Laplace problem where J objects O j (holes) are present in a simply connected domain Ω ⊂ R n ; for an example, see Figure 1 on the left. We assume that Ω and O j are simply connected with sufficiently smooth boundaries ∂Ω and ∂O j , e.g., of class C 1 , and that

(1) O j ∩ O k = ∅ for any j = k.
Our goal is finding the (weak) solution u ∈ H 1 0 (Ω) to the problem ∆u = 0 in Ω \ ∪ j O j , u = 0 on ∂Ω, u = g j on ∂O j for j = 1, . . . , J,

where we assume that the functions g j are bounded on ∂O j and in H 1/2 (∂O j ) for j = 1, . . . , J. This problem is well-posed and uniquely solved by u ∈ H 1 0 (Ω \ ∪ j O j ) ∩ C 2 (Ω); see, e.g., [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF][START_REF] Kress | Linear Integral Equations[END_REF]. In order to introduce the method of reflections, we modify (2) by extending harmonically the solution u into the objects O j . We define Ω := Ω \ ∪ J j=1 ∂O j and then (2) becomes ∆u = 0 in Ω, u = 0 on ∂Ω, u = g j on ∂O j for j = 1, . . . , J.

(

The solution u to (3) belongs to the Hilbert space (4)

H := v ∈ H 1 0 (Ω) : ∆v = 0 in Ω , endowed with the H 1 0 scalar product •, • H 1 0 := ∇(•), ∇(•) L 2 ,
and by Weyl's theorem [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF], we have that H ⊂ C ∞ ( Ω).

We introduce the subspaces H j of H defined by

(5)

H j := v ∈ H : ∆v = 0 in Ω \ ∂O j .
If Ω is unbounded, one needs to consider the completion of H 1 0 for the semi-norm induced by •, • H 1 0 , for instance weighted Sobolev spaces, see [START_REF] Deny | Les espaces du type de Beppo Levi[END_REF][START_REF] Hanouzet | Espaces de Sobolev avec poids. Application au problème de Dirichlet dans un demi espace[END_REF], but similar results could be obtained.

It has been proved in [START_REF] Laurent | On the method of reflections[END_REF] that H = J j=1 H j . For our purposes, we need to refine this result. We have the following lemma. Lemma 1. Under the Assumption (1), we have that H = J j=1 H j , which means that the subspaces H j are linearly independent, that is

(H 1 + • • • + H j-1 ) ∩ H j = {0} for any j ≥ 2.
Proof. Consider any function v in H j ∩ H k for any j = k. Then v is harmonic everywhere and zero on ∂Ω, thus v = 0.

This and the fact that H = J j=1 H j proved in [START_REF] Laurent | On the method of reflections[END_REF] imply the claim. To explain Lemma 1 we consider a one-dimensional example, and then show that if Assumption (1) is not satisfied, then the result does not hold in general. The domain Ω is the interval (0, 1) and the objects are subintervals

O j = (a j , b j ) for 1 ≤ j ≤ J such that 0 < a 1 < b 1 < a 2 < • • • < a J < b J < 1.
In this case, the space H j is spanned by the hat-functions

ϕ a j (x) :=      x aj if x ∈ [0, a j ], bj -x bj -aj if x ∈ [a j , b j ], 0 if x ∈ [b j , 1], ϕ b j (x) :=      0 if x ∈ [0, a j ],
x-aj bj -aj

if x ∈ [a j , b j ], 1-x 1-bj if x ∈ [b j , 1],
that have value 1 at a j and b j and are supported in [0, b j ] and [a j , 1], for an example, see Figure 1 (right). Notice that dim H j = 2. Therefore, since all the points a j and b j are distinct, we have that all the functions ϕ a j and ϕ b j are linearly independent,

H = span{ϕ a 1 , ϕ b 1 , . . . , ϕ a J , ϕ b J }, and dim H = 2J. The hypothesis a j < b j < a j+1 < b j+1 then clearly implies that (H 1 + • • • + H k-1 ) ∩ H k = {0} for 2 ≤ k ≤ J,
which is the result of Lemma 1. On the other hand, if we assume that two objects are not distinct, that is, e.g., b j = a j+1 , then one can verify that c 0 ϕ a j+1 + c

1 ϕ b j+1 + c 2 ϕ a j + c 3 ϕ b j = 0, with c 0 = 1, c 1 = 1-bj+1 1-bj , c 2 = - aj aj+1
, and c 3 = -1. Hence the functions ϕ a j , ϕ b j , ϕ a j+1 , and ϕ b j+1 are not linearly independent and dim H = 2J -1 < 2J = J j=1 dim H j . Hence H cannot be written as the direct sum of the subspaces H j . Iterative methods suitable for solving problem (2)-(3) sequentially or in parallel are the methods of reflections (MR), which we will present and study in the following sections.

The alternating method of reflections

In this section, we study the alternating method of reflections (AltMR). We first recall in Section 3.1 the classical AltMR formulation present in the literature. In Section 3.2, we introduce a new, equivalent (volume) formulation of the AltMR. This is then used in Section 3.3 where we present the AltMR in substructured form, and study its relations with the classical block Gauss-Seidel iterative method. Finally, in Section 3.4 we study the relations between the AltMR and the alternating Schwarz method.

3.1.

The classical AltMR formulations. The alternating method of reflections was invented for two objects by Smoluchowski in 1911 [START_REF] Smoluchowski | Über die Wechselwirkung von Kugeln, die sich in einer zähen Flüssigkeit bewegen[END_REF], and then extended to the general case of J objects by Luke [START_REF] Luke | Convergence of a multiple reflection method for calculating Stokes flow in a suspension[END_REF]. To understand this method, we begin by quoting a sentence from Luke to apply it to our problem: "The strategy of the method of reflections is to repeatedly correct the boundary values on the various objects. The reflection procedure is formalized as follows": it starts with u 0 ∈ H, that is [START_REF] Chaouqui | On the scalability of classical one-level domain-decomposition methods[END_REF] ∆u 0 = 0 in Ω \ ∪ j ∂O j , u 0 = 0 on ∂Ω.

Note that u 0 does not satisfy the boundary conditions of (2) on the holes, which would require u 0 = g j on ∂O j for j = 1, . . . , J.

The sequence of approximate solutions {u k } k∈N is defined as follows: given u 0 , one computes for iteration index k = 1, 2, . . . from u k-1

u k-1+ 1 J = u k-1 + d k 1 , u k-1+ 2 J = u k-1 + d k 1 + d k 2 , u k-1+ 3 J = u k-1 + d k 1 + d k 2 + d k 3 , . . .
and then obtains the new approximation ( 7)

u k = u k-1 + J j=1 d k j ,
where d k j ∈ H j is computed in such a way that, when added to u k-1+ j-1 J , the value on the boundary ∂O j is corrected to its given boundary value g j , that is

u k-1+ j J = u k-1+ j-1 J
+ d k j = g j on ∂O j . This means that d k j , for j = 1, . . . , J, must be the solution to ∆d k j = 0 in Ω \ ∂O j , d k j = 0 on ∂Ω,

d k j = g j -u k-1+ j-1 J = g j -u k-1 - j-1 =1 d k on ∂O j . (8)
A simple and intuitive explanation of the iteration described in ( 7)-( 8) can be found in [START_REF] Ciaramella | Review of the methods of reflections[END_REF].

In [START_REF] Laurent | On the method of reflections[END_REF][START_REF] Balabane | Décomposition de domaine pour un calcul hybride de l'équation de Helmholtz[END_REF][START_REF] Ciaramella | Review of the methods of reflections[END_REF] an equivalent form of ( 8) is presented, which we derive now for completeness: we consider the boundary condition of (8) on ∂O j , for j = 1, . . . , J, and manipulate it as follows:

d k j = g j -u k-1 - j-1 =1 d k = g j -u k-2 - J =1 d k-1 - j-1 =1 d k = g j -u k-2 - j =1 d k-1 - J =j+1 d k-1 - j-1 =1 d k . Now, we notice that g j -u k-2 - j =1 d k-1 = 0, since d k-1 j solves at iteration k -1 problem (8), that is ∆d k-1 j = 0 in Ω \ ∂O j , d k-1 j = 0 on ∂Ω, d k-1 j = g j -u k-2 - j-1 =1 d k-1 on ∂O j .
Therefore, we have obtained that d k j can be expressed on ∂O j as combinations of other differences only, (9)

d k j = - J =j+1 d k-1 - j-1 =1 d k ,
and the explicit dependence on g j and u k-1 disappeared. Hence, iteration [START_REF] Ciaramella | Analysis of the parallel Schwarz method for growing chains of fixed-sized subdomains: Part II[END_REF] becomes

∆d k j = 0 in Ω \ ∂O j , d k j = 0 on ∂Ω, d k j = - j-1 =1 d k - J =j+1 d k-1 on ∂O j , (10)
for j = 1, . . . , J, that is the form of the sequential method of reflections presented in [START_REF] Laurent | On the method of reflections[END_REF][START_REF] Balabane | Décomposition de domaine pour un calcul hybride de l'équation de Helmholtz[END_REF]. Now the sequences {d k j } k∈N + have to be initialized for all j. To this end, it is sufficient, for example, to consider (8) for k = 1:

∆d 1 j = 0 in Ω \ ∂O j , d 1 j = 0 on ∂Ω, d 1 j = g j - j-1 =1 d 1 -u 0 on ∂O j , (11)
for j = 1, . . . , J.

Notice that the AltMR iteration in the form (8) cannot be initialized with an arbitrary function: the initial guess u 0 must be in H (recall (6)). A choice of a function u 0 / ∈ H will produce a sequence {u k } k that does not necessarily converge to the solution u to (2). In fact, (8) produces corrections d k j which are harmonic in Ω \ ∪ J j=1 ∂O j . Hence, if u 0 can be decomposed into the sum of harmonic and non-harmonic components, namely u 0 = u 0 harm + u 0 non-harm , then only u 0 harm is corrected by the MR procedure. The AltMR as a program is given in Algorithm 1.

Algorithm 1 Alternating Method of Reflections (AltMR)

Input: K (maximum number of iterations), tol (tolerance).

1: Set u 0 ∈ H and k = 1.

2: for j = 1:J do

3:

Initialize d 1 j solving problem [START_REF] Deny | Les espaces du type de Beppo Levi[END_REF]. 4: end for 5: Compute u 1 = u 0 + J j=1 d 1 j . 6: while k < K and u k -u k-1 > tol do Update k = k + 1.

8:

for j = 1:J do Compute the approximation u k = u k-1 + J j=1 d k j . 12: end while 3.2. A new formulation of the AltMR. In this section, we formulate the AltMR in a third equivalent way. To do so, we recall the d k j that solve (8), define J sequences {v k j } k∈N for any j ∈ {1, . . . , J} as ( 12)

v k j := v k-1 j + d k j with v 0 j ∈ H j ,
and u 0 := J j=1 v 0 j (recall Lemma 1). Notice that since d k j ∈ H j , v k j ∈ H j as well. Now, we notice that [START_REF] Ciaramella | Analysis of the parallel Schwarz method for growing chains of fixed-sized subdomains: Part I[END_REF] implies that u k = u 0 + J j=1 k n=0 d n j , which allows us to compute

u k = u 0 + J j=1 k n=1 d n j = u 0 + J j=1 k n=0 v n j -v n-1 j = u 0 + J j=1 v k j -v 0 j ,
and recalling that

J j=1 v 0 j = u 0 we obtain (13) u k = J j=1 v k j .
Notice that equation ( 13) allows us to express the approximate solution at the kth iteration as a decomposition (unique according to Lemma 1) into J components, each of them being associated with an object and belonging to one of the spaces H j . Starting from (12), using [START_REF] Efstathiou | Why restricted additive Schwarz converges faster than additive Schwarz[END_REF] and the boundary condition in [START_REF] Ciaramella | Analysis of the parallel Schwarz method for growing chains of fixed-sized subdomains: Part II[END_REF], that is

d k j = g j -u k-1 - j-1 =1 d k on ∂O j ,
we have on ∂O j that

v k j = v k-1 j + d k j = v k-1 j + g j -u k-1 - j-1 =1 v k -v k-1 = g j + v k-1 j - J =1 v k-1 - j-1 =1 v k + j-1 =1 v k-1 = g j - J =j+1 v k-1 - j-1 =1 v k . Recalling that v k j ∈ H j , we obtain that v k j is solution to ∆v k j = 0 in Ω \ ∂O j , v k j = 0 on ∂Ω, v k j = g j - J =j+1 v k-1 - j-1 =1 v k on ∂O j . ( 14
)
The equivalence between ( 14) and ( 10)- [START_REF] Deny | Les espaces du type de Beppo Levi[END_REF] is proved in the following theorem.

Theorem 1. Consider the sequences {u k } k∈N ⊂ H, {d k j } k∈N + with d k j ∈ H j , and {v k j } k∈N with v k j ∈ H j for j = 1, . . . , J, and assume that v k j = v k-1 j + d k j for k ≥ 1 and j = 1, . . . , J. Then the following statements are equivalent:

(a) the d k j solve (10)- [START_REF] Deny | Les espaces du type de Beppo Levi[END_REF] and

u k = u k-1 + J j=1 d n j = k n=1 J j=1 d n j . (b) the v k j solve (14) and u k = J j=1 v k j .
Proof. The implication (a) ⇒ (b) is proved in the discussion before the theorem. Hence, we need to show that (b) ⇒ (a). First, using

v k j = v k-1 j + d k j we get (15) u k = J j=1 v k j = u 0 + k n=1 J j=1 v n j -v n-1 j = u 0 + k n=1 J j=1 d n j = u k-1 + J j=1 d k j .
Notice that

d k j = v k j -v k-1 j . Clearly, d k j ∈ H j because v k j , v k-1 j ∈ H j .
Hence we have to prove that d k j satisfies the transmission condition of [START_REF] Ciaramella | Review of the methods of reflections[END_REF]. For k = 1, we use that d 1 j = v 1 j -v 0 j and the boundary condition in [START_REF] Gander | Optimized Schwarz methods[END_REF] for k = 1 to write

d 1 j = v 1 j -v 0 j = g j - j-1 =1 v 1 j - J =j v 0 j = g j - j-1 =1 d 1 j -u 0 .
In the case k > 1, we write the transmission condition on ∂O j of ( 14) for k and k -1:

v k j = g j - j-1 =1 v k - J =j+1 v k-1 , v k-1 j = g j - j-1 =1 v k-1 - J =j+1 v k-2 .
Subtracting term by term these equations, we get the transmission condition of [START_REF] Ciaramella | Review of the methods of reflections[END_REF].

Using [START_REF] Gander | Optimized Schwarz methods[END_REF], we now rewrite Algorithm 1 in terms of v k j to obtain Algorithm 2.

Algorithm 2 Alternating Method of Reflections (AltMR)

Input: K (maximum number of iterations), tol (tolerance).

1: Set v 0 j ∈ Hj for j = 1, . . . , J, and k = 1. 2: for j = 1:J do

3:

Compute v 1 j solving problem [START_REF] Gander | Optimized Schwarz methods[END_REF]. 4: end for 5: Compute the approximation u 1 = J j=1 v 1 j . 6: while k < K and u k -u k-1 > tol do 7:

Update k = k + 1.

8:

for j = 1:J do 9:

Compute v k j solving problem [START_REF] Gander | Optimized Schwarz methods[END_REF]. Compute the approximation u k = J j=1 v k j . 12: end while

We next show that this new formulation of the AltMR in terms of the decomposition functions v k j allows us to easily obtain a substructured formulation of the AltMR procedure.

3.3. Substructured AltMR as block Gauss-Seidel method. In this section, we first write the AltMR in substructured form and then show that it can be interpreted as a block Gauss-Seidel method. To do so, we recall Lemma 1, which plays a very important role in our interpretation of the method of reflections, since it ensures that the unique solution u ∈ H to (3) can be uniquely decomposed as

v 1 + v 2 + • • • + v J with v j ∈ H j .
For each of these v j , let us denote by g j the corresponding Dirichlet trace on ∂O j so that, because of Definition (5), we have

∆v j = 0 in Ω \ ∂O j , v j = 0 on ∂Ω, v j = g j on ∂O j . (16)
Now, we introduce the operator G j that maps Dirichlet boundary data into functions that are harmonic in Ω \ ∂O j , and such that the solution to ( 16) can be written as v j = G j ( g j ). The operator G j can be explicitly expressed in terms of Green's representation formulas, and, in particular, as the composition of a single-layer potential integral operator and Dirichlet-to-Neumann operators: for ḡ ∈ H 1/2 (∂O j ), we have

G j (ḡ)(x) := ∂Oj G(x, y) DtN j,e (ḡ) + DtN j,i (ḡ) (y)ds(y),
where G(x, y) is the Green's function associated to the problem, DtN j,i : H 1/2 (∂O j ) → H -1/2 (∂O j ) and DtN j,e : H 1/2 (∂O j ) → H -1/2 (∂O j ). These Dirichlet-to-Neumann operators are obtained by first solving the interior and exterior problems for O j ,

∆v e = 0 in Ω \ O j , v e = ḡ on ∂O j , v e = 0 on ∂Ω, ∆v i = 0 in O j , v i = ḡ on ∂O j ,
and then extracting the Neumann trace on ∂O j of their solutions, that is DtN j,e (ḡ) :=

∂ n v e (ḡ)| ∂Oj and DtN j,i (ḡ) := -∂ n v i (ḡ)| ∂Oj , with ∂ n = n • ∇ the unit outward nor- mal derivative with respect to Ω \ O j .
In addition, we also need the trace operator

τ j : H 1 (Ω) → H 1/2 (∂O j ) such that τ j v is the Dirichlet trace of v on ∂O j , τ j v := v| ∂Oj .
With this framework, we can rewrite problem (2) in integral form: we look for boundary data g j , for j = 1, . . . , J, such that the function u := J j=1 G j ( g j ) satisfies [START_REF] Gander | A new coarse grid correction for RAS/AS[END_REF] τ

j u = τ j   J j=1 G j ( g j )   = g j , on ∂O j for j = 1, . . . , J.
Equation ( 17) is equivalent to the linear system (18)

A g = g, where (19)
A :=        I 1 τ 1 G 2 τ 1 G 3 • • • τ 1 G J τ 2 G 1 I 2 τ 2 G 3 • • • τ 2 G J τ 3 G 1 τ 3 G 2 I 3 • • • τ 3 G J . . . . . . . . . . . . . . . τ J G 1 τ J G 2 τ J G 3 • • • I J        , g :=        g 1 g 2 g 3 . . . g J       
, and g :=

       g 1 g 2 g 3 . . . g J       
, and I j are identity operators, and we used the fact that τ j G j = I j .

The following proposition shows that the linear system ( 18) is an equivalent formulation of the problem (3).

Proposition 1. The system (18) is equivalent to (3). Moreover, the solution u to (3) is uniquely decomposed as u = J j=1 G j ( g j ), where g := g 1 • • • g J is the solution to [START_REF] Gander | SHEM: An optimal coarse space for RAS and its multiscale approximation[END_REF].

Proof. The equivalence between (3) and [START_REF] Gander | SHEM: An optimal coarse space for RAS and its multiscale approximation[END_REF] follows by Lemma 1 and the fact that each component v j of u is uniquely determined by the Dirichlet boundary data g j on ∂O j . The second statement is obtained by Lemma 1 and the following argument: from [START_REF] Gander | A new coarse grid correction for RAS/AS[END_REF] we have that τ j u = g j on ∂O j for j = 1, . . . , J, and we have that ∆u = ∆ J j=1 G j ( g j ) = 0 on Ω \ ∪ j ∂O j .

Remark 1. The variables g j are not necessarily of Dirichlet type: one could formulate the problem with, e.g., Neumann data. In this case the definition of G j would not require Dirichlet-to-Neumann operators, and τ j G j would be a Neumannto-Dirichlet operator.

The linear system (18) can be solved by classical iterative methods as discussed in the next subsections. For this purpose, we consider the splitting

A = D + L + U , where D := diag(A) = I, L :=        0 0 0 • • • 0 τ 2 G 1 0 0 • • • 0 τ 3 G 1 τ 3 G 2 0 • • • 0 . . . . . . . . . . . . . . . τ J G 1 τ J G 2 τ J G 3 • • • 0        , U :=        0 τ 1 G 2 τ 1 G 3 • • • τ 1 G J 0 0 τ 2 G 3 • • • τ 2 G J 0 0 0 • • • τ 3 G J . . . . . . . . . . . . . . . 0 0 0 • • • 0        .
This decomposition leads to the classical iterative methods based on the splitting A = M -N , and we can write the iterative methods in the standard form

M g k+1 = N g k + g and in the difference form Mδ δ δ k+1 = Nδ δ δ k , with δ δ δ k := g k -g k-1
. The convergence of this class of methods is related to the iteration operator G := M -1 N .

Let us now consider the block Gauss-Seidel method for the solution of system [START_REF] Gander | SHEM: An optimal coarse space for RAS and its multiscale approximation[END_REF], which is obtained via the splitting A = M -N with M := D+L and N := -U , [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF] (D + L) g k+1 = -U g k + g, and the difference form is given by

(21) (D + L)δ δ δ k+1 = -Uδ δ δ k .
More explicitly, the standard form is

(22)        I 1 0 0 ••• 0 τ 2 G 1 I 2 0 ••• 0 τ 3 G 1 τ 3 G 2 I 3 ••• 0 . . . . . . . . . . . . . . . τ J G 1 τ J G 2 τ J G 3 ••• I J               g k+1 1 g k+1 2 g k+1 3 . . . g k+1 J        =        g 1 g 2 g 3 . . . g J        -        0 τ 1 G 2 τ 1 G 3 ••• τ 1 G J 0 0 τ 2 G 3 ••• τ 2 G J 0 0 0 ••• τ 3 G J . . . . . . . . . . . . . . . 0 0 0 ••• 0               g k 1 g k 2 g k 3 . . . g k J        . Denoting the components of the difference δ δ δ k by δ k j := g k j -g k-1 j
, we can write the difference relation [START_REF] Golusin | Auflösung des dreidimensionalen Dirichleteschen Problems für die Laplacesche Gleichung und Gebiete, die durch endlich viele Sphären ohne gemeinsame Punkte begrenzt sind[END_REF] in the explicit form

(23)        I 1 0 0 ••• 0 τ 2 G 1 I 2 0 ••• 0 τ 3 G 1 τ 3 G 2 I 3 ••• 0 . . . . . . . . . . . . . . . τ J G 1 τ J G 2 τ J G 3 ••• I J               δ k+1 1 δ k+1 2 δ k+1 3 . . . δ k+1 J        =-        0 τ 1 G 2 τ 1 G 3 ••• τ 1 G J 0 0 τ 2 G 3 ••• τ 2 G J 0 0 0 ••• τ 3 G J . . . . . . . . . . . . . . . 0 0 0 ••• 0               δ k 1 δ k 2 δ k 3 . . . δ k J       
. Now we show that ( 20) and ( 21) are equivalent to ( 14) and [START_REF] Ciaramella | Review of the methods of reflections[END_REF].

Theorem 2. Assume that v 0 ∈ H and g 0 j = τ j v 0 . Then the AltMR methods ( 14) and (10) are equivalent to [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF] and [START_REF] Golusin | Auflösung des dreidimensionalen Dirichleteschen Problems für die Laplacesche Gleichung und Gebiete, die durch endlich viele Sphären ohne gemeinsame Punkte begrenzt sind[END_REF].

Proof. Consider the transmission condition of [START_REF] Gander | Optimized Schwarz methods[END_REF], that is

(24) τ j v k j = g j - J =j+1 τ j v k-1 - j-1 =1 τ j v k ,
for j = 1, . . . , J. We denote by g k j = τ j v k j and, recalling that v k-1 j solves ( 14), we can use the operator G j defined above to write

v k-1 j = G j (τ j v k-1 j ) = G j ( g k-1 j ).
Therefore, ( 24) can be equivalently written as [START_REF] Hengstenberg | Messen und Regeln in der chemischen Technik[END_REF] g k j = g j -

j-1 =1 τ j G j ( g k ) - J =j+1 τ j G j ( g k-1
).

If g 0 j = g 0 j , then equation ( 25) is the equivalent pointwise form of ( 20)-( 22) with g k j = g k j . In a similar way, using the last equation in [START_REF] Ciaramella | Review of the methods of reflections[END_REF], we obtain the formula

δ k j = - j-1 =1 τ j G j (δ k ) - J =j+1 τ j G j (δ k-1 ),
which is equivalent to ( 21)-( 23).

3.4.

Analogy with the alternating Schwarz method. In this section, we study the relations between the AltMR and the alternating Schwarz method (AltSM). In particular, we prove that if J = 2, the AltMR is equivalent to the AltSM in the sense of Theorem 3. This identification does not hold in general for J > 2.

To do so, consider problem (2), and define for j = 1, . . . , J, Ω j := Ω \ O j , and a smooth enough function g : Ω → R such that g |∂Oj = g j , g |∂Ω = 0, and we assume that there exists a smooth extension of u inside the objects that we denote by E j u. Then, we set

u := u -g in Ω j , E j u -g in O j .
Notice that the extension E j u must be smooth enough, e.g. twice differentiable, such that its Laplacian is well defined. Next, we define f : Ω → R as

f := ∆g in Ω j , ∆E j u + ∆g in O j .
Therefore, u solves the problem [START_REF] Höfer | The method of reflections, homogenization and screening for Poisson and Stokes equations in perforated domains[END_REF] ∆ u = f in Ω with u = 0 on ∂Ω, and by construction it satisfies u = 0 on J j=1 ∂O j . We can now apply the AltSM to solve [START_REF] Höfer | The method of reflections, homogenization and screening for Poisson and Stokes equations in perforated domains[END_REF]. To do so, we consider the domain decomposition Ω = ∪ J j=1 Ω j , and denote by u k j the approximate solution at the kth iteration on the subdomain Ω j . Starting with some initial approximation u 0 j for all j, the AltSM is defined as

∆ u k j = f in Ω j , u k j = 0 on ∂Ω, u k j = j-1 =1 α j, u k + J =j+1 α j, u k-1 on ∂O j ,
where the coefficients α j, are non-negative and such that J =1 α j, = 1; see, e.g., [START_REF] Gander | Schwarz methods over the course of time[END_REF][START_REF] Toselli | Domain Decomposition Methods: Algorithms and Theory[END_REF][START_REF] Efstathiou | Why restricted additive Schwarz converges faster than additive Schwarz[END_REF]. This definition implies that the corresponding differences

d k j := u k j -u k-1 j solve ∆ d k j = 0 in Ω j , d k j = 0 on ∂Ω, d k j = j-1 =1 α j, d k + J =j+1 α j, d k-1 on ∂O j . (27)
Next, we prove in Theorem 3 the analogy between the AltMR and the AltSM for J = 2, and afterwards we provide a counterexample to show that this analogy does not hold in general for J > 2. Proof. Since J = 2, we have that α 1,2 = α 2,1 = 1. We first prove the relation

d k 2 = d k 2 by induction.
The result is true for k = 0 by assumption. Now, we assume that the relation holds for k and prove that it is true for k + 1 as well. Recalling the transmission conditions of ( 27) and [START_REF] Ciaramella | Review of the methods of reflections[END_REF]

(for J = 2 and α 1,2 = α 2,1 = 1), it holds that (28) d k+1 1 = d k 2 = d k 2 = -d k+1 1 on ∂O 1 ,
where we used the induction hypothesis. Equation ( 28), together with the existence of unique solutions to [START_REF] Ichiki | Many-body effects and matrix inversion in low-Reynolds-number hydrodynamics[END_REF] and [START_REF] Ciaramella | Review of the methods of reflections[END_REF], implies that

d k+1 1 = -d k+1 1
on Ω 1 . Using this equality and the transmission conditions of ( 27) and [START_REF] Ciaramella | Review of the methods of reflections[END_REF], we have that

d k+1 2 = d k+1 1 = -d k+1 1 = d k+1 2 on ∂O 2 .
Hence the well-posedness of ( 27) and [START_REF] Ciaramella | Review of the methods of reflections[END_REF]

implies that d k 2 = d k 2 on Ω 2
, which is our first claim. The second relation follows directly by the first one together with [START_REF] Jabin | Identification of the dilute regime in particle sedimentation[END_REF].

The equivalence proved in Theorem 3 does not hold in general for J > 2. In fact, we now show that classical choices of Schwarz methods (like the classical AltSM) does not lead to the result proved in Theorem 3. This negative result is suggested by the transmission conditions of ( 27) and [START_REF] Ciaramella | Review of the methods of reflections[END_REF]. In fact, if we set α ,j = -1 for all j, , the transmission condition of ( 27) coincides with the ones of [START_REF] Ciaramella | Review of the methods of reflections[END_REF]. However, this is not possible in a Schwarz method framework, where the hypothesis J =1 α j, = 1, with α j, non-negative, is required. Nevertheless, we provide the following example to show that the iterates d k j and d k j are unrelated for J > 2. To do so, consider a domain Ω = (0, 1) and three holes O j = (a j , b j ) for j = 1, 2, 3 with a j = 2j-1 7 and b j = 2j 7 . The alternating Schwarz method [START_REF] Ichiki | Many-body effects and matrix inversion in low-Reynolds-number hydrodynamics[END_REF] becomes

∆ d k 1 = 0 in Ω 1 , d k 1 (0) = d k 1 (1) = 0, d k 1 ( x) = α 1,2 d k-1 2 ( x) + α 1,3 d k-1 3 ( x) for x = a 1 , b 1 , ∆ d k 2 = 0 in Ω 2 , d k 2 (0) = d k 2 (1) = 0, d k 2 ( x) = α 2,1 d k 1 ( x) + α 2,3 d k-1 3 ( x) for x = a 2 , b 2 , ∆ d k 3 = 0 in Ω 3 , d k 3 (0) = d k 3 (1) = 0, d k 3 ( x) = α 3,1 d k 1 ( x) + α 3,2 d k-1 2 ( x) for x = a 3 , b 3 .
Notice that by setting α ,j = -1, the above problems coincide with the AltMR [START_REF] Ciaramella | Review of the methods of reflections[END_REF]. The general solutions to these three problems are

d k 1 (x) = A k 1 x a1
x ∈ [0, a 1 ],

B k 1 (1-x) 1-b1 x ∈ [b 1 , 1], d k 2 (x) = A k 2 x a2
x ∈ [0, a 2 ],

B k 2 (1-x) 1-b2 x ∈ [b 2 , 1], d k 3 (x) = A k 3 x a3
x ∈ [0, a 3 ],

B k 3 (1-x) 1-b3 x ∈ [b 3 , 1],
where A k j and B k j are constants depending on the transmission conditions. Defining

v k := [A k 1 , B k 1 , A k 2 , B k 2 , A k 3 , B k 3 ]
and using the transmission conditions, we obtain the iteration relation

(I + L)v k = -U v k-1 , where L = -           0 0 0 0 0 0 0 0 0 0 0 0 0 α2,1(1-a2) (1-b1) 0 0 0 0 0 α2,1(1-b2) (1-b1) 0 0 0 0 0 α3,1(1-a3) (1-b1) 0 α3,2(1-a3) (1-b2) 0 0 0 α3,1(1-b3) (1-b1) 0 α3,2(1-b3) (1-b2) 0 0           and U = -          0 0 α1,2a1 a2 0 α1,3a1 a3
0 0 0 α1,2b1 a2 0 α1,3b1 a3 0 0 0 0 0 α2,3a2 a3 0 0 0 0 0 α2,3b2 a3 0 0 0 0 0 0 0 0 0 0 0 0 0          .
Let us define the matrix G = -(I 6 + L) -1 U , with I 6 the 6×6 identity. The iteration relation above then reads v k = Gv k-1 , where the matrix G depends on the weights α j,k . Notice that the same relation (with appropriately chosen weights in G) can be obtained for the AltMR. We consider three different cases. The first case is α 1,2 = α 2,3 = α 3,1 = 0 and α 1,3 = α 2,1 = α 3,2 = 1, which corresponds to the classical AltSM, the second case is α j,k = 1 2 for any j, k, which leads to a weighted AltSM, and the third case α j,k = -1 for any j, k, which corresponds to the AltMR. Computing explicitly the iteration matrix G is these cases, we get

G 1 =         0 0 0 0 1 5 0 0 0 0 0 2 5 0 0 0 0 0 8 25 0 0 0 0 0 6 25 0 0 0 0 0 4 25 0 0 0 0 0 2 25 0         , G 2 =         0 0 1 6 0 1 10 0 0 0 1 3 0 1 5 0 0 0 2 15 0 19 50 0 0 0 1 10 0 23 50 0 0 0 1 10 0 29 150 0 0 0 1 20 0 29 300 0         , G 3 =         0 0 -1 3 0 -1 5 0 0 0 -2 3 0 -2 5 0 0 0 8 15 0 -7 25 0 0 0 2 5 0 -14 25 0 0 0 0 0 8 15 0 0 0 0 0 4 15 0         . Their spectra are σ(G 1 ) = 4 25 , 0 , σ(G 2 ) = -(3 √ 389 -49) 300 , 49 + 3 √ 389 300 , 0 , σ(G 3 ) = 8 15
, 0 , and we thus obtain for the spectral radii

ρ(G 1 ) = 4 25 = 0.16, ρ(G 2 ) = 49 + 3 √ 389 300 ≈ 0.36, ρ(G 3 ) = 8 15 ≈ 0.53.
This shows that the iterations generated by G 1 and G 2 corresponding to alternating Schwarz methods must be different from the iterations generated by G 3 corresponding to the AltMR, and we have thus shown that there is in general no relation between the AltMR and the AltSM for J > 2.

The reader may ask himself how the iteration matrices look like for J = 2. To illustrate this, we consider a domain Ω = (0, 1) with two holes O j = (a j , b j ) for j = 1, 2 with a j = 2j-1 5 and b j = 2j 5 , and recall that α 1,2 = α 2,1 = 1. The same arguments as above allow us to obtain the iteration matrices

G AltSM =     0 0 0 0 1 3 0 0 0 0 0 2 3 0 0 0 0 0 4 9 0 0 0 0 0 2 9 0     , G AltMR =     0 0 0 0 -1 3 0 0 0 0 0 -2 3 0 0 0 0 0 4 9 0 0 0 0 0 2 9 0     ,
whose spectra coincide and are given by

σ(G AltMR ) = σ(G PSM ) = 4 9 , 0 .
It is clear that the two matrices generate similar iterates. Notice also the negative signs in G AltMR that produce d k 1 with opposite sign to d k 1 , in agreement with Theorem 3.

Convergence analysis for

J = 2 objects. Consider problem (2) for J = 2: ∆u = 0 in Ω \ ∂O 1 ∪ ∂O 2 , u = 0 on ∂Ω, u = g 1 on ∂O 1 , u = g 2 on ∂O 2 . ( 29
)
The AltMR [START_REF] Gander | Optimized Schwarz methods[END_REF] for the solution of ( 29) is

∆v k 1 = 0 in Ω \ ∂O 1 , v k 1 = 0 on ∂Ω, v k 1 = g 1 -v k-1 2 on ∂O 1 , ∆v k 2 = 0 in Ω \ ∂O 2 , v k 2 = 0 on ∂Ω, v k 2 = g 2 -v k 1 on ∂O 2 . (30)
Let us define the error at the kth iteration by e k j := v j -v k j for j = 1, . . . , J. In terms of the errors, the AltMR (30) reads

∆e k 1 = 0 in Ω \ ∂O 1 , e k 1 = 0 on ∂Ω, e k 1 = -e k-1 2 on ∂O 1 , ∆e k 2 = 0 in Ω \ ∂O 2 , e k 2 = 0 on ∂Ω, e k 2 = -e k 1 on ∂O 2 . (31)
We can now prove the following theorem using similar techniques considered in [START_REF] Ciaramella | Analysis of the parallel Schwarz method for growing chains of fixed-sized subdomains: Part II[END_REF].

Theorem 4. The AltMR for the solution to problem (29) converges geometrically, in the sense that max

Ω |e k j | ≤ ρ k AltMR max Ω |e 0 j |,
for j = 1, 2, where ρ AltMR = max ∂O1 w 2 max ∂O2 w 1 < 1, and w j solves (32) ∆w j = 0 in Ω \ ∂O j , w j = 0 on ∂Ω, w j = 1 on ∂O j .

Proof. The functions w j , for j = 1, 2 which are solutions to [START_REF] Kynch | The slow motion of two or more spheres through a viscous fluid[END_REF], satisfy because of the maximum principle that w j ≥ 0 and |e k j | ≤ w j max ∂Oj |e k j |. Now, since e k and e k 2 are harmonic functions and solve [START_REF] Kress | Linear Integral Equations[END_REF], the maximum principle allows us to write max

Ω |e k 1 | = max ∂O1 |e k 1 | = max ∂O1 |e k-1 2 | ≤ max ∂O1 w 2 max ∂O2 |e k-1 2 | = max ∂O1 w 2 max Ω |e k-1 2 |, max Ω |e k-1 2 | = max ∂O2 |e k-1 2 | = max ∂O2 |e k-1 1 | ≤ max ∂O2 w 1 max ∂O1 |e k-1 1 | = max ∂O2 w 1 max Ω |e k-1 1 |.
By combining these two estimates we get max

Ω |e k 1 | ≤ max ∂O1 w 2 max Ω |e k-1 2 | ≤ max ∂O1 w 2 max ∂O2 w 1 max Ω |e k-1 1 | = ρ AltMR max Ω |e k-1 1 |.
Since w 1 and w 2 are harmonic functions in Ω \ ∂O 1 and Ω \ ∂O 2 , Assumption (1) and the maximum principle imply that max ∂O1 w 2 < 1 and max ∂O2 w 1 < 1. Hence ρ AltMR < 1 and our proof is complete.

This theorem allows us to compute explicitly the contraction factor ρ AltMR as a function of the geometry of the domain. For example, consider a domain whose geometry is shown in Figure 2, where Ω ⊂ R 2 is the unit disc, the objects O 1 and O 2 are two discs whose centers are aligned on a straight line passing trough the center of Ω ⊂ R 2 . The two discs have the same radius r and the distance between the center of each of them and the center of Ω is denoted by δ. We can prove the following result.

Corollary 1. Consider the problem (29) defined on a domain whose geometry is depicted in Figure 2. We have that

ρ AltMR (r, δ) = log (-a-δ+r) 2 (a(r-δ)-1) 2 log (-a+δ-r) 2 (a(δ-r)-1) 2 2 ,
where Proof. To prove the result, we first need to solve problems [START_REF] Kynch | The slow motion of two or more spheres through a viscous fluid[END_REF] and compute the w j . To do so, we recall the Möbius transformation h : C → C:

(33) a = 1 + δ 2 -r 2 -r 4 + (-2 -2δ 2 ) r 2 + δ 4 -2δ 2 + 1 2δ .
(34) h(z) := z -a āz -1 ,
which maps the unit circle into itself and circles into circles. In particular, by imposing the conditions h(δ -r) = R and h(δ + r) = -R, one obtains a (real) coefficient a, as in [START_REF] Lamb | Hydrodynamics[END_REF], such that h maps Ω \ O j into an annulus; see Figure 3. We consider now the problem (35) ∆ w 1 = 0 in h(Ω) \ h(∂O 1 ), w 1 = 0 on h(∂Ω), w 1 = 1 on h(∂O 1 ).

Since [START_REF] Lions | On the Schwarz alternating method[END_REF] is radially symmetric, a separation of variables argument allows us to compute its solution,

w 1 (ζ) = log Re(ζ) 2 + Im(ζ) 2 log(R) .
By transforming back from ζ to z, we obtain

w 1 (Re(z), Im(z)) = w 1 (h(z)) = log(|h(z)|) log(R) = log(|h(z)| 2 ) log(|R| 2 ) = log z-a az-1 2 log δ-r-a a(δ-r)-1 2 ,
which solves (32) for j = 1 since h is a conformal mapping. A further simplification leads to

w 1 (x, y) = log (x-a) 2 +y 2 (ax-1) 2 +(ay) 2 log (δ-r-a) 2 (a(δ-r)-1) 2
. Now, since w 1 is symmetric with respect to the Re(ζ)-axis and decays monotonically in any radial direction, the same holds for w 1 . Hence, the maximum of w 1 along ∂O 2 is attained at (x, y) = (-δ + r, 0) and has the value

max ∂O2 w 1 = w 1 (-δ + r, 0) = log (-δ+r-a) 2 (a(-δ+r)-1) 2 log (δ-r-a) 2 (a(δ-r)-1) 2 .
Since the solution w 2 to problem (32) for j = 2 can be obtained by rotating w 1 by π around the origin, it holds that max

∂O1 w 2 = w 2 (δ -r, 0) = w 1 (-δ + r, 0) = max ∂O2 w 1 .
Recalling from Theorem 4 that ρ AltMR = max ∂O1 w 2 max ∂O2 w 1 , the claim follows.

The parallel method of reflections

The parallel method of reflections (PMR) was introduced by Golusin in 1934 [START_REF] Golusin | Auflösung des dreidimensionalen Dirichleteschen Problems für die Laplacesche Gleichung und Gebiete, die durch endlich viele Sphären ohne gemeinsame Punkte begrenzt sind[END_REF] and formally recalled by Traytak [START_REF] Traytak | Convergence of a reflection method for diffusion-controlled reactions on static sinks[END_REF] for the Laplace equation. Ichiki and Brady [START_REF] Ichiki | Many-body effects and matrix inversion in low-Reynolds-number hydrodynamics[END_REF] present exactly the parallel method of Golusin, and they state: "It is easy to extend this procedure to the N body problem by superposing distances by other particles". In practice, the parallel version is obtained by replacing in the righthand side of the boundary condition ( 9) the differences at the iteration k with the corresponding ones at the iteration k -1. Hence, problem [START_REF] Ciaramella | Review of the methods of reflections[END_REF] becomes

∆d k j = 0 in Ω \ ∂O j , d k j = 0 on ∂Ω, d k j = - J =1, =j d k-1 on ∂O j . (36)
The sequences {d k j } k∈N + are initialized by solving for each j = 1, . . . , J the problem

∆d 1 j = 0 in Ω \ ∂O j , d 1 j = 0 on ∂Ω, d 1 j = g j -u 0 on ∂O j , (37)
for some u 0 ∈ H, and the approximate solution at the kth iteration is defined by ( 38)

u k = u k-1 + J j=1 d k j .
The PMR ( 36)-( 37), as presented in [START_REF] Laurent | On the method of reflections[END_REF], leads to Algorithm 3.

Algorithm 3 Parallel Method of Reflection (PMR)

Input: K (maximum number of iterations), tol (tolerance).

1: Set u 0 ∈ H and k = 1.

2: for j = 1:J (this loop is executed in parallel) do

3:

Compute d 1 j solving problem (37). 4: end for 5: Compute the approximation u 1 = J j=1 d 1 j . 6: while k < K and u k -u k-1 > tol do Update k = k + 1.

8:

for j = 2:J (this loop is executed in parallel) do Compute the approximation u k = u k-1 + J j=1 d k j . 12: end while Now, as in Section 3, we introduce the variable v k j defined in [START_REF] Dhont | An introduction to dynamics of colloids[END_REF], and we write the PMR in terms of v k j . To do so, we first use the transmission condition of [START_REF] Lions | On the Schwarz alternating method. II. Stochastic interpretation and other properties[END_REF] and [START_REF] Luke | Convergence of a multiple reflection method for calculating Stokes flow in a suspension[END_REF] to compute on ∂O j that

d k j = - J =1, =j d k-1 = d k-1 j - J =1 d k-1 = d k-1 j -(u k-1 -u k-2 ),
which implies that

d k j + u k-1 = d k-1 j + u k-2 = • • • = d 1 j + u 0 . Recalling that d 1 j = g j -u 0 on ∂O j , we obtain (39)
d k j = g j -u k-1 on ∂O j . Therefore, we have that ( 36) is equivalent to

∆d k j = 0 in Ω \ ∂O j , d k j = 0 on ∂Ω, d k j = g j -u k-1 on ∂O j , (40)
which is the parallel version of ( 8). Now, starting with the definition (12) and using ( 39) and ( 13), we compute on ∂O j that

v k j = v k-1 j + d k j = v k-1 j + g j -u k-1 = v k-1 j + g j - J =1 v k-1 = g j - J =1, =j v k-1 . Recalling that v k j ∈ H j , we conclude that v k j solves ∆v k j = 0 in Ω \ ∂O j , v k j = 0 on ∂Ω, v k j = g j - J =j+1 v k-1 - j-1 =1 v k-1 on ∂O j . (41)
The equivalence between ( 41) and ( 36)-( 37) is proved in the following theorem.

Theorem 5. Consider the sequences {u k } k∈N ⊂ H, {d k j } k∈N + with d k j ∈ H j , and {v k j } k∈N with v k j ∈ H j for j = 1, . . . , J, and assume that v k j = v k-1 j + d k j for k ≥ 1 and j = 1, . . . , J. Then the following statements are equivalent: The PMR algorithm in terms of v k j is given in Algorithm 4.

4.1.

Substructured PMR as a block Jacobi method. In this section, we prove that the PMR is equivalent to the block Jacobi method. We recall that Ichiki and Brady [27, page 351] mention, without motivation or any rigorous argument, that "this iterative method (the PMR) is equivalent to the block Jacobi method". This sentence suggests a block Jacobi substructured interpretation of the PMR, as the block Gauss-Seidel one for the AltMR proved in Section 3.3. To see it, consider the operator A defined in Section 3. Input: K (maximum number of iteration), tol (tolerance).

1: Set v 0 j ∈ Hj for j = 1, . . . , J, and k = 1. 2: for j = 1:J (this loop is executed in parallel) do

3:

Compute v 1 j solving problem (41). 4: end for 5: Compute the approximation u 1 = J j=1 v 1 j . 6: while k < K and u k -u k-1 > tol do 7:

Update k = k + 1.

8:

for j = 1:J (this loop is executed in parallel) do

9:

Compute v k j solving problem [START_REF] Richardson | The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam[END_REF]. Compute the approximation u k = J j=1 v k j . 12: end while and in the difference form it is

(43) Dδ δ δ k+1 = -(L + U )δ δ δ k ,
where δ δ δ k := g kg k-1 , see also [START_REF] Hanouzet | Espaces de Sobolev avec poids. Application au problème de Dirichlet dans un demi espace[END_REF]. More explicitly, the standard form is

(44)        g k+1 1 g k+1 2 g k+1 3 . . . g k+1 J        =        g 1 g 2 g 3 . . . g J        -        0 τ 1 G 2 τ 1 G 3 • • • τ 1 G J τ 2 G 1 0 τ 2 G 3 • • • τ 2 G J τ 3 G 1 τ 3 G 2 0 • • • τ 3 G J . . . . . . . . . . . . . . . τ J G 1 τ J G 2 τ J G 3 • • • 0               g k 1 g k 2 g k 3 . . . g k J       
, whereas the difference form ( 43) is

(45)        δ k+1 1 δ k+1 2 δ k+1 3 . . . δ k+1 J        = -        0 τ 1 G 2 τ 1 G 3 • • • τ 1 G J τ 2 G 1 0 τ 2 G 3 • • • τ 2 G J τ 3 G 1 τ 3 G 2 0 • • • τ 3 G J . . . . . . . . . . . . . . . τ J G 1 τ J G 2 τ J G 3 • • • 0               δ k 1 δ k 2 δ k 3 . . . δ k J        .
The next theorem shows that ( 42)-( 44) and ( 43)-( 45) are equivalent to the PMR.

Theorem 6. Assume that v 0 ∈ H and g 0 j = τ j v 0 . Then the PMR forms (41) and (40) are equivalent to [START_REF] Smoluchowski | Über die Wechselwirkung von Kugeln, die sich in einer zähen Flüssigkeit bewegen[END_REF] and [START_REF] Toselli | Domain Decomposition Methods: Algorithms and Theory[END_REF].

Proof. We proceed as in Section 3.3 and consider the transmission condition of [START_REF] Richardson | The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam[END_REF], that is

(46) τ j v k j = g j - J =j+1 τ j v k-1 - j-1 =1 τ j v k-1 .
We define g k j := τ j v k j and recall that v k-1 j can be written as

v k-1 j = G j (τ j v k-1 j ) = G j ( g k-1 j
). Hence, (46) becomes

g k j = g j - J =j+1 τ j G j ( g k-1 ) - j-1 =1 τ j G j ( g k-1
), and α 1,3 = α 2,1 = α 3,2 = 1, which corresponds to the classical PSM, the second case is α j,k = 1 2 for any j, k, which leads to a weighted PSM, and the third case α j,k = -1 for any j, k, which corresponds to the PMR. Computing explicitly the iteration matrix G is these cases, we get

G 1 =         0 0 0 0 1 5 0 0 0 0 0 2 5 0 0 4 5 0 0 0 0 0 3 5 0 0 0 0 0 0 0 2 3 0 0 0 0 0 1 3 0 0         , G 2 =         0 0 1 6 0 1 10 0 0 0 1 3 0 1 5 0 0 2 5 0 0 3 10 0 0 3 10 0 0 2 5 0 0 1 5 0 1 3 0 0 0 1 10 0 1 6 0 0         , G 3 = -         0 0 1 3 0 1 5 0 0 0 2 3 0 2 5 0 0 4 5 0 0 3 5 0 0 3 5 0 0 4 5 0 0 2 5 0 2 3 0 0 0 1 5 0 1 3 0 0         ,
where G 1 , G 2 , and G 3 correspond to the three cases. Their spectra are

σ(G 1 ) = i √ 3 4 1/3 -4 1/3 2 25 1/3 , - 4 1/3 + i √ 3 4 1/3 2 25 1/3 , 4 1/3 25 1/3 , 0 , σ(G 2 ) = - 3 + √ 39
ρ(G 1 ) = 4 1/3 25 1/3 ≈ 0.54, ρ(G 2 ) = √ 219 + 3 30 ≈ 0.59, ρ(G 3 ) = √ 219 + 3 15 ≈ 1.18.
We see again that the iterations generated by G 1 and G 2 corresponding to parallel Schwarz methods are different from the iterations generated by G 3 corresponding to the PMR, and we have therefore shown that there is no relation in general between PMR and PSM for J > 2.

4.3.

Convergence analysis for J = 2 objects. Next, we prove for the PSM the same convergence result obtained for the AltMR in Theorem 4. The PMR for the errors e k 1 and e k 2 is ∆e Proof. As in Corollary 1, the functions w j solving (32) satisfy by the maximum principle that w j ≥ 0 and |e k j | ≤ w j max ∂Oj |e k j |. Now, since e k 1 and e k 2 are harmonic functions and solve [START_REF] Kress | Linear Integral Equations[END_REF], the maximum principle allows us to write max

k 1 = 0 in Ω \ ∂O 1 , e k 1 = 0 on ∂Ω, e k 1 = -e k-1 2 on ∂O 1 , ∆e k 2 = 0 in Ω \ ∂O 2 , e k 2 = 0 on ∂Ω, e k 2 = -e k-
Ω |e k 1 | = max ∂O1 |e k 1 | = max ∂O1 |e k-1 2 | ≤ max ∂O1 w 2 max ∂O2 |e k-1 2 | = max ∂O1 w 2 max Ω |e k-1 2 |, max Ω |e k-1 2 | = max ∂O2 |e k-1 2 | = max ∂O2 |e k-2 1 | ≤ max ∂O2 w 1 max ∂O1 |e k-2 1 | = max ∂O2 w 1 max Ω |e k-2 1 |.
By combining these two estimates we get max

Ω |e k 1 | ≤ max ∂O1 w 2 max Ω |e k-1 2 | ≤ max ∂O1 w 2 max ∂O2 w 1 max Ω |e k-2 1 | = ρ AltMR max Ω |e k-2 1 |.
Using this inequality recursively leads to max Ω |e 2k 1 | ≤ ρ k AltMR max Ω |e 0 1 |, which is the claim for j = 1. The same arguments allow us to prove this result also for j = 2.

Notice that Theorem 4 and Theorem 8 are proved using maximum principle arguments that are classical techniques for proving convergence of classical Schwarz methods; see, e.g., [START_REF] Lions | On the Schwarz alternating method. II. Stochastic interpretation and other properties[END_REF][START_REF] Ciaramella | Analysis of the parallel Schwarz method for growing chains of fixed-sized subdomains: Part II[END_REF]. Notice also that the estimated contraction factors for AltMR and PMR are one the square of the other. This resembles classical results for Schwarz methods; see, e.g., [START_REF] Gander | Optimized Schwarz methods[END_REF][START_REF] Gander | Schwarz methods over the course of time[END_REF] and references therein.

Methods of reflections with relaxation

In the following subsections, we present further forms of MR, and prove that they are equivalent to classical stationary iterative methods with relaxation. In particular, in Section 5.1, we present an averaged version of the parallel method introduced in [START_REF] Laurent | On the method of reflections[END_REF] and show its equivalence with the damped block Jacobi iterations. In Section 5.2, we introduce a new MR that corresponds to a block variant of the well-known method of successive over relaxation (SOR).

5.1.

The averaged parallel method of reflection. Golusin already said in [START_REF] Golusin | Auflösung des dreidimensionalen Dirichleteschen Problems für die Laplacesche Gleichung und Gebiete, die durch endlich viele Sphären ohne gemeinsame Punkte begrenzt sind[END_REF] that the parallel method of reflections always converges for J = 2 objects, but not for J > 2, "but by modifying [START_REF] Gander | A new coarse grid correction for RAS/AS[END_REF] (or ( 42)-( 43)) appropriately one could increase the area of applicability of the preceding result". This fact has been mentioned in several publications; see, e.g., [START_REF] Golusin | Auflösung des dreidimensionalen Dirichleteschen Problems für die Laplacesche Gleichung und Gebiete, die durch endlich viele Sphären ohne gemeinsame Punkte begrenzt sind[END_REF][START_REF] Ichiki | Many-body effects and matrix inversion in low-Reynolds-number hydrodynamics[END_REF]. In order to improve the convergence behavior of the PMR, Laurent et al. have recently proposed in [START_REF] Laurent | On the method of reflections[END_REF] a modified version that is obtained (as mentioned by the authors) by averaging the different components d k j . To do so, the problem ( 36) is modified by adding a weight, ∆d k j = 0 in Ω \ ∂O j , d k j = 0 on ∂Ω,

d k j = 1 - 1 J d k-1 j - 1 J J =1, =j d k-1 on ∂O j , (50)
for j = 1, . . . , J, with the initialization problems

∆d 1 j = 0 in Ω \ ∂O j , d 1 j = 0 on ∂Ω, d 1 j = g j -u 0 on ∂O j , (51)
for j = 1, . . . , J and u 0 ∈ H. The approximate solution u k can be obtained from u k-1 by ( 52)

u k = u k-1 + 1 J J j=1 d k j ,
assuming that u 0 = 0. This new formulation of the method, that we call averaged parallel method of reflections (APMR), is proved to be always convergent in [START_REF] Laurent | On the method of reflections[END_REF]. Now, we want to formulate this new version in terms of v k j . To do so, starting with the last boundary condition in (50) and using (52), we obtain on ∂O j that

d k j = 1 - 1 J d k-1 j - 1 J J =1, =j d k-1 = d k-1 j - 1 J J =1 d k-1 = d k-1 j -(u k-1 -u k-2 ),
which implies that

d k j + u k-1 = d k-1 j + u k-2 = • • • = d 1 j + u 0 = g j ,
where we used that d 1 j = g j -u 0 on ∂O j . Hence, on ∂O j it holds that (53) d k j = g j -u k-1 . Now, we need to define the variable v k j in a slightly different way than [START_REF] Dhont | An introduction to dynamics of colloids[END_REF], namely

(54) v k j := v k-1 j + 1 J d k j with v 0 j ∈ H j ,
where J j=1 v 0 j = u 0 . Then, we can use (52) to obtain, like for ( 13), the relation

(55) u k = u k-1 + 1 J J j=1 d k j = u 0 + 1 J k n=1 J j=1 d n j = u 0 + k n=1 J j=1 v n j -v n-1 j = J j=1 v k j .
Now, using ( 53), (54), and (55), we get

v k j = v k-1 j + 1 J d k j = v k-1 j + 1 J g j - 1 J u k-1 = v k-1 j + 1 J g j - 1 J J =1 v k-1 = 1 - 1 J v k-1 j + 1 J g j - 1 J J =1, =j v k-1 , which holds on ∂O j . Recalling that v k j ∈ H j , it solves the problem ∆v k j = 0 in Ω \ ∂O j , v k j = 0 on ∂Ω, v k j = 1 - 1 J v k-1 j + 1 J g j - 1 J J =1, =j v k-1 on ∂O j . (56)
The equivalence between (56) and ( 50)-( 51) is proved in the following theorem. The block Jacobi method [START_REF] Smoluchowski | Über die Wechselwirkung von Kugeln, die sich in einer zähen Flüssigkeit bewegen[END_REF] we have presented earlier might not converge if the off-diagonal blocks are too heavy, i.e. the objects are too close (in agreement with the results of Golusin [START_REF] Golusin | Auflösung des dreidimensionalen Dirichleteschen Problems für die Laplacesche Gleichung und Gebiete, die durch endlich viele Sphären ohne gemeinsame Punkte begrenzt sind[END_REF] and Traytak [START_REF] Traytak | Convergence of a reflection method for diffusion-controlled reactions on static sinks[END_REF]); see also [START_REF] Ciaramella | Review of the methods of reflections[END_REF]Section 3] for detailed one-dimensional examples. However, one can consider a relaxation, (57)

g k = (1 -ω) g k-1 + ωD -1 -(L + U ) g k-1 + g ,
where the parameter ω has to be chosen in a proper way. This is the relaxed block Jacobi method, and we have the following result.

Theorem 10. We have the following equivalences:

• The relaxed block Jacobi method (57) is equivalent to the damped block Jacobi method, that is

(58) g k = g k-1 + ωD -1 g -A g k-1 .
• If one considers a sequence of relaxation parameters {ω k } k , then the damped block Jacobi method is equivalent to Richardson's method [START_REF] Richardson | The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam[END_REF], namely

g k = g k-1 + ω k g -A g k-1 .
• Assume that v 0 ∈ H and g 0 j = τ j v 0 and ω = 1 J . Then the APMR (56) is equivalent to (58).

Proof. To show the first statement, we compute

g k = (1 -ω) g k-1 + ωD -1 -(L + U ) g k-1 + g = g k-1 -ω g k-1 + ωD -1 g -ωD -1 (L + U ) g k-1 = g k-1 -ω g k-1 + ωD -1 g -ωD -1 (A -D) g k-1 = g k-1 + ωD -1 g -A g k-1 ,
which is the damped block Jacobi iteration (58). The second statement follows easily by recalling that D = I and considering the parameter ω depending on the iteration k. Setting ω = 1 J and recalling the operators D, L, and U given in Section 3.3, we can write (57) in the pointwise form

g k j = 1 - 1 J g k-1 j + 1 J g j - J =1, =j τ j G ( g k-1
) .

Then similar arguments as in Theorem 2 imply that this is equivalent to the APMR in (56).

We have seen that the APMR can be regarded as a relaxed (or damped) block Jacobi method or as Richardson's method with the specific choice ω k = ω = 1 J . This choice, motivated in [START_REF] Laurent | On the method of reflections[END_REF] by a Hilbert projection analysis, is not guaranteed to be optimal. In what follows we derive under suitable conditions, an explicit formula for the optimal parameter ω in discrete settings. Denoting by G J (ω) the damped-block Jacobi iteration matrix obtained, for example, by a boundaryelement discretization of ( 18) and (57), the optimal parameter ω is then given by ω = arg min

ω∈[0,1] ρ(G J (ω)),
where ρ(G J (ω)) is the spectral radius of G J (ω). We have the following results.

Theorem 11. Let λ k (G J (1)), for k = 1, 2, . . . , be the eigenvalues of the block Jacobi iteration matrix. Assume that λ k (G J (1)) ∈ R for any k and that ρ(G J (1)) < 1, that is the block Jacobi method converges. Then

ω = 2 2 -(λ max + λ min ) ,
where λ max and λ min are the maximum and minimum eigenvalues of G J (1).

Proof. Recalling (57), we observe that the damped block Jacobi iteration matrix has the form G J (ω) = (1 -ω)I + ωG J (1). Hence its eigenvalues are

λ k (ω) = 1 -ω(1 -λ k (G J (1)))
for k = 1, 2, . . . . Now, the optimal parameter ω solves the problem min

ω∈[0,1] max k |1 -ω(1 -λ k (G J (1)))|.
Since ρ(G J (1)) < 1, then λ k (G J (1))) ∈ (0, 1). Therefore, the spectral radius ρ(G J (ω)) has the form

ρ(G J (ω)) = max |1 -ω(1 -λ max )| , |1 -ω(1 -λ min )| ,
and the optimal parameter ω is the point where the two straight functions |1 -

ω(1 -λ max )| and |1 -ω(1 -λ min )| intersect, that is ω = 2 2-(λmax+λmin)
, which is our claim.

Theorem 11 provides the optimal parameter under the assumption that the block Jacobi method converges. If block Jacobi does however not converge, can ω be used to make it convergent? In other words, if an eigenvalue of G J (1) is such that |λ k (G J (1))| ≥ 1, is there a choice of ω capable to correct this behavior? To study this unfavorable case, we have to distinguish different cases:

• if λ k (G J (1)) = 1, then λ k (ω) = 1 -ω(1 -λ k (G J (1))) = 1 for any ω. Hence, there exists no ω ∈ R such that |λ k (ω)| < 1. • if λ k (G J (1)) = -1, then λ k (ω) = 1 -2ω, which implies that |λ k (ω)| < 1 if and only if 0 < ω < 1. • if λ k (G J (1)) < -1, then |λ k (ω)| = |1 -ω(1 -λ k (G J (1)))| < 1 ⇔ 0 < ω < 2 1 -λ k (G J (1)) . • if λ k (G J (1)) > 1, then |λ k (ω)| = |1 -ω(1 -λ k (G J (1)))| < 1 ⇔ 2 1 -λ k (G J (1)) < ω < 0.
We can summarize these facts in the following result.

Theorem 12. Let G J (1) and G J (ω) be the block Jacobi and damped block Jacobi iteration matrices and assume that |λ k (G J (1))| ≥ 1 for a given k. We have that

• if λ k (G J (1)) = -1, then |λ k (ω)| < 1 for any ω ∈ (0, 1). • if λ k (G J (1)) = 1, then |λ k (ω)| = 1 for any ω ∈ R. • if λ k (G J (1)) < 0, then |λ k (ω)| < 1 if and only if 0 < ω < 2 1-λ k (GJ(1)) . • if λ k (G J (1)) > 1, then |λ k (ω)| < 1 if and only if 2 1-λ k (GJ(1)) < ω < 0.
Theorem 12 shows that if block Jacobi does not converge, it is not always possible to find a parameter ω ∈ (0, 1] that makes the method convergent. If the eigenvalues λ k (G J (1)) are such that λ k (G J (1)) < 1, i.e. there can be arbitrarily large negative eigenvalues, a case that we have observed numerically, then a similar proof as for Theorem 11 allows us to obtain the following result.

Theorem 13. Let G J (1) and G J (ω) be the block Jacobi and damped-block Jacobi iteration matrices and assume that the eigenvalues λ k (G J (1)) are real and such that λ k (G J (1)) < 1, then the optimal parameter is ω = 2 2-(λmax+λmin) , where λ max and λ min are the maximum and minimum eigenvalues of G J (1). 5.2. Successive-Over-Relaxation method (SOR). In the previous sections, we have seen that the different forms of the method of reflections are equivalent to block Gauss-Seidel, block Jacobi, damped block Jacobi, and Richardson's method. This suggests that a possible improvement of the method of reflections would be to rewrite the AltMR as an SOR method, which converges in general much faster for a well chosen relaxation parameter [START_REF] Young | Iterative Methods for Solving Partial Difference Equations of Elliptic Type[END_REF]. SOR for the system ( 18) is (59)

g k = (D + ωL) -1 ωg + -ωU + (1 -ω)D g k-1 .
Recalling the form of D, L, and U given in Section 3.3, (59) written component-wise is

(60) g k j = (1 -ω) g k-1 j + ω g j - j-1 =1 τ j G ( g k ) - J =j+1 τ j G ( g k-1
) .

Defining v k j := G j ( g k j ), this iteration can be rewritten in the volume form,

∆v k j = 0 in Ω \ ∂O j , v k j = 0 on ∂Ω, v k j = 1 -ω v k-1 j + ω g j - j-1 =1 v k - J =j+1 v k-1 on ∂O j , (61)
with u k = J j=1 v k j . This is a new formulation of the method of reflections that can be also written in terms of the (more usual) d k j , that is

∆d k j = 0 in Ω \ ∂O j , d k j = 0 on ∂Ω, d k j = 1 -ω d k-1 j -ω j-1 =1 d k + J =j+1 d k-1 on ∂O j , (62)
with

∆d 1 j = 0 in Ω \ ∂O j , d 1 j = 0 on ∂Ω, d 1 j = g j -u 0 -ω j-1 =1 d 1 on ∂O j , (63)
and

u k = u k-1 + ω J j=1 d k j = u 0 + ω k n=1 J j=1 d n j .
Notice that (62)-( 63) and (60) are equivalent, because (60) is the equivalent substructured form of (61) and the equivalence between (61) and ( 62)-( 63) is proved in the following theorem. Theorem 14. Consider the sequences {u k } k ⊂ H with u 0 ∈ H, {d k j } k with d k j ∈ H j , and {v k j } k with v k j ∈ H j for j = 1, . . . , J, and assume that v k j = v k-1 j + ωd k j for k ≥ 1 and ω > 0 and v 0 j ∈ H j for j = 1, . . . , J. Then the following statements are equivalent:

(a) the d k j , j = 1, . . . , J, solve (62)-(63) and u k = u k-1 +ω J j=1 d n j = ω k n=1 J j=1 d n j . (b) the v k j , j = 1, . . . , J, solve (61) and u k = J j=1 v k j . Proof. (a) ⇒ (b): first, using v k j = v k-1 j + ωd k j , we write that (64) u k = u k-1 + ω J j=1 d k j = u 0 + ω k n=1 J j=1 d n j = u 0 + k n=1 J j=1 v n j -v n-1 j = J j=1 v k j .
Notice that v k j ∈ H j , hence it is sufficient to prove that the boundary condition on ∂O j in (61) holds. Using the transmission condition of (62) and

u k-1 = u k-2 + ω J =1 d k-1 , we obtain d k j = d k-1 j -ω j-1 =1 d k -ω J =j d k-1 = d k-1 j -ω j-1 =1 d k -ω J =1 d k-1 + ω j-1 =1 d k-1 = d k-1 j -ω j-1 =1 d k -u k-1 -u k-2 + ω j-1 =1 d k-1 ,
which implies that

d k j + u k-1 + ω j-1 =1 d k = d k-1 j + u k-2 + ω j-1 =1 d k-1 = • • • = d 1 j + u 0 + ω j-1 =1 d 1 = g j ,
where we used the transmission condition of (63). Hence, we have obtained that

d k j = g j -u k-1 -ω j-1 =1 d k .
Replacing this into the transmission condition of (62) and using v k j = v k-1 j + ωd k j and (64), we obtain

1 ω v k j -v k-1 j = d k j = d k-1 j -ω j-1 =1 d k -ω J =j d k-1 = g j -u k-2 -ω J =1 d k-1 -ω j-1 =1 d k = g j -u k-1 -ω j-1 =1 d k = g j -u k-1 - j-1 =1 v k -v k-1 = g j - j-1 =1 v k - J =j v k-1 ,
which implies that

v k j = v k-1 j +ωg j -ω j-1 =1 v k -ω J =j v k-1 = (1-ω)v k-1 j +ω g j - j-1 =1 v k - J =j+1 v k-1 ,
which is our claim. (b) ⇒ (a): using

d k j = 1 ω v k j -v k-1 j
, we obtain

u k = J j=1 v k j = u 0 + J j=1 v k j -v 0 j = u 0 + k n=1 J j=1 v n j -v n-1 j = u 0 + ω k n=1 J j=1 d n j , which implies that u k = u k-1 + ω J j=1 d k j . Now, since d k j = 1 ω v k j -v k-1 j
we have that d k j ∈ H j . It remains to show that the transmission conditions of (62)-(63) hold. For k = 1 we have

d 1 j = 1 ω (v 1 j -v 0 j ) = g j - j-1 =1 v 1 - J =j v 0 = g j -ω j-1 =1 d 1 - j-1 =1 v 0 - J =j v 0 = g j -ω j-1 =1 d 1 -u 0 .
In the case k > 1, we write the transmission condition of (61) on ∂O j for k and k -1:

v k j = 1 -ω v k-1 j + ω g j - j-1 =1 v k - J =j+1 v k-1 , v k-1 j = 1 -ω v k-2 j + ω g j - j-1 =1 v k-1 - J =j+1 v k-2 .
Subtracting term by term these equations, dividing by ω, and recalling that

d k j = 1 ω v k j -v k-1 j
we obtain the transmission condition of (62).

A good choice of the parameter ω is not trivial. In general one would consider the optimal parameter ω , that is the parameter that minimizes the contraction rate of the SOR method. From a discrete point of view, namely when one is interested in solving a linear system A g = g using a stationary iterative method of the form g k+1 = M -1 N g k + M -1 g, there are several results in the literature for point relaxation (in contrast to block relaxation here). A famous result proved by Kahan [START_REF] Kahan | Gauss-Seidel Methods of Solving Large Systems of Linear Equations[END_REF] says that for the SOR method one has ρ(M -1 N ) ≥ |1 -ω|, which implies that a necessary condition for convergence is ω ∈ (0, 2). In the case the matrix A has the so called Property A, David Young proved in [START_REF] Young | Iterative Methods for Solving Partial Difference Equations of Elliptic Type[END_REF] that the optimal parameter, that is the parameter that minimizes the spectral radius ρ(G SOR (ω)) with respect to ω , where G SOR (ω) is the SOR-iteration matrix, is given by ω

= 2 1+ √ 1-ρ(GJ) 2
, where G J is the iteration matrix of the corresponding point-Jacobi method applied to the same linear system. A discrete form of our problem [START_REF] Gander | SHEM: An optimal coarse space for RAS and its multiscale approximation[END_REF] for J = 2 is characterized by a matrix A of the form

A = I 1 G 2 G 1 I 2 ,
where I 1 and I 2 are identity matrices. This shows that A is exactly in the form required to have the Property A. Therefore, we can apply the theory developed by David Young to get an explicit formula for the optimal parameter. This observation is not in general true for J > 2. Moreover, to compute ω one would need the spectral radius of the block Jacobi method (the PMR). This is not an easy task, even for J = 2. However, we can estimate ρ(G J ) at a continuous level using maximum principle arguments, as we show in Section 3.5. Numerical experiments in Section 7 show that this leads to a very good estimate of ω for J = 2. In the case that J > 2 the optimal parameter depends strongly on the geometry and we observed that when the distance between the objects is sufficiently large the optimal choice is ω ≈ 1.

Scalability analysis and coarse correction

An iterative method is said to be scalable, if its contraction factor does not deteriorate when the number of unknowns grows; see, e.g., [START_REF] Toselli | Domain Decomposition Methods: Algorithms and Theory[END_REF] and references therein. In the field of domain-decomposition methods for the solution of partial differential equations, we distinguish different scalability problems depending on the way the number of unknowns increases: The first property is widely studied in the literature, see, e.g., [START_REF] Toselli | Domain Decomposition Methods: Algorithms and Theory[END_REF] and references therein. The case (b) has been only recently observed in [START_REF] Cancès | Domain decomposition for implicit solvation models[END_REF] and theoretically investigated in [START_REF] Ciaramella | Analysis of the parallel Schwarz method for growing chains of fixed-sized subdomains: Part I[END_REF][START_REF] Ciaramella | Analysis of the parallel Schwarz method for growing chains of fixed-sized subdomains: Part II[END_REF]9] for classical Schwarz methods and for classical one-level Studying methods of reflections, we observed numerically that scalability in the case (a) holds, that is the convergence of the method is not affected by the number of discretization points of the boundaries of the objects ∂O j , for j = 1, . . . , J. On the other hand, scalability in the sense (b) does not hold, as we will show below. Moreover, if the domain Ω is fixed and one increases the number of objects, then obviously the distance between them reduces, and we have seen in Section 9 for J = 2 that the contraction factor deteriorates in this case. The same behavior has been observed, e.g., in [START_REF] Laurent | On the method of reflections[END_REF][START_REF] Ciaramella | Review of the methods of reflections[END_REF]. This implies the need of a coarse-correction to address these unfavorable behaviors.

x = 0 1 2 3 J -2 J -1 J x = J δ δ δ δ δ δ δ
We consider a one-dimensional problem characterized by J equidistant objects (subintervals); see Figure 4. This is a finite-dimensional problem because the unknowns are the Dirichlet data g ∈ R 2J on the 2J extrema of the objects. The problem can then be written in a form A g = g, where g ∈ R 2J and A ∈ R 2J×2J . Using the equivalence result proved in Theorem 2, one can easily construct the block Gauss-Seidel (AltMR) iteration matrix, that is G GS = -(D + L) -1 U and numerically compute its spectral radius for increasing number of objects J. The result is given in Figure 5, which clearly shows that the AltMR does not scale because ρ(G GS ) deteriorates for growing J. The heuristic reason for this behavior is that, once the approximate solution is corrected on one object, say the jth one, then the procedure continues to correct it sequentially on all the other objects. Every correction induces an additional error in the approximate solution on the jth object. Therefore, when J increases, more and more objects have to be corrected, which adds more and more error induced by the alternating correction procedure on the jth object. Similar arguments apply also for the other methods of reflections. We thus need a coarse correction to obtain a scalable method. A spectral analysis of the one-dimensional iteration matrix G GS reveals that many of its eigenvalues are clustered around 1, see Figure 5 (right), and that the corresponding eigenvectors have a special structure: they globally oscillate with (local) maxima and minima attained on the objects, and many of them have constant value on the objects; see Figure 6. This suggests that a coarse-correction has mainly to deal with errors on the boundaries of the objects. Therefore, to design our coarse-space we denote by ψ j,n for n = 0, 1, . . . the eigenfunction of the Laplace-Beltrami operator on the jth object ∂O j , for j = 1, . . . , J, and introduce functions ϕ j,n as solution to the Dirichlet problem

∆ϕ j,n = 0 in Ω \ ∪ J =1 ∂O , ϕ j,n = ψ j,n on ∂O j , ϕ j,n = 0 on ∂Ω ∪ ∪ J =1, =j ∂O , (65)
for j = 1, . . . , J and n = 0, 1, . . . . Notice that ϕ j,n ∈ H for any j and n, and the functions ψ j,n are classical Fourier-basis functions in two dimensions and spherical harmonics in three dimensions. Our coarse-space is then defined as

V N c := span j=1,...,J, n=0,1,...,N {ϕ j,n },
where N is the number of eigenfunctions considered. Notice that the dimension of V N c is proportional to N . It is clear that the coarse space V N c has mainly information condensed on the boundaries of the objects, similar to the Spectral Harmonically Enriched Multiscale coarse space SHEM in domain decomposition [START_REF] Gander | Analysis of a new harmonically enriched multiscale coarse space for domain decomposition methods[END_REF][START_REF] Gander | SHEM: An optimal coarse space for RAS and its multiscale approximation[END_REF], which contains mainly information on the interfaces between subdomains, see also [START_REF] Gander | A new coarse grid correction for RAS/AS[END_REF][START_REF] Gander | Discontinuous coarse spaces for DDmethods with discontinuous iterates[END_REF]. It is important to remark that the construction of each function ϕ j,n would require the solution of problem (65), which requires the same computational effort of the original problem (3). However, this is not needed because the substructured formulation introduced in Section 3.3 allows us to work directly on the boundary of the objects: our coarse space contains thus spectral approximations of the substructured problem. For this reason, we introduce the restriction operator and its corresponding adjoint (prolongation or extension) operator as

R :=      F N 1 F N 2 . . . F N J      and R * =      (F N 1 ) * (F N 2 ) * . . . (F N J ) *     
, where F N j : L 2 (∂O j ) → R N is given by

F K j (w) = c 1,j . . . c N,j =: c j ,
with c k,j = ∂Oj wf k,j for any w ∈ L 2 (∂O j ) and j = 1, . . . , J, and (F

N j ) * : R K → L 2 (∂O j ) is (F N j ) * (c j ) = K k=1 c k,j ψ k,j ,
for c j ∈ R K . These operators allow us to restrict the operator A introduced in ( 18) and ( 19) on the coarse space V N c as A c = RAR * ∈ R JN ×JN . We are now ready to state our two-level method of reflection: given an approximation g k to g at the iteration k, a coarse-corrected method of reflections step is defined as

g k+1/2 := G MR g k + M -1 MR g (one-level MR step), r k+1/2 := g -A g k+1/2
(compute the residual),

g k+1/2 c := A -1 c R r k+1/2
(compute the correction),

g k+1 := g k+1/2 + R * g k+1/2 c (correct g k+1/2 ),
where G MR is the one-level method of reflection operator and M MR the corresponding preconditioning matrix, which can be, e.g., G GS (with M MR = (D + L)) for the AltMR or G J (with M MR = D) for the PMR. A direct calculation reveals that the two-level method of reflections iteration operator is

G MR-c.c. = I -R * A -1 c RA G MR .
A coarse corrected method of reflection is given in Algorithm 5.

Algorithm 5 Method of Reflections with Coarse Correction

Input: K (maximum number of iterations), tol (tolerance), g 0 (initial guess).

Input: fMR(v, w) := GMR v + M -1 MR w (function that performs one step of a method of reflection MR), Input: g (data of the problem).

1: Compute r 0 = g -A g 0 and set k = 0. 2: while k < K and r k > tol do 3:

Set k = k + 1.

4:

g k = fMR( g k-1 , g).

5:

r k = g -A g k . 6: gc = A -1 c R r k . 7:
g k+1 = g k + R gc. 8: end while Notice that like in domain decomposition, Algorithm 5 requires the solution of the problem projected onto the coarse space. If the dimension of the coarse space is proportional to J, then this solution requires in principle O(J 3 ) flops. Therefore, for very large J the cost of the coarse correction becomes comparable to the cost of solving the full problem by a direct method. In this case, our two-level framework would need to be extended to a multi-level setting, which is however beyond the scope of this paper.

Numerical experiments

We now present some numerical tests to illustrate the results we obtained, and to explore cases not covered by the theory we developed. We start with the case of two objects, and then investigate the case of three objects. We then study the scalability of the method with respect to the number of objects and finally consider the method as a preconditioner for the GMRES algorithm. We used a publicly available package 1 of Matlab functions to solve the integral equations of the problem at collocation points by the Nyström method. 7.1. Two objects. We first consider the case of two objects of radius r = 0.2, in the unit disc, for which we have theoretical results. The theoretical value of the optimal ω for the SOR variant of the method of reflections is given by ( 66)

ω (r, δ) = 2 1 + 1 -ρ PMR (r, δ) 2 .
With this choice of ω, the theoretical convergence factor is

(67) ρ SOR (r, δ) ω=ω = ω (r, δ) -1
In a first experiment, we consider two cases where the distance between the objects is either relatively large, that is δ = 0.5, or relatively small that is δ = 0.25, see Figure 7. We consider both SOR and damped PMR. In the case of SOR, ω is given by Equation (66), where, in our 2 object case, ρ PMR (r, δ) can be computed numerically and estimated theoretically. Alternatively the latter quantity can be estimated using formulas provided in Theorem 8 and Corollary 1. Results are presented in Figure 8. We see that the agreement between the theoretical and numerical values of ω is good when the objects are far from each other, and less accurate when the objects are close to each other. We also see that the relaxation parameter ω is only close to 1 in the first case, which means that the SOR variant provides a real improvement with respect to the standard AltMR when the objects are relatively close to each other. In the case of PMR, we observe that the damping does not improve the contraction factor, see Figure 9.

The theoretical estimate of the contraction factor is rather accurate and does not seem to depend on the distance between the objects. Note also that the eigenvalues of the PMR iteration matrix form a symmetric set with respect to zero. Indeed, one has As a consequence, if λ is an eigenvalue, then also -λ is an eigenvalue. Because of Theorem 11, we find that in this case ω = 1, as observed in Figure 9. Finally, we compare theoretical and numerical values of the contraction factors both for AltMR and PMR when the distance δ between the objects varies; the results are shown in Figure 10, and we see a good agreement between theory and numerics. 7.2. Three objects. We now consider the case of three objects of radius r = 0.2 in the unit disc, for which we do not have theoretical results. For the sake of simplicity, we restrict ourselves to SOR, and repeat the experiments done to obtain Figure 8. The results are presented in Figure 11. We observe again that SOR provides a real improvement compared to the standard AltMR, in the sense that the observed optimal ω differs from 1. We however also see that the estimates obtained with Equation (66) are not accurate now. 7.3. Scalability with respect to the number of objects. In this test, we study the scalability of the AltMR and PMR when the number of objects varies, and the effect of the coarse correction on both methods. We consider for a given L ∈ N + a square of size L × L with rounded corners containing J = L 2 objects of radius r = 0.1, see Figure 12. We then evaluate numerically the contraction factor for both methods, with and without coarse correction, see Figure 13. We consider here the coarse space V 0 c = span j=1,...,J {ϕ j,0 } of dimension J, where we use only the first (constant) Fourier mode ψ j,0 for each object, see (65). We observe that the coarse correction significantly reduces the contraction factor, by approximately one order of magnitude for the AltMR and even more for the PMR, where the coarse correction leads to a convergence factor below 1, so that the PMR with coarse correction remains convergent even in cases where the standard PMR diverges. We next repeat these two tests with the method of reflection variants with relaxation, namely, with the SOR variant and the APMR, see Figure 14. We see that the coarse correction significantly improves the performance of the SOR variant, but not the APMR. It even seems that the convergence rate remains unchanged in the case of APMR. This is due to the relaxation parameter ω = 1 J that decays as J grows.

G J = -D -1 (L + U ) = - 0 G G 0 ,
Next, we solve the problem using GMRES and employing the methods of reflections as preconditioners, both with and without coarse correction. For reference, we compare the performance of these methods with the coarse-corrected methods of reflections used as stationary iterative methods; the required iteration numbers are summarized in Table 1. We observe that with GMRES now the PMR preconditioner leads to a convergent method, even though as a stationary iterative method it was not convergent. This is very similar to the additive Schwarz method which is Table 2. Same experiment setting as in Table 1, but for a larger object radius r = 0.3.

also not converging as a stationary iterative method [START_REF] Efstathiou | Why restricted additive Schwarz converges faster than additive Schwarz[END_REF][START_REF] Gander | Schwarz methods over the course of time[END_REF]. We also see that with our coarse correction one obtains scalability with respect to the number of objects. This scalability also holds for the stationary iterative variants. We next repeat the experiment using the same geometry but increasing the radius of the objects to r = 0.3, which implies that the distance between the objects is smaller, and thus a slower convergence rate of the methods of reflections. The number of iterations are summarized in Table 2. We see that indeed iteration numbers are now larger, but the methods of reflections with our coarse correction are still scalable.

9 :

 9 Compute d k j solving problem[START_REF] Ciaramella | Review of the methods of reflections[END_REF].

Theorem 3 .

 3 Consider J = 2 and assume that d 0 2 = d 0 2 . Then for any k ≥ 1 we have that d k 2 = d k 2 on Ω 1 and d k 1 = -d k 1 on Ω 2 , where d k j and d k j solve (27) and (10).

Figure 2 .

 2 Figure 2. A unit disc Ω with two objects O 1 and O 2 , which are two discs having the same radius denoted by r. The distance between the centers of O 1 and O 2 and the center of Ω is denoted by δ.

Figure 3 .

 3 Figure 3. Möbius transformation of the unit disc Ω (dashed line) with a hole O 1 (solid line) into an annulus.

9 :

 9 Compute d k j solving problem[START_REF] Lions | On the Schwarz alternating method. II. Stochastic interpretation and other properties[END_REF].

  (a) the d k j solve (36)-(37) with u k = u k-1 + the v k j solve (41) with u k = J j=1 v k j . Proof. The implication (a) ⇒ (b) is proved in the discussion before the theorem, whereas (b) ⇒ (a) follows by similar arguments as in Theorem 1.

Algorithm 4

 4 3 and the decomposition A = M -N with M := D and N := -(L + U ). The block Jacobi method in the standard form is given by (42) D g k+1 = -(L + U ) g k + g, Parallel Method of Reflection (PMR)

Theorem 9 . 1 j + d k j for k ≥ 1

 911 Consider the sequences {u k } k∈N ⊂ H, {d k j } k∈N + with d k j ∈ H j , and {v k j } k∈N with v k j ∈ H j for j = 1, . . . , J, and assume that v k j = v k-and j = 1, . . . , J. Then the following statements are equivalent: (a) the d k j solve (50)-(51) with u k = u k-1 + the v k j solve (56) with u k = J j=1 v k j . Proof. The implication (a) ⇒ (b) is proved in the discussion before the theorem, whereas (b) ⇒ (a) follows by similar arguments as in Theorem 1.

  (a) the number of discretization points increases, but number and size of subdomains is unchanged. (b) the subdomains have fixed size and their number increases making also the size of the overall domain increase.

Figure 4 .Figure 5 .

 45 Figure 4. Geometry of the one-dimensional example used to study scalability of the AltMR. In this figure, δ is the distance between the centers of the J equidistant objects that are represented by the thick segments.

Figure 6 .

 6 Figure 6. Example of an eigenvector of the one-dimensional iteration matrix G GS corresponding to an eigenvalue λ ≈ 1. The right plot shows a zoom into the highly oscillating area of the eigenvector depicted in the left plot.

Figure 7 .

 7 Figure 7. Geometry of the problem under consideration, domain with some normal unit vectors. Left: objects relatively far from each other, right: objects relatively close to each other.

Figure 8 .

 8 Figure 8. Contraction factor of the SOR variant of the method of reflections with respect to ω (stars), theoretical values of ω given by (66) where ρ PMR (r, δ) is obtained either using Theorem 8 and Corollary 1 (dashed line), or computed numerically (dotted line). In the latter case, the optimal contraction factor is computed by Equation (67) (solid line). Left: the objects are far away from each other. Right: the objects are close to each other.

Figure 9 .

 9 Figure 9. Contraction factor of the damped PMR with respect to ω (stars), theoretical estimate of the contraction factor (solid line) given by Theorem 8 and Corollary 1. Left: the objects are far from each other. Right: the objects are close to each other.

  so that det(G J -λI) = det -λI -G -G -λI = det (-λI) det -λI + 1 λ G 2 = det -λ 2 I + G 2 .

Figure 10 .

 10 Figure 10. Contraction factors of AltMR (left) and PMR (right) with respect to δ and theoretical contraction factors (solid lines) given by Theorems 4 and 8 and Corollary 1, compared with numerical contraction factors (stars).

Figure 11 .

 11 Figure 11. Left: geometry of the problem, domain with some normal unit vectors. Right: Contraction factor of the SOR variant with respect to ω (stars), theoretical values of ω given by (66)where ρ PMR (r, δ) is obtained either from Theorem 8 and Corollary 1 (dashed line) or computed numerically (dot line). In the latter case, the optimal contraction factor is computed by Equation (67) (solid line).

Figure 12 .

 12 Figure 12. Geometry of the problem under consideration for L = 3 (left) and L = 4 (right).

Figure 13 .

 13 Figure 13. Contraction factors of the AltMR (left) and PMR (right) with respect to the number of objects. Contraction factor without (crosses) and with (circles) coarse correction.

Figure 14 .

 14 Figure 14. Contraction factors of the SOR variant with ω = 1.3 (left) and APMR (right) with respect to the number of objects. Contraction factor without (crosses) and with (circles) coarse correction.

  Theorem 8. The PMR for the solution to problem (29) converges geometrically, in the sense that max Ω |e 2k j | ≤ ρ 2k PMR max Ω |e 0 j |, for j = 1, 2, where ρ PMR = √ ρ AltMR and ρ AltMR is given in Corollary 1.

	(49)	1	1	on ∂O 2 .

Table 1 .

 1 Number of GMRES iterations necessary to get a residual smaller than 10 -10 . The different columns correspond to different numbers of objects and show the number of iterations performed by GMRES and various preconditioners: AltMR, PMR with coarse correction, and AltMR with coarse correction. In the two last lines, we show the number of iterations required by coarse corrected AltMR and PMR used as stationary iterations.

	# Objects (r = 0.3)	4	9	16 25 36 49 64 81 100 121 144
	GMRES-PMR	10 15 19 22 25 28 30 33	34	37	39
	GMRES-AltMR	7	11 15 19 23 27 32 36	40	44	49
	GMRES-PMR c.c.	5	6	6	6	6	6	7	6	7	7	7
	GMRES-AltMR c.c.	5	6	6	7	7	8	8	8	8	8	8
	AltMR c.c.	18 17 21 19 22 20 22 22	22	23	23
	PMR c.c.	12 15 18 20 22 23 23 23	24	25	26

which is equivalent to ( 42)- [START_REF] Traytak | Convergence of a reflection method for diffusion-controlled reactions on static sinks[END_REF] if g 0 j = g 0 j . The same arguments can be used to show that [START_REF] Lions | On the Schwarz alternating method. II. Stochastic interpretation and other properties[END_REF] is equivalent to (43)- [START_REF] Young | Iterative Methods for Solving Partial Difference Equations of Elliptic Type[END_REF]. [START_REF] Burgers | Hydrodynamics. -On the influence of the concentration of a suspension upon the sedimentation velocity[END_REF].2. Analogies with Schwarz methods. Similar to Section 3.4, we show now that if J = 2, then the PMR is the analog of the classical parallel Schwarz method (PSM), and if J > 2, the PMR and the PSM define different iterates. To do so, we recall [START_REF] Höfer | The method of reflections, homogenization and screening for Poisson and Stokes equations in perforated domains[END_REF] and apply the PSM,

where the coefficients α j, are non-negative and such that α j, = 1. This definition implies that the corresponding differences

Theorem 7. Consider J = 2 and assume that d 0 2 = d 0 2 . Then, for any k ≥ 0 we have that

, where d k j and d k j solve (48) and (41).

Proof. Since J = 2, we have that α 1,2 = α 2,1 = 1. We only prove the first relation

, since the second follows by the same arguments. We proceed by induction. The statement is true for k = 0 by assumption. Now, we assume that d k 1 = (-1) k d k 1 holds and we show that it remains true for k + 1. Recalling the transmission condition of (48) and ( 41), we have

where we used the induction hypothesis. Since (48) and ( 41) admits unique solutions, we have that

in Ω 1 , which is our claim. Now, to show that for J > 2 the analogy proved in Theorem 7 is not in general true, we consider the same example provided in Section 3.4 and use the same notation. In this case, the PSM is

The corresponding iteration in terms of the constants A k j and B k j is given by v k = Gv k-1 , where G := -( L + U ), with G depending on the weights α ,j . Notice that the same relation (with appropriately chosen weights in G) can be obtained for the PMR. We consider three different cases. The first case is α 1,2 = α 2,3 = α 3,1 = 0 Table 3. Same experiment setting as in Table 1, but for a fixed number of objects when the mesh is refined.

We finally repeat this test with a fixed number of J = L 3 = 9 objects of radius r = 0.1, but using various step sizes for the boundary discretization. The results are shown in Table 3. We observe scalability with respect to the mesh size in all cases.

Conclusions

We presented an extensive analysis of the alternating and parallel methods of reflections using domain decomposition techniques. We proved that for two objects, the methods of reflections can be identified with Schwarz domain decomposition methods, which led to new, sharp convergence estimates in this case. We also showed that for more than two objects, the methods of reflections are different from Schwarz domain decomposition methods. We then used substructuring techniques from domain decomposition to rewrite the methods of reflections only iterating on traces, which allowed us to identify these methods with block Gauss-Seidel and block Jacobi methods. Using this insight, we derived new, relaxed variants of the methods of reflections, which converge faster that the classical variants, and also introduced for the first time a coarse correction for the methods of reflections to make them scalable when the number of objects becomes large. The substructured formulation allowed us also naturally to use these methods as preconditioners for Krylov methods. All our results were obtained for the case of the Laplacian, but other operators could be used as well in the relations we found.

University of Konstanz, Germany

Email address: gabriele.ciaramella@uni-konstanz.de University of Geneva, Switzerland Email address: martin.gander@unige.ch Université Paris 13, France Email address: halpern@math.univ-paris13.fr INRIA Paris & Laboratoire Jacques-Louis Lions, Sorbonne Université Email address: julien.salomon@inria.fr