# Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients. Part 1: the constant coefficients case.

1. Introduction. Numerous geophysical phenomena with a strong societal impact involve the coupled ocean-atmosphere system; e.g., those for climate change, tropical cyclones, or sea-level rise predictions. To get a good depiction of the complex air-sea dynamics, it is often necessary to couple atmospheric and oceanic simulation models. However, connecting the two model solutions at the air-sea interface is a difficult problem, which is presently often addressed in a simplified way from a mathematical point of view. Indeed, with the ad-hoc coupling methods currently in use, the fluxes exchanged by the two models are generally not in exact balance [17]. This may be one factor explaining the generally observed strong sensitivity of coupled solutions to the initial conditions or parameter values [23]. This kind of coupling raises a number of challenges in terms of numerical simulation since we are considering two highly turbulent fluids with widely different scales in time and space. It is thus natural to use some specific numerical treatment to match the physics of the two fluids at their interface. It is known that, even if numerical models are much more complicated, a simple one-dimensional diffusion equation is relevant to locally represent the turbulent mixing in the boundary layers encompassing the air-sea interface. The corresponding diffusion coefficients are spatially variable and their values are given by a so-called eddy-viscosity closure [21]. To perform this coupling in a more consistent way than ad-hoc methods, we propose here to adapt a global-in-time domain decomposition based on an optimized Schwarz method. This type of method is thoroughly described in [9] and designed thanks to the pioneering work in [12, 14]. Schwarz-like domain decomposition methods provide flexible and efficient tools for coupling models with non-conforming time and space discretizations [3, 10]. Transmission conditions of Robin type have been proposed in [19] to circumvent the inability of the classical Schwarz method (i.e., with the exchange of Dirichlet data) to solve coupling problems in the case of non-overlapping subdomains. Then, thanks to the free parameters in the Robin conditions, an optimization of the convergence speed has been proposed in [12, 15]: this is the basis of the so called optimized Schwarz methods (OSM). This kind of method, originally introduced for stationary problems, has been extended to the unsteady cases by adapting the waveform relaxation algorithms to provide a global-in-time Schwarz method [14, 16] (sometimes referred to as Schwarz waveform relaxation). This notion of optimization of the convergence speed is critical in the context of ocean-atmosphere coupling as the numerical codes involved are very expensive from a computational point of view. In the present series of two papers, we derive interface conditions leading to an efficient Schwarz coupling algorithm between two unsteady diffusion equations defined on non-overlapping subdomains. The convergence properties of this kind of problem have already been extensively studied in the case of a constant diffusion coefficient having the same value in all subdomains [8]. There are some asymptotic results in the case of coefficients with different constant values in the different subdomains [10] (in the more general case of advection-diffusion-reaction equations). In the present papers, we extend these studies to the general case of diffusion coefficients which vary in each subdomain and whose values are different on both sides of the interface. In this first part, we consider the case of diffusion coefficients that do not vary spatially in each medium. We study a zeroth-order two-sided optimized method by considering two different Robin conditions on both sides of the interface. In the second paper [18], the emphasis is on the impact of the spatial variability of the coefficients on the convergence speed.This first paper is organized as follows. In Section 2, we recall the basics of optimized Schwarz methods in the framework of time evolution problems. Sections 3 and 4 are dedicated to the study of a diffusion problem with discontinuous but piecewise constant coefficients. In Section 3, we analytically determine the solution of an optimization problem to improve the convergence speed of a simplified algorithm with only one Robin condition combined with a Neumann condition. In Section 4, we address the more general case of two-sided optimized Robin-Robin transmission conditions determined through a thorough study of the behavior of the convergence factor. Finally in Section 5, some numerical results are shown to prove the efficacy of the optimized algorithms derived in the previous sections.

2. Model problem and optimized Schwarz methods. Our guiding example is the one-dimensional diffusion equation of a scalar u

(2.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [OMEGA] is a bounded domain defined as [OMEGA] =] - [L.sub.1], [L.sub.2] [, ([L.sub.1], [L.sub.2] [member of] [R.sup.+]) and D(x) > 0, for x [member of] Q. In practical applications, Li denotes the bottom of the ocean (of the order of 5 km in the open ocean), while [L.sub.2] is typically the top of the troposphere (of the order of 15 km). This problem is supplemented by an initial condition

u(x,0) = [u.sub.0] (x), x [member of] [OMEGA],

and boundary conditions

[B.sub.1]u(-[L.sub.1], t)= [g.sub.1], [B.sub.2]u([L.sub.2],t) = [g.sub.2], t [member of] [0,T],

where [B.sub.1] and [B.sub.2] are two partial differential operators. In this paper, we always assume that [u.sub.0] [member of] [H.sup.1] ([OMEGA]), f [member of] [L.sup.2] (0,T; [L.sup.2]([OMEGA])) and that D(x) is bounded in the [L.sup.[infinity]]-norm. Note that in actual applications such assumptions are generally fulfilled. Existence and uniqueness results for this problem can be proved following [10] and are not discussed here.

2.1. Formulation of the global-in-time Schwarz method. In the present study, we consider a case where the diffusion coefficient D(x) has one discontinuity in [OMEGA]. This discontinuity represents the transition between two media with heterogeneous physical properties. In this case we define two subdomains, each subdomain having its own diffusion profile [D.sub.j] (x), (j = 1, 2). This amounts to splitting [OMEGA] into two non-overlapping domains [[OMEGA].sub.1] and [[OMEGA].sub.2]; see Figure 2.1. These subdomains communicate through their common interface at [GAMMA] = {x = 0}. (Note that there can be various reasons for such a splitting: different physics, parallelization and/or different numerical treatment requirements.) We propose to use a nonoverlapping global-in-time Schwarz algorithm to solve the corresponding coupling problem. This method consists in iteratively solving subproblems in [[OMEGA].sub.1] x [0, T] and [[OMEGA].sub.2] x [0, T] using the values computed at the previous iteration in the other subdomain as an interface condition at x = 0. The operator L introduced in (2.1) is split into two operators [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] restricted to [[OMEGA].sub.j] (j = 1, 2). Introducing the operators [F.sub.1], [F.sub.2], [G.sub.1], and [G.sub.2] to define the interface conditions, the algorithm reads

(2.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where k = 1, 2,... is the iteration number, and where the initial guess [u.sup.0.sub.2] (0,t) is given. Algorithm (2.2) corresponds to the so-called "multiplicative" form of the Schwarz method. If we replace the interface condition [G.sub.2][u.sup.k.sub.2] = [G.sub.1][u.sup.k.sub.1] on [[OMEGA].sub.2] by [G.sub.2][u.sup.k.sub.2] = [G.sub.1][u.sup.k-1.sub.2], we obtain the "parallel" version of the algorithm. The multiplicative form converges more rapidly than the parallel one but prevents us from solving subproblems in parallel (this problem can, however, be circumvented when we consider more than two subdomains). The interested readers may refer to [7] for further details regarding the different variants of the Schwarz method. Although the present study uses the multiplicative form of the algorithm, the theoretical results regarding the determination of optimized transmission conditions are also valid for the parallel form. Note that the usual algorithmic approach used in ocean-atmosphere climate models as described in [4] generally corresponds to one (and only one) iteration of the algorithm in (2.2) (with [F.sub.j] = [G.sub.j] = [D.sub.j](0)[[partial derivative].sub.x], j = 1, 2).

The primary role of the operators [F.sub.j] and [G.sub.j] (j = 1, 2) in (2.2) is to ensure a given consistency of the solution on the interface r. In our context we require the equality of the subproblems solutions and of their fluxes. The most natural choice to obtain such a connection consists in choosing

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

However, as proposed in [19], the same consistency can be obtained using mixed boundary conditions of Robin type, leading to

(2.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The advantage of (2.3) is that if the operators [[LAMBDA].sub.1] and [[LAMBDA].sub.2] are correctly chosen, then we can greatly improve the convergence speed of the corresponding algorithm [12]. Note that [[LAMBDA].sub.j] must also be carefully chosen to ensure the well-posedness of the problem. In this paper we focus on Robin-type transmission conditions since Dirichlet-Neumann-type algorithms generally converge quite slowly except for large discontinuities between the coefficients [D.sub.2] and [D.sub.1]. (It can easily be shown that the convergence rate is given by the square root of the ratio between [D.sub.1] and [D.sub.2].)

At this point, we have formulated the coupling problem we want to address. The convergence properties of this kind of problem have been extensively studied in the case of constant and continuous diffusion coefficients [8]. There are also some results in the case of constant and discontinuous coefficients [10] in the more general case of an advection-diffusion-reaction problem. This latter study provides results for specific asymptotic cases that are discussed later in Section 4.4. In this paper, we propose to investigate the problem with diffusion coefficients being constant in each subdomain and discontinuous at the interface, i.e., [D.sub.j] (x) = [D.sub.j], with [D.sub.j] > 0 and [D.sub.1] [not equal to] [D.sub.2]. We prove convergence of the algorithm (2.2) and we determine optimal choices for the operators [[LAMBDA].sub.j] under some constraints on the parameters of the problem.

2.2. Convergence of the algorithm. A classical approach to demonstrate the convergence of algorithm (2.2) consists of introducing the error [e.sup.k.sub.j] between the exact solution u* and the iterates [u.sup.k.sub.j] = 1, 2. By linearity, those errors satisfy homogeneous diffusion equations with homogeneous initial conditions. We denote the Fourier transform in time by [??] = F(g) for any g [member of] [L.sup.2] (R). Assuming that T [right arrow] [infinity] and that all the functions are equal to zero for negative times, it can easily be shown that the errors [[??].sup.k.sub.j] in Fourier space satisfy a second- order ordinary differential equation in x

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with characteristic roots [[sigma].sup.[+ or -].sub.j] = [+ or -] [square root of [absolute value of [omega]] / 2Dj] (1 + [[absolute value of [omega]] / [omega]]i). Note that the particular case [omega] = 0 would correspond to the existence of a stationary part in the error. However, since the error is initially zero, such a stationary part is also necessarily zero. To study the convergence of algorithm (2.2), it is usually assumed that [L.sub.1], [L.sub.2] [right arrow] [infinity], thus leading to

(2.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The validity of this assumption is discussed in [17]. The functions [alpha]([omega]) and [beta]([omega]) are determined using the Robin interface conditions at x = 0

(2.5) ([D.sub.1] [[sigma].sup.+.sub.1] + [[lambda].sub.1]) [[alpha].sup.k] (w) = ([D.sub.2] [[sigma].sup.-.sub.2] + [[lambda].sub.1]) [[beta].sup.k-1] (w),

(-[D.sub.2] [[sigma].sup.-.sub.2] + [[lambda].sub.2]) [[beta].sup.k] (w) = (-[D.sub.1] [[sigma].sup.+.sub.1] + [[lambda].sub.2]) [[alpha].sup.k] (w),

where [[lambda].sub.j] is defined as the symbol of the operator [A.sub.j] (j = 1, 2). A convergence factor [rho] of the Schwarz algorithm (2.2) can be defined as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Given (2.4), this amounts to [rho]([omega]) = [absolute value of [[alpha].sup.k] / [[alpha].sup.k-1]] = [absolute value of [[beta].sup.k]/[[beta].sup.k-1]]. Using (2.5) we obtain

(26) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

A more general derivation of the convergence factor for the case of an advection-diffusion-reaction problem with discontinuous coefficients can be found in [10]. At this point, we are not able to infer the convergence (or the divergence) of the corresponding algorithm because the operators [[LAMBDA].sub.j] have not been explicitly determined. It is often a difficult task to choose them in an appropriate way. The main difficulty comes from the fact that the convergence factor is formulated in the Fourier space, meaning that we can only act on symbols [[LAMBDA].sub.j] and not directly on pseudo-differential operators [[LAMBDA].sub.j] in the physical space.

2.3. The optimized Schwarz method. It is possible to find values [[lambda].sub.1] and [[lambda].sub.2] canceling the convergence factor (2.6) and therefore ensuring a convergence in exactly two iterations. Their expressions are

(2.7) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

These symbols correspond to so-called absorbing conditions. Unfortunately, since these optimal symbols are not polynomials in i[omega], the absorbing conditions are nonlocal in time in the physical space. The problem is thus to determine local operators providing a good approximation of nonlocal ones by finding a polynomial form in i[omega] to approximate [[lambda].sup.opt.sub.j]. There are mainly two approaches for such an approximation [12]. The first approach is a low frequency approximation, namely a Taylor expansion for a small [omega]. We decided not to adopt this approach because we want to be able to consider a wide range of frequencies. The second and more sophisticated approach is to solve a minmax problem to determine local operators that optimize the convergence speed over the full range of admissible frequencies [[[omega].sub.min], [[omega].sub.max]]. For a zeroth-order approximation, we look for values [[lambda].sup.0.sub.j] [member of] R such that [[lambda].sup.0.sub.j] [approximately equal to] [[lambda].sup.opt.sub.j]. The numbers [[lambda].sup.0.sub.j] can be defined as the solution of the optimization problem

(2.8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Since we work on a discrete problem, the frequencies allowed by our temporal grid range from [[omega].sub.min] = [pi]/T to [[omega].sub.max] = [pi]/[DELTA]t, where [DELTA]t is the time step of the temporal discretization. The analytical solution of problem (2.8) is not an easy task: the minimization of a maximum is known to be one of the most difficult problems in optimization theory [5]. Moreover, we are faced with an optimization problem for two parameters [[lambda].sup.0.sub.1] and [[lambda].sup.0.sub.2] which substantially enlarges the difficulty. Some analytical results exist in the case of a two-sided optimization for the 2D stationary diffusion equation [6, 20] and for the 2D Helmholtz equation [11]. In [10], the asymptotic solution of (2.8) for an advection-diffusion-reaction problem for [DELTA]t [right arrow] 0, [w.sub.min] = 0, and a positive advection is found in two particular cases: first for [[lambda].sup.0.sub.1] = [[lambda].sup.0.sub.2] (one-sided) and second for [[lambda].sup.0.sub.1] [not equal to] [[lambda].sup.0.sub.2] (two-sided) but [D.sub.1] = [D.sub.2]. In this paper, we study the complete minmax problem (2.8) in the general case [[lambda].sup.0.sub.1] [not equal to] [[lambda].sup.0.sub.2] and [D.sub.1] [not equal to] [D.sub.2]. Solving numerically the minmax problem (2.8) is quite expensive from a computational point of view. Moreover, this optimization must be performed for any change in the values of [D.sub.1] and [D.sub.2]. That is why we look for an analytical solution in the case of a zeroth-order approximation of the absorbing conditions. This is done with two different sets of interface conditions, the Neumann-Robin case and the Robin-Robin case.

Algorithm (2.2) with two-sided Robin conditions is well-posed for any choice of [[lambda].sup.0.sub.1] and [[lambda].sup.0.sub.2] such that [[lambda].sup.0.sub.1] + [[lambda].sup.0.sub.2] > 0. This result can be shown following the methodology based on a priori energy estimates as described in [1] and [8].

3. The optimized Schwarz method with Neumann-Robin interface conditions. In this section, we assume that the solution in [[OMEGA].sub.2] is subject to a Neumann boundary condition. The convergence speed of the Neumann-Robin algorithm is expected to be slower than that one obtained by a Robin-Robin algorithm. However, this easier case is treated explicitly because it introduces several methodological aspects useful for the determination of the general Robin-Robin optimized interface conditions. Imposing a Neumann boundary condition on the solution [u.sub.2] on [GAMMA] corresponds to having [[LAMBDA].sub.2] =0 in (2.3). The convergence factor [[rho].sub.NR] (NR stands for "Neumann-Robin") obtained from (2.6) reduces to

(3.1) = [[rho].sub.NR] = [absolute value of [D.sub.1][[sigma].sup.+.sub.1] ([D.sub.2][[sigma].sup.-.sub.2] + [[lambda].sub.1]) / [D.sub.2][[sigma].sup.-.sub.2] ([D.sub.1][[sigma].sup.+.sub.1] + [[lambda].sub.1])]

THEOREM 3.1 (Optimized Robin parameter). The analytical solution [[lambda].sup.0,*.sub.1] of the mini-max problem

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Proof. Introducing [zeta] = [square root of [absolute value of w][D.sub.1]], [gamma] = [square root of [D.sub.2]/[D.sub.1]], [[lambda].sup.0.sub.1] = ([square root of [[zeta].sub.min] [[zeta].sub.max]/2]) p, for p [member of] R, and making [[sigma].sup.+.sub.1] and [[sigma].sup.-.sub.2] in (3.1) explicit, we obtain

[[rho].sub.NR] (p, [zeta]) = 1 / [gamma] [square root of [(p - [gamma][zeta]).sup.2] + [[gamma].sup.2][[zeta].sup.2] / [(p + [zeta]).sup.2] + [[zeta].sup.2]],

with [zeta] = [zeta]/ [square root of [[zeta].sub.max] [[zeta].sub.min]]. Moreover, to ensure the well-posedness of the algorithm, we consider [[lambda].sup.0.sub.1] > 0, i.e., p > 0. Defining an additional parameter [mu] = [square root of [[zeta].sub.max] [[zeta].sub.min]], we obtain that [zeta] is bounded by [[zeta].sub.min] = [[mu].sup.-1] and [[zeta].sub.max] = [mu]. The aim is to optimize the convergence speed by finding p*, the solution of the minimax problem

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We first study the behavior of the derivative of [[rho].sub.NR] with respect to [zeta] and p with [zeta] [greater than or equal to] 0 and p [greater than or equal to] 0. For simplicity we introduce q = p ([gamma] - 1 + [square root of 1 + [[gamma].sup.2]]).

We first derive some restriction of the parameter's range. We can easily show that

(3.2) Sign ([partial derivative][[rho].sub.NR] / [[partial derivative].sub.p] = Sign (q - [zeta]).

Looking at the sign of the derivative of [[rho].sub.NR] with respect to p, we see that for all values of [zeta], the convergence factor [[rho].sub.NR] is a decreasing function of p for q < [[zeta].sub.min] = [[mu].sup.- 1], proving that q* [greater than or equal to] [[zeta].sub.min]. A similar argument shows that q* [less than or equal to] [[zeta].sub.max]. This proves that the optimized parameter q* satisfies

1 / [mu] [less than or equal to] q* [less than or equal to] [mu].

Along with (3.2), this shows that the convergence factor has to be an increasing function of p at [zeta] = 1 /[mu] and a decreasing function of p at [zeta] = [mu].

Next we show an equioscillation property of the optimal parameter. The sign of the derivative of [[rho].sub.NR] with respect to [zeta] is given by

Sign ([partial derivative][[rho].sub.NR] / [[partial derivative].sub.[zeta]]) = Sign ([zeta] - q).

This relation implies that [[rho].sub.NR] has a local minimum between 1 / [mu] and [mu]. The maximum value of the convergence factor is thus attained either at [zeta] = 1/[mu] or at [zeta] = [mu] (or both). If we assume [[rho].sub.NR] (p, 1/[mu]) < [[rho].sub.NR] (p, [mu]), it is always possible to decrease the maximum value of [[rho].sub.NR] (p, [zeta]) by increasing the value of p so that we have [[rho].sub.NR](p, 1/[mu]) [greater than or equal to] [[rho].sub.NR] (p, [mu]). A similar argument shows that [[rho].sub.NR] (p, [mu]) [greater than or equal to] [[rho].sub.NR] (p, 1/[mu]). The optimal parameter must thus satisfy the equioscillation property [[rho].sub.NR] (p*, 1/[mu]) = [[rho].sub.NR] (p*,[mu]). After some simple algebra, we find that p* is a solution of

([gamma] - 1)([mu] +1/[mu]) + 2[gamma]/p* - p* = 0.

The unique positive solution of the equation v* = 2[gamma]/p* - p* with v* = (1 - [gamma]) ([mu] + 1/[mu]) is given by p* = 1/2 (-v* + [square root of 8[gamma] + [(v*).sup.2]]). After a substitution of [gamma] and [mu] and a multiplication of p* by [square root of [[zeta].sub.min] [[zeta].sub.max] /2], we retrieve the stated result for [[lambda].sup.0,*.sub.1].

We find that the optimized convergence factor satisfies an equioscillation property. This concept of equioscillation property comes from the Chebyshev's alternation theorem (or equioscillation theorem). The similarities between the Chebyshev's theorem and the optimized Schwarz method are clearly exposed in [2, 6]. Two typical optimized convergence factors [[rho].sup.*.sub.NR] = [[rho].sub.NR] ([[lambda].sup.0,*.sub.1]) are shown in Figure 3.1 (left) for [mu] = 2 and [mu] = 6 with [gamma] = 5. Note that the performance of the optimized algorithm is only a function of the ratio [gamma] between [D.sub.1] and [D.sub.2] and not of the actual values of those parameters. The same remark applies to the temporal frequencies [[omega].sub.min] and [[omega].sub.max]: [[rho].sup.*.sub.NR] is only a function of their ratio [mu].

It is also instructive to look at three particular cases: [gamma] [right arrow] [0.sup.+], [gamma] = 1, and [gamma] [right arrow] [infinity].

* [gamma] [right arrow] [0.sup.+] ([D.sub.1] >> [D.sub.2]):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The minimum value of the convergence factor is attained at [mu] = 1 and is equal to [square root of 2]/2. When [mu] is increased, the convergence is very slow. Indeed, we tend towards a Neumann-Neumann algorithm in this case.

* [gamma] = 1 ([D.sub.1] = [D.sub.2] = D):

[[rho].sup.*.sub.NR] = [square root of 1 - 2 [square root of 2][mu] / 1 + [mu] ([mu] + [square root of 2])], [[lambda].sup.0,*/1] [([[omega].sub.max][ [omega].sub.min]).sup.1/4].

The convergence factor [[rho].sup.*.sub.NR] approaches 1 when [mu] is increased. One can also remark that the optimal parameter [[lambda].sup.0,*.sub.1] is the same as that one found in [8] in the Robin-Robin one-sided case.

* [gamma] [right arrow] +[infinity] ([D.sub.1] << [D.sub.2]):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

When [gamma] tends to +[infinity], the convergence is very fast (the convergence factor approaches 0) and the optimal boundary condition tends towards a Neumann-Dirichlet operator.

These results are illustrated in Figure 3.1 (right). The efficiency of the Neumann-Robin algorithm is greatly improved when [gamma] becomes large and [mu] becomes small. We continue this section by studying the asymptotic convergence rate for the discretized algorithm when the time step [DELTA]t tends to 0.

THEOREM 3.2 (Asymptotic performance). For [D.sub.2] > [D.sub.1] (i.e., [gamma] > 1), [[omega].sub.max] = [pi]/[[DELTA].sub.2] and for [DELTA]t tending to zero, the optimal Robin parameter given by Theorem 3.1 is

[[lambda].sup.0,*.sub.1] [approximately equal to] [square root of 2[D.sub.1]] ([gamma] - 1 / 2 [square root of [pi][DELTA][t.sup.-1/2]] + [[gamma].sup.2] + 1 / 2([gamma] - 1) [square root of [[omega].sub.min]])

and the asymptotic convergence of the optimized Neumann-Robin algorithm is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We conclude that the zeroth-order optimized Neumann-Robin boundary conditions are efficient when the Robin condition is imposed at the boundary of the domain with the smaller diffusion coefficient ([[OMEGA].sub.1] here). In this case, the asymptotic convergence factor [[rho].sup.*.sub.NR] is of the form [square root of [D.sub.1]/[D.sub.2]] (1 - [OMICRON]([DELTA][t.sup.1/2])) for small [DELTA]t. In the next section, we study the zeroth-order two-sided Robin-Robin boundary conditions.

4. The optimized Schwarz method for a diffusion problem with discontinuous (but constant) coefficients: two-sided Robin transmission conditions. In this section we optimize the conditions on both sides of the interface to get a faster convergence speed regardless of the value of [gamma]. By keeping the notations [zeta], [zeta], [mu], and [gamma] defined in the previous section and by approximating [[lambda].sup.opt.sub.1] and [[lambda].sup.opt.sub.2] respectively by [[lambda].sup.0.sub.2] = [square root of [[zeta].sub.min][[zeta].sub.max] /2] [p.sub.2] and [[lambda].sup.0.sub.2] = [square root of [[zeta].sub.min][[zeta].sub.max]] [p.sub.1], the convergence factor [[rho].sub.RR] reads

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We can easily demonstrate that for nonnegative fixed values of [zeta] and [gamma] and for [p.sub.1], [p.sub.2] > 0, we find the three inequalities [[rho].sub.RR] ([p.sub.1], [p.sub.2], [zeta]) < [[rho].sub.RR] (-[p.sub.1], -[p.sub.2], [zeta]), [[rho].sub.RR] ([p.sub.1],[p.sub.2], [zeta]) < [[rho].sub.RR] ([p.sub.1], -[p.sub.2], [zeta]), and [[rho].sub.RR] ([p.sub.1], [p.sub.2], [zeta]) < [[rho].sub.RR] (-[p.sub.1], [p.sub.2], [zeta]). This shows that we can restrict our study to strictly positive values of [p.sub.1] and [p.sub.2] (note that [p.sub.1] = 0 or [p.sub.2] = 0 corresponds to the Neumann-Robin case). The restriction of the parameter range to strictly positive values ensures that [[lambda].sup.0.sub.1] + [[lambda].sup.0.sub.2] > 0, and hence the corresponding problem is well-posed. In the following, we assume that [gamma] [greater than or equal to] 1. The optimal parameters [p.sub.1] and [p.sub.2] for the case [gamma] [less than or equal to] 1 can be obtained by switching (i.e. [p.sub.1] becomes [p.sub.2] and [p.sub.2] becomes [p.sub.1]) the optimal values for the case [gamma] [greater than or equal to] 1. As it was done previously, we choose the values [p.sub.1] and [p.sub.2] by solving the optimization problem

(4.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

4.1. Behavior of the convergence factor with respect to the Robin parameters. First, we study the behavior of [[rho].sub.RR] with respect to the parameters [p.sub.1] and [p.sub.2]. We introduce two new parameters [q.sub.1] and [q.sub.2] defined by

[q.sub.1] = [p.sub.1] / 1 - [gamma] + [square root of 1 + [[gamma].sup.2]] and [q.sub.2] = [p.sub.2] / [gamma] - 1 + [square root of 1 + [[gamma].sup.2]]

We can demonstrate that for [gamma] [greater than or equal to] 1 and [q.sub.1] [less than or equal to] [q.sub.2], we have that [[rho].sub.RR] ([p.sub.1], [p.sub.2], [zeta]) [less than or equal to] [[rho].sub.RR] ([p.sub.2], [p.sub.1], [zeta]). This proves that the optimal parameters satisfy [q*.sub.1] < [q*.sub.2]. This implies that in turn [p.sub.1] [less than or equal to] [p.sub.2] and that [p.sub.1] < [p.sub.2] if [gamma] > 1. This immediately proves that one-sided ([p.sub.1] = [p.sub.2]) Robin-Robin boundary conditions are not optimal as soon as [gamma] > 1.

Note that Sign ([partial derivative][[rho].sub.RR] /[partial derivative][p.sub.1]) = Sign([q.sub.1] - [zeta]) and Sign ([partial derivative][[rho].sub.RR] / [partial derivative][p.sub.2]) = Sign([q.sub.2] - [zeta]) implies

(4.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Looking at the signs of the derivatives of [[rho].sub.RR] with respect to [p.sub.1] and [p.sub.2], we find that if we choose [q.sub.1] < [[zeta].sub.min] = [[mu].sup.-1], we can decrease the convergence factor by increasing [p.sub.1] because [partial derivative][[rho].sub.RR] / [partial derivative][p.sub.1] < 0 holds for all [q.sub.1] > [[zeta].sub.min]. A similar argument shows that [q.sub.2] [less than or equal to] [[zeta].sub.max]. This means that the optimized parameters [q*.sub.1] and [q*.sub.2] must satisfy

(4.3) [[mu].sup.-1] [less than or equal to] [q*.sub.1] < [q*.sub.2] [less than or equal to] [mu].

The inequalities (4.2) and (4.3) imply that at [zeta] = 1/[mu] is an increasing function of [p.sub.1] and [p.sub.2] (or [q.sub.1] and [q.sub.2]) while at [zeta] = [mu], [[rho].sub.RR] is a decreasing function of [p.sub.1] and [p.sub.2] (or [q.sub.1] and [q.sub.2]).

4.2. Extrema of [[rho].sub.RR] with respect to [zeta]. The next step to solve (4.1) analytically is to find the location of the extrema of [[rho].sub.RR] ([p.sub.1], [p.sub.2], [zeta], [gamma]) with respect to [zeta].

Theorem 4.1 (Extrema of [[rho].sub.RR] ([zeta])). The function [zeta] [right arrow] [[rho].sub.RR] ([p.sub.1],[p.sub.2], [zeta]) has one or three positive local extrema. In the case of one extremum, it corresponds to a minimum and is located at [zeta] = x := [square root of [p.sub.1][p.sub.2] / 2[gamma]].

Proof. We first introduce the following property that can easily be verified:

[[rho].sub.RR] ([p.sub.1], [p.sub.2], [zeta]) = [[rho].sub.RR] ([p.sub.1], [p.sub.2]/[zeta]), where X = [square root of [p.sub.1][p.sub.2] / 2[gamma]].

Differentiating with respect to [zeta] leads to

(4.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which shows that [partial derivative][[rho].sub.RR] / [partial derivative][zeta] ([p.sub.1], [p.sub.2], [+ or -]x) = 0. The derivative [partial derivative][[rho].sub.RR]([p.sub.1], [p.sub.2], [zeta]) / [partial derivative][zeta] has the same sign as a (unitary) sixth-order polynomial P([zeta]) (the full expression of P is complicated and not given here). P([zeta]) has thus either two or six real roots, among them [zeta] = X is positive and [zeta] = -x is negative. Let us suppose that P([zeta]) has six real roots. We can show that only three of these six roots (including [zeta] = x) are positive. From (4.4) we see that if [[zeta].sup.0] is a root of P([zeta]), then [[zeta].sup.1] = [chi square]/[[zeta].sup.0] is another one. Assuming that the four other roots are positive, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and the sum of the six roots must be greater than [2.sub.x] and is therefore positive. However, the sum of the six roots of P([zeta]) is given by -[a.sub.5] where a5 is the coefficient of the term [[zeta].sup.5] and is equal to [a.sub.5] = ([gamma]-1)([p.sub.2]-[p.sub.1]) / [gamma]. Using the facts [gamma] [greater than or equal to] 1 and [p.sub.2] [greater than or equal to] [p.sub.1] (from (4.3)), -[a.sub.5] cannot be positive. We conclude that we have at most three positive roots for P([zeta]). It can be verified that P(0) < 0 and P(+[infinity]) > 0 so that if only one positive root exists (at [zeta] = x), it is a local minimum.

4.3. Equioscillation of [[rho].sub.RR] at the end points.

Theorem 4.2 (Equioscillation at the end points). The optimized convergence factor [[rho].sub.RR] ([p*.sub.1], [p*.sub.2], [zeta]) satisfies

* [[rho].sub.RR] ([p*.sub.1], [p*.sub.2], x) [less than or equal to] max ([[rho].sub.RR] ([p*.sub.1], [p*.sub.2], [[mu].sup.-1]), [[rho].sub.RR] ([p*.sub.1], [p*.sub.2], [mu])) with x = [square root of [p*.sub.1][p*.sub.2] / 2[gamma]],

* the equioscillation property [[rho].sub.RR] ([p*.sub.1], [p*.sub.2], [[mu].sup.-1]) = [[rho].sub.RR] ([p*.sub.1], [p*.sub.2], [mu]), which holds only for [p*.sub.1] [p*.sub.2] = 2[gamma].

Proof. We first show that [[rho].sub.RR] ([p*.sub.1], [p*.sub.2], x) [less than or equal to] max ([[rho].sub.RR]([p*.sub.1], [p*.sub.2], [[mu].sup.-1]), [[rho].sub.RR] ([p*.sub.1], [p*.sub.2], [mu])). This is straightforward when x is the only positive root of [partial derivative][[rho].sub.RR] / [partial derivative][zeta] because x is a local minimum. In the case when there are three positive roots, x is a local maximum. From the identity x = [square root of [p.sub.1][p.sub.2] / 2[gamma]] = [square root of [q.sub.1][q.sub.2]] and (4.3), we get

(4.5) 1/[mu] [less than or equal to] [q.sub.1] [less than or equal to] x = [square root of [q.sub.1][q.sub.2]] [less than or equal to] [q.sub.2] [less than or equal to] [mu].

We already know that at [zeta] = 1/[mu], [[rho].sub.RR] is a decreasing function of [q.sub.1] and that at [zeta] = [mu], [[rho].sub.RR] is an increasing function of [q.sub.1]. The inequality (4.5) shows that at [zeta] = x, [[rho].sub.RR] is an increasing function of [q.sub.1] since [q.sub.1] [less than or equal to] x. If we assume that [[rho].sub.RR] ([p*.sub.1], [p*.sub.2], x) [greater than or equal to] [[rho].sub.RR] ([p*.sub.1], [p*.sub.2], [[mu].sup.-1]), then we can always decrease [q.sub.1] (or [p.sub.1]) such that it improves the convergence factor (by reducing the values both at [zeta] = x and at [zeta] = [mu]). Playing with [q.sub.2], we can similarly prove that [[rho].sub.RR] ([p*.sub.1], [p*.sub.2], x) [less than or equal to] [[rho].sub.RR] ([p*.sub.1], [p*.sub.2], [mu]). Note that this also proves that [[zeta].sub.1] > 1/[mu] and [[zeta].sub.3] [less than or equal to] [mu]. This is sufficient to fully describe the behavior of the convergence factor with respect to [q.sub.1], [q.sub.2], and z, as shown in Figure 4.1. In practice, the two cases will differ by the sign of the second-order derivative of [[rho].sub.RR] ([p.sub.1], [p.sub.2], [zeta]) at [zeta] = x. The following argument proves that the values taken by [[rho].sub.RR] ([p.sub.1],[p.sub.2], [zeta]) at the two end points [zeta] = 1/[mu], and [zeta] = [mu] are equal. Indeed, if we assume that [[rho].sub.RR] ([p.sub.1], [p.sub.2], 1/[mu]) < [[rho].sub.RR] ([p.sub.1], [p.sub.2], [mu]) (or [[rho].sub.RR] ([p.sub.1], [p.sub.2], 1/[mu]) > [[rho].sub.RR] ([p.sub.1], [p.sub.2], [mu])) holds, it is always possible to decrease the maximum value of [[rho].sub.RR]([zeta]) by increasing (respectively decreasing) the values of [p.sub.1] (respectively [p.sub.2]). The optimal parameters must thus satisfy the equioscillation property [[rho].sub.RR] ([p*.sub.1], [p*.sub.2], [[mu].sup.-1]) = [[rho].sub.RR]([p*.sub.1], [p*.sub.2], [mu]). This holds when

(4.6) ([p.sub.1] + [p.sub.2])(2[gamma] - [p.sub.1][p.sub.2])S([p.sub.1],[p.sub.2], [mu], [gamma]) = 0,

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Obviously every pair ([p.sub.1], [p.sub.2]) that satisfies the relation [p.sub.1][p.sub.2] = 2[gamma] is a solution to (4.6). We now show that there are no other admissible values. Other potential solutions of the problem are the solutions of S([p.sub.1],[p.sub.2], [mu]) =0. S is a second-order polynomial in [p.sub.2] and thus has two real solutions:

(4.7) [p.sub.2] = [f.sub.1] ([p.sub.1]), [p.sub.2] = [f.sub.2]([p.sub.1]).

If we assume that [p.sub.2] is related to [p.sub.1] with one of the relations (4.7), looking at Figure 4.1, we can argue that for any pair ([p.sub.1], [p.sub.2]) we must have d[p.sub.2] / d[p.sub.1] < 0 to satisfy an equioscillation property. Indeed, let [[rho].sup.[dagger].sub.RR]([p.sub.1], [zeta]) be defined as

[[rho].sup.[dagger].sub.RR]([p.sub.1], [zeta]) = [[rho].sub.RR] ([p.sub.1], [p.sub.2] ([p.sub.1]), [zeta])

Then

(4.8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We have already proved that the following properties must hold

(4.9) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

If we suppose d[p.sub.2] / d[p.sub.1] > 0, then (4.8) and (4.9) show that [[rho].sup.[dagger].sub.RR]([p.sub.1],1/[mu]) is an increasing function of [p.sub.1] while [[rho].sup.[dagger].sub.RR] ([p.sub.1], [mu]) is a decreasing function of [p.sub.1]. Hence, (4.9) and the equioscillation property cannot be satisfied at the same time if d[p.sub.2]/d[p.sub.1] > 0. It can be shown that the two solutions given by (4.7) do not satisfy this latter condition. Indeed, we can prove that we have d[f.sub.i]/d[p.sub.1] > 0 and d[f.sub.2] / d[p.sub.1] > 0. Details of the computations are omitted here but we mention that the only conditions necessary to find this result are [gamma] > 0, [mu] > 1. We can conclude that [p.sub.1] [p.sub.2] = 2[gamma] is the only solution leading to an equioscillation property.

It is worth mentioning that [chi] = [square root of [p.sub.1][p.sub.2]/2[gamma]] = 1 and that

[[rho].sub.RR] ([p*.sub.1], [p*.sub.2], [zeta]) = [[rho].sub.RR] ([p*.sub.1], [p*.sub.2], 1/[zeta]), [for all][zeta] [member of] [1/[mu], [mu]].

4.4. Solution of the minmax problem. The convergence factor is now a function of [p.sub.1] and [zeta] only:

[[rho].sup.[dagger].sub.RR]([p.sub.1], [zeta]) = [[rho].sub.RR] ([p.sub.1], 2[gamma]/[p.sub.1], [zeta]).

Lemma 4.3. The solution of the minmax problem (4.1) is given by the solution of the constraint minimization problem

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [p.sup.*,equi.sub.1] is the solution of the three-point equioscillation problem

[[rho].sup.[dagger].sub.RR]([p.sub.1], 1) = [[rho].sup.[dagger].sub.RR] ([p.sub.1], 1/[mu]) = [[rho].sup.[dagger].sub.RR]([p.sub.1], [mu]).

Proof. Thanks to Figure 4.1, we can remark that the resolution of the minmax problem corresponds to the minimization of [[rho].sup.[dagger].sub.RR]([p.sub.1], 1/[mu]) (or [[rho].sup.[dagger].sub.RR]([p.sub.1], [mu])) with respect to [p.sub.1]. In the case where x =1 is a local maximum, the additional constraint given by Theorem 4.2 must be imposed

(4.10) [[rho].sup.[dagger].sub.RR]([p.sub.1], 1) [less than or equal to] [[rho].sup.[dagger].sub.RR]([p.sub.1], 1/[mu])

Knowing that [p.sub.1][p.sub.2] = 2[gamma] or equivalently [q.sub.1][q.sub.2] = 1, the range of admissible values given by the inequality (4.3) can now be written as 1 / [mu] [less than or equal to] [q.sub.1] [less than or equal to] 1 and translates in terms of the variable [p.sub.1] :

(4.11) [p.sub.1] [member of] [[p.sub.1,min], [p.sub.1,max]], where

[p.sub.1,min] = (1 - [gamma] + [square root of 1 + [[gamma].sup.2])/[mu], [p.sub.1,max] = (1 - [gamma] [square root of 1 + [[gamma].sub.2]).

Moreover, it can be shown that [[rho].sup.[dagger].sub.RR]([p.sub.1], 1) is a decreasing function of [p.sub.1] and therefore the constraint (4.10) can also be written as [p*.sub.1] [greater than or equal to] [p.sup.*,equi.sub.1] where [p.sup.*,equi.sub.1] is the solution of a three-point equioscillation problem [[rho].sup.[dagger].sub.RR]([p.sup.*,equi.sub.1], 1) = [[rho].sup.[dagger].sub.RR]([p.sup.*,equi.sub.1], 1/[mu]) (=[[rho].sup.[dagger].sub.RR]([p.sup.*,equi.sub.1], [mu]))).

We now look at the minimization of [[rho].sup.[dagger].sub.RR]([p.sub.1], 1/[mu]) for [p.sub.1] [member of] [[p.sub.1,min], [p.sub.1,max]].

Lemma 4.4. For [gamma] > 1, the derivative of [[rho].sup.[dagger].sub.RR]([p.sub.1], 1/[mu]) has exactly one root in the range [[p.sub.1,min], [p.sub.1,max]]. This root corresponds to a local minimum of [[rho].sup.[dagger].sub.RR] ([p.sub.1], 1/[mu]). In the special case [gamma] = 1, [p.sub.1] = [p.sub.1,max] (= [square root of 2]) is always a root of [partial derivative][[rho].sup.[dagger].sub.RR] / [partial derivative][p.sub.1] ([p.sub.1],1/[mu]).

Proof. The derivative of [[rho].sup.[dagger].sub.RR] ([p.sub.1], 1/[mu]) can be written as

[partial derivative][[rho].sup.[dagger].sub.RR] / [partial derivative][p.sub.1] ([p.sub.1], 1/[mu]) = g([p.sub.1], [mu]) N ([p.sub.1], [mu]),

where g is a strictly positive function and N([p.sub.1], [mu]) a sixth-order polynomial in [p.sub.1]. The change of variable v = 2[gamma]/[p.sub.1] - [p.sub.1] transforms N([p.sub.1], [mu]) into

N ([p.sub.1], [mu]) = [p.sup.3.sub.1] Q(v),

where Q(v) is the third-order polynomial given by

(4.12) Q(v) = 8([gamma] - 1) (1 + [[gamma].sup.2]) + 2[beta]([gamma][[beta].sup.2] - 3(1 + [[gamma].sup.2])) v + 2(y - 1) [[beta].sup.2][v.sup.2] - [beta][v.sup.3],

with [beta] = 1/[mu] + [mu].

It can be shown that, for [gamma] > 1, this polynomial has only one root in [[v.sub.min], [v.sub.max]], where, according to (4.11), [v.sub.min] and [v.sub.max] are given by

[v.sub.min] = 2([gamma] - 1), [v.sub.max] = ([gamma] - 1) [beta] + [square root of 1 + [[gamma].sup.2]] [square root of [[beta].sup.2] - 4].

This root corresponds to a minimum of the functional [[rho].sup.[dagger].sub.RR]([p.sub.1], 1/[mu]) because we can show that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. For [gamma] = 1, the value v = [v.sub.min] = 0, i.e., [p.sub.1] = [p.sub.1,max] = [square root of 2], is always a root of Q(v). Figure 4.2 illustrates the variations of [[rho].sup.[dagger].sub.RR]([p.sub.1], 1/[mu]) with respect to [p.sub.1]. [p.sup.min.sub.1] is the location of the minimum of [[rho].sup.[dagger].sub.RR]([p.sub.1], 1/[mu]) over the interval [[p.sub.1,min], [p.sub.1,max]]. The solution of the constrained minimization problem is now easily handled: if [p.sup.min.sub.1] [less than or equal to] [p.sup.*,equi.sub.1] then the solution of the minmax problem is given by [p.sup.*,equi.sub.1], otherwise the solution of the minmax problem is given by [p.sup.min.sub.1] .

The inequality [p.sup.min.sub.1] [less than or equal to] [p.sup.*,equi.sub.1] is satisfied if and only if [partial derivative][[rho].sup.[dagger].sub.RR] / [partial derivative][p.sub.1] ([p.sup.*,equi.sub.1], [mu]) [greater than or equal to] 0 or equivalently if Q([v.sup.*,equi]) [greater than or equal to] 0, where [v.sup.*,equi] = 2[gamma]/[p.sup.*,equi.sub.1] - [p.sup.*,equi.sub.1].

Finally, the following result is useful: for v [greater than or equal to] [v.sub.max] (or equivalently [p.sub.1] [less than or equal to] [p.sub.1,min]), we have Q(v) [less than or equal to] 0 ([MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]). Indeed, using the relation (4.8) at [zeta] = 1/[mu], we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

If [p.sub.1] [less than or equal to] [p.sub.1,min] holds, then we get [partial derivative][[rho].sub.RR] ([p.sub.1], [p.sub.2] ([p.sub.1]), 1/[mu]) / [partial derivative][p.sub.1] < 0, but the relation [p.sub.2] = 2[gamma]/[p.sub.1] implies that [p.sub.2] [greater than or equal to] [p.sub.2,max] (= ([gamma] - 1 + [square root of 1 + [[gamma].sup.2]]) [mu]) so that [partial derivative][[rho].sub.RR] ([p.sub.1], [p.sub.2] ([p.sub.1]), 1/[mu]) / [partial derivative][p.sub.2] [greater than or equal to] 0. With the help of d[p.sub.2]/d[p.sub.1] = -2[gamma]/[p.sup.2.sub.1] < 0, this proves that [partial derivative][[rho].sup.[dagger].sub.RR] ([p.sub.1], 1/[mu]) / [partial derivative][p.sub.1] [less than or equal to] 0.

We are now finished with the problem of finding the solution of the three-point equioscillation problem.

Theorem 4.5 (Equioscillation between three points). The only parameters [p.sup.*,equi.sub.1] and [p.sup.*,equi.sub.2], such that [p.sup.*,equi.sub.1] [less than or equal to] [p.sub.1,max], that satisfy an equioscillation of the convergence factor [[rho].sub.RR] between the three points (1 /[mu], 1, [mu]) are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where

(4.13) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Proof. We have to find the solution of the problem [[rho].sup.[dagger].sub.RR]([p.sub.1],1/[mu]) = [[rho].sup.[dagger].sub.RR]([p.sub.1],1). It can be shown that this is equivalent to the search for the roots of a fourth-order polynomial R([p.sub.1]) that can be written as

R([p.sub.1]) = [p.sup.2.sub.1]T (v), T (v) = 2(1 + [[gamma].sup.2]) - 4[gamma][beta] + (1 - [gamma])(2 + [beta])v + [v.sup.2],

where v is again defined by v = 2[gamma]/[p.sub.1] - [p.sub.1]. The unique root of T(v) that satisfies v [greater than or equal to] [v.sub.min] (i.e., [p.sub.1] [less than or equal to] [p.sub.1,max]) is given by

[v.sup.*,equi] = 1/2 [(2 + [beta])([gamma] - 1) + [square root of 4[(1 + [gamma]).sup.2] ([beta] - 1) + [[beta].sup.2] [([gamma] - 1).sup.2]],

and the expression for [p.sup.*,equi.sub.1] is deduced from the relation between [p.sub.1] and v.

Gathering the results, the solution of the minmax problem is given in the following theorem.

Theorem 4.6. The analytical solution [[lambda].sup.0,*.sub.1] and [[lambda].sup.0,*.sub.2] of the minmax problem

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with [v.sup.*,equi] given by (4.13) and [v.sup.*,mini] is the unique solution of Q(v) = 0 in [[v.sub.min], [v.sub.max]].

Proof. All the ingredients for the proof are stated before. Note that [v.sup.*,equi] may be larger than [v.sub.max]. However, since we have proved that Q(v [greater than or equal to] [v.sub.max]) [less than or equal to] 0, this case does not have to be considered explicitly. A substitution of [gamma] and [mu] by their respective expressions, and a multiplication of [p*.sub.1] and [p*.sub.2] by [square root of [[zeta].sub.min][[zeta].sub.max]/2] lead to the result for [[lambda].sup.0,*.sub.1] and [[lambda].sup.0,*.sub.2] with respect to [D.sub.1], [D.sub.2], [[omega].sub.min], and [[omega].sub.max].

Note that the following additional result can be shown as well:

(4.14) Q([v.sup.*,equi]) [greater than or equal to] 0 [??] [beta] [greater than or equal to] 1 + [square root of 5] or ([[beta].sub.0] < [beta] < 1 + [square root of 5] and [gamma] [greater than or equal to] f([beta])),

where [[beta].sup.0] is the root of the fourth-order polynomial 16 - 16X -4[X.sup.2] + [X.sup.4] whose approximate value is given by [[beta].sup.0] [approximately equal to] 2.77294 and f is given by

f ([beta]) = 1 / 16 - 16[beta] - 4[[beta].sup.2] + [[beta].sup.4] ([([beta] - 2).sup.3] [beta] ([beta] + 2)

+ (4 + 2[beta] - [[beta].sup.2]) [square root of -16 + 48[beta] - 44[[beta].sup.2] + 12[[beta].sup.3] + 3[[beta].sup.4] - 4[[beta].sup.5] + [[beta].sup.6]).

The function f([beta]) is plotted in Figure 4.3 for [[beta].sup.0] < [beta] < 1 + [square root of 5]. We remark that f([beta]) [greater than or equal to] 1 for all [beta] so that the condition [gamma] [greater than or equal to] f([beta]) is always false for [gamma] = 1 (the continuous case).

It is also interesting to know if x = [square root of [p.sub.1][p.sub.2] / 2[gamma]] = 1 is either a local minimum or a local maximum of the optimized convergence factor by looking at the sign of [[partial derivative].sup.2] [[rho].sup.[dagger].sub.RR] ([p.sub.1], x). It can be proved that in terms of the variable v = 2[gamma]/[p.sub.1] - [p.sub.1], the inequality [[partial derivative].sup.2] [[rho].sup.[dagger].sub.RR] ([p.sub.1], x) > 0 can be written as

v [greater than or equal to] [v.sub.0], where [v.sub.0] = 2([gamma] - 1) + [square root of 2(1 + [[gamma].sup.2])].

We deduce that [zeta] = x = 1 is a local minimum only if [v.sup.*,mini] [less than or equal to] [v.sub.0]. This can be verified by evaluating the polynomial Q(v) at v = [v.sub.0] and looking at the sign of the result: if Q([v.sub.0]) [less than or equal to] 0, then [v.sup.*,mini] [less than or equal to] [v.sub.0] and we have a local minimum at [zeta] = x = 1. It can be found that

Q([v.sub.0]) < 0 [??] 2 < [beta] < [[beta].sub.0] or ([[beta].sub.0] [less than or equal to] [beta] [less than or equal to] 2[square root of 2] and [gamma] < g([beta])), where [[beta].sub.0] = 8+5 [square root of 2]/2(3+2 [square root of 2]) + [square root of 90+64 [square root of 2]] / 2(3+2[square root of 2]) [approximately equal to] 2.44547. The analytical expression of g([beta]) is complicated and not given here. Note that g([beta]) [greater than or equal to] 1 for all [beta] so that for the special case [gamma] = 1, the inequality Q([v.sub.0]) < 0 is equivalent to 2 < [beta] < 2 [square root of 2].

Figure 4.4 summarizes the three different domains: three-point equioscillation, two-point equioscillation with x as a local maximum and two-point equioscillation with x as a local minimum. The resulting optimized convergence factor is shown in Figure 4.5 with respect to [mu] and [gamma].

We make the following remarks about the convergence properties of the Schwarz algorithms: the convergence speed increases when the discontinuities of the coefficients ([gamma]) is increased and the convergence speed decreases when [mu], an increasing function of the ratio [[omega].sub.max] / [[omega].sub.min], is increased. In Figure 4.6 we compare, for [mu] = 2 and [mu] = 6, the results found in the optimized two-sided case with the optimized Robin-Neumann transmission conditions (found in Section 3). The Robin-Robin approach is significantly more efficient than the Robin-Neumann approach when [gamma] is close to one. When [gamma] is increased, both tend towards a Dirichlet-Neumann operator.

THEOREM 4.7 (Asymptotic performance). For [D.sub.2] > [D.sub.1] (i.e., [gamma] > 1), [[omega].sub.max] = [pi]/[DELTA]t, and for [DELTA]t tending to zero, the optimal Robin parameters given by Theorem 4.6 are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and the asymptotic convergence of the optimized two-sided Robin-Robin algorithm is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Note that these asymptotic results are obtained by assuming that [v.sup.*] = [v.sup.*,equi], which is always the case when [DELTA]t [right arrow] 0 (i.e., [mu] [right arrow] [infinity]), as shown by (4.14). The optimized Robin-Robin conditions lead to an asymptotic convergence factor [square root of [D.sub.1]/[D.sub.2]] (1 - [OMICRON]([DELTA][t.sup.1/4])) for small [DELTA]t and [D.sub.1] < [D.sub.2]. The associated algorithm is thus less sensitive to [DELTA]t than the Neumann-Robin algorithm. However, the asymptotic Robin parameters given in Theorem 4.7 must be used with caution as they degenerate when [gamma] [right arrow] 1 as well as when [DELTA]t >> 0 (in this case [[lambda].sup.0,(as).sub.1] (as) can become negative). It is worth mentioning that the asymptotic bound on the optimized convergence factor given in Theorem 4.7 shows that the optimized Robin-Robin conditions will always be more efficient than Dirichlet-Neumann conditions. Indeed, it can easily be checked that the multiplicative factor 1/[gamma] in front of the bound corresponds to the convergence factor of the Dirichlet-Neumann algorithm.

Furthermore, we can not directly compare this result with the one obtained in [10] for the advection-diffusion-reaction equation. The latter study is done by assuming [[omega].sub.min] = 0 and as a result of this assumption their optimized parameters, when canceling the advection and reaction coefficients, are simply [[lambda].sup.0,*.sub.1] = [[lambda].sup.0,*.sub.2] = 0. Indeed, one can easily find that for a diffusion problem, the low frequency approximation [[lambda].sup.low.sub.j] of the absorbing conditions [[lambda].sup.opt.sub.j] given in (2.7) for [[omega].sub.min] [right arrow] 0 is indeed [[lambda].sup.low.sub.j] = 0.

4.5. The continuous case. Because the two-sided Robin-Robin case with continuous diffusion coefficients has never been studied in the literature, we now provide the results in this particular case.

THEOREM 4.8 (Continuous case). Under the assumption [D.sub.1] = [D.sub.2] = D, the optimal parameters [[lambda].sup.0,*.sub.1] and [[lambda].sup.0,*.sub.2] are given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Proof. We use Theorem 4.6, which gives the optimality conditions in the general case. As already mentioned, the condition Q([v.sup.*,equi]) [greater than or equal to] 0 reduces to [beta] [greater than or equal to] 1 + [square root of 5] for [gamma] = 1. In that case, the solution of the minmax problem is given by v* = [v.sup.*,equi] = 2 [square root of [beta] - 1]. If [beta] < 1 + [square root of 5], we have to compute [v.sup.*,min], the value that cancels Q(v) in [[v.sub.min], [v.sub.max]], where [v.sub.min] = 0, [v.sub.max] = 2 [square root of [[beta].sup.2] - 4]. For [gamma] = 1, the expression (4.12) of the polynomial Q(v) is

Q(v) = -[beta]v ([v.sup.2] - (2[[beta].sup.2] - 12)).

We find that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Note that when [beta] [less than or equal to] [square root of 6], we get v* = 0. This implies [[lambda].sup.0,*.sub.1]] = [[lambda].sup.0,*.sub.2]] = [square root of [D.sub.1]] [([[omega].sub.min][[omega].sub.max]).sup.1/4], which corresponds to the zeroth-order one-sided optimal parameters found in [8].

5. Numerical experiments with two subdomains. The model problem (2.2) is discretized using a backward Euler scheme in time and a second-order scheme on a staggered grid in space. For the interior points, the scheme is

[u.sup.n+1.sub.k] - [u.sup.n.sub.k] / [DELTA]t = 1 / [x.sub.k+1/2] - [x.sub.k-1/2] [[F.sup.n+1.sub.k+1/2] - [F.sup.n+1.sub.k-1/2]],

with [F.sup.n+1.sub.k+1/2] = [D.sub.k+1/2] [u.sup.n+1.sub.k] - [u.sup.n.sub.k] /[x.sub.k+1] - [x.sub.k]. Note that for practical applications, the use of the Crank-Nicolson scheme in time is avoided because this leads to unphysical behavior. Indeed, unlike the backward Euler scheme, the Crank-Nicolson scheme fails to satisfy the so-called monotonic damping property [22]. We decompose the computational domain [OMEGA] into two non-overlapping subdomains [[OMEGA].sub.1] = [-[L.sub.1], 0] and [[OMEGA].sub.2] = [0, [L.sub.2]], with [L.sub.1] = [L.sub.2] = 500 m. A homogeneous Neumann boundary condition is imposed at x = -[L.sub.1] and x = [L.sub.2]. As it is usually done in numerical models, the resolution [DELTA][x.sub.k] is progressively refined to enhance the resolution in the boundary layers in the vicinity of the air-sea interface. We use N = 75 points in each subdomain and the resolution varies from [DELTA][x.sub.k] = 25 m at x = [L.sub.1] (respectively x = [L.sub.2]) to [DELTA][x.sub.k] = 1 m at x = 0. The Robin condition [g.sub.N+1/2] on the interface [GAMMA] (located at x = [x.sub.N+1/2] on [[OMEGA].sub.1] and at x = [x.sub.1/2] on [[OMEGA].sub.2]) is discretized by assuming that the flux F is constant on the first cell near [GAMMA]. This leads to

[g.sub.N+1/2] = [D.sub.N-1/2] [u.sub.N] - [u.sub.N-1] / [x.sub.N-1/2] - [x.sub.N-3/2] + [lambda][u.sub.N],

where [lambda] is the Robin parameter. We simulate directly the error equations, i.e., [f.sub.1] = [f.sub.2] = 0 in (2.2) and [u.sub.0] (x) = 0. We start the iteration with a random initial guess [u.sup.0.sub.2] (0, t), t [member of] [0, T], so that it contains a wide range of the temporal frequencies that can be resolved by the computational grid. We perform simulations for four different types of transmission conditions at x = 0: Dirichlet-Neumann (DN), optimized Neumann-Robin (NR*), optimized Robin-Robin (RR*), and asymptotically optimized Robin-Robin ([RR.sup.(as)]). In Figure 5.1 we show the evolution of the [L.sup.[infinity]]-norm of the error obtained for those four cases for [gamma] = [10.sup.1/4] [approximately equal to] 1.7783, [gamma] = [square root of 10] [approximately equal to] 3.1623, and [gamma] = 10, with [mu] = 6 and [mu] = 12. We choose the time steps [DELTA][t.sub.1] = [DELTA][t.sub.2] = [DELTA]t = 100 s, [D.sub.2] = 0.5 [m.sup.2][s.sub.-1], and [D.sub.1] is then deduced depending on the value of [gamma]. As expected, we get the best results with the two-sided Robin conditions. Consistent with Figure 4.5, the convergence is faster when [gamma] is large and when [mu] is small. Moreover, when the discontinuity [gamma] between the diffusion coefficients is increased, the algorithm becomes less and less sensitive to the choice of transmission conditions and to the parameter [mu]. The asymptotic optimized Robin-Robin conditions provide a good approximation of the optimized Robin-Robin conditions, even for [DELTA]t = 100 s >> 0. Those conditions are especially efficient when [gamma] is sufficiently larger than 1. We remark that the optimized Neumann-Robin conditions provide only a slight improvement compared to the classical Dirichlet-Neumann conditions.

Because we consider a problem with discontinuous coefficients, the time step used to advance the diffusion equation may be different in each subdomain. It is thus instructive to look at the impact of non-conformities in time on the performance of the optimized algorithm. For the two cases [gamma] = [10.sup.1/4] and [gamma] = 10, which we considered so far, we adapt the time step in each subdomain so that the ratio r = [DELTA][t.sub.1]/[DELTA][t.sub.2] between the time steps varies from 100 to 1/100. To handle the exchange of boundary data between the non-conforming grids we use a linear method for the interpolation step and an averaging for the restriction step, both steps are conservative. Note that we got very similar results using the [L.sup.2] projection algorithm described in [13, Appendix A]. We consider [[omega].sub.max] = [pi]/ min([DELTA][t.sub.1], [DELTA][t.sub.2]) for the optimization of the Robin conditions.

As shown in Figure 5.2, the performance of the algorithm is degraded as long as r [not equal to] 1. Indeed, the interpolation/restriction step modifies the frequency spectrum of the error and thus affects the convergence speed of the algorithm. For the cases r [not equal to] 1, we have investigated a wide range of values for the parameters [p.sub.j] in the Robin transmission conditions. Even if the algorithm is slower than the one with r = 1, it turned out that the optimal Robin parameters found in Theorem 4.6 are still optimal in the case r [not equal to] 1. Note that the results after one iteration can be quite dependent on the value of r. This is to be expected since the frequency spectrum in the initial error is very different whether the (random) initial guess is initialized on the grid with the smaller time step or the larger time step.

Conclusion. In this paper, we obtain new results for an optimized Schwarz method defined for non-overlapping diffusion problems with discontinuous coefficients. This method uses zeroth-order two-sided Robin transmission conditions, i.e., we consider two different Robin conditions on each side of the interface. We base our approach on a model problem with two subdomains and we prove the convergence of the corresponding algorithm. Then we analytically study the behavior of the convergence factor with respect to the parameters of the problem. We show that the optimized convergence factor satisfies an equioscillation property between two or three points depending on the parameter values. In comparison with other methods using the Neumann-Robin or Dirichlet-Neumann conditions, these two-sided Robin-Robin conditions are significantly more efficient, especially when the ratio between the discontinuous coefficients is close to one. Asymptotic results for [DELTA]t small are given. Numerical results show the performance of the different type of transmission conditions introduced in this paper. Those results are consistent with the analytical study.

Acknowledgements. This research was partially supported by the French National Research Agency (ANR) project "COMMA" and by the INRIA project-team MOISE. The authors would also like to thank two anonymous reviewers, whose comments and suggestions during the review process helped to clarify earlier version of the manuscript.

REFERENCES

[1] D. BENNEQUIN, M. GANDER, AND L. HALPERN, Optimized Schwarz waveform relaxation methods for convection reaction diffusion problems, Technical Report 2004-24, LAGA, Universite Paris 13, 2004.

[2] --, A homographic best approximation problem with application to optimized Schwarz waveform relaxation, Math. Comp., 78 (2009), pp. 185-223.

[3] E. BLAYO, L. HALPERN, AND C. JAPHET, Optimized Schwarz waveform relaxation algorithms with nonconforming time discretization for coupling convection-diffusion problems with discontinuous coefficients, in Domain Decomposition Methods in Science and Engineering XVI, O. Widlund and D. E. Keyes, eds., Lect. Notes Comput. Sci. Eng., 55, Springer, Berlin, 2007, pp. 267-274.

[4] G. Danabasoglu, W. G. Large, J.J. Tribbia, P. R. Gent, B. P. Briegleb, and J. C. McWilliams, Diurnal coupling in the tropical oceans of CCSM3, J. Climate, 19 (2006), pp. 2347-2365.

[5] V. F. DEM'YANOV AND V. N. MALOZEMOV, On the theory of non-linear minimax problems, Russian Math. Surveys, 26 (1971), pp. 57-115.

[6] O. DUBOIS, Optimized Schwarz Methods for the Advection-Diffusion Equation and for Problems with Discontinuous Coefficients, Ph.D. Thesis, Dept. of Math. and Stat., McGill University, Montreal, 2007.

[7] M. J. GANDER, Schwarz methods over the course of time, Electron. Trans. Numer. Anal., 31 (2008), pp. 228-255. http://etna.math.kent.edu/vol.31.2 008/pp228-255.dir

[8] M.J. GANDER AND L. HALPERN, Methodes de relaxation d'ondes pour l''equation de la chaleur en dimension 1, C. R. Acad. Sci. Paris, 336 (2003), pp. 519-524.

[9] --, Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems, SIAM J. Numer. Anal., 45 (2007), pp. 666-697.

[10] M. J. GANDER, L. HALPERN, AND M. KERN, A Schwarz waveform relaxation method for advection-diffusion-reaction problems with discontinuous coefficients and non-matching grids, in Domain Decomposition Methods in Science and Engineering XVI, O. Widlund and D. E. Keyes, eds., Lect. Notes Comput. Sci. Eng., 55, Springer, Berlin, 2007, pp. 283-290.

[11] M. J. GANDER, L. HALPERN, AND F. MAGOULES, An optimized Schwarz method with two-sided Robin transmission conditions for the Helmholtz equation, Internat. J. Numer. Methods Fluids, 55 (2007), pp. 163-175.

[12] M. J. GANDER, L. HALPERN, AND F. NATAF, Optimal convergence for overlapping and non-overlapping Schwarz waveform relaxation, in Eleventh International Conference on Domain Decomposition Methods London 1998, C. H. Lai, P. E. Bjorstad, M. Cross, and O. Widlund, eds., DDM.org, Augsburg, 1999, pp. 27-36.

[13] --, Optimal Schwarz waveform relaxation for the one dimensional wave equation, SIAM J. Numer. Anal., 41 (2003), pp. 1643-1681.

[14] M. J. GANDER AND A. M. STUART, Space-time continuous analysis of waveform relaxation for the heat equation, SIAM J. Sci. Comput., 19 (1998), pp. 2014-2031.

[15] C. JAPHET AND F. NATAF, The best interface conditions for domain decomposition methods: absorbing boundary conditions, in Absorbing Boundaries and Layers, Domain Decomposition Methods, L. Tourrette and L. Halpern, eds., Nova Sci. Publ., Huntington, NY, 2001, pp. 348-373.

[16] E. LELARASMEE, A. RUEHLI, AND A. SANGIOVANNI-VINCENTELLI, The waveform relaxation method for time-domain analysis of large scale integrated circuits, IEEE Trans. Computer-Aided Design Integr. Circuits Syst., 1 (1982), pp. 131-145.

[17] F. LEMARlE, Algorithmes de Schwarz et couplage ocean-atmosphere, Ph.D. Thesis, Laboratoire Jean Kuntzmann, Joseph Fourier University, Grenoble, 2008.

[18] F. LEMARIE, L. DEBREU, AND E. BLAYO, Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients. Part 2: the variable coefficients case, Electron. Trans. Numer. Anal., 40 (2013), pp. 170-186. http://etna.math.kent.edu/vol.40.2 013/pp17 0-18 6.dir

[19] P.-L. LIONS, On the Schwarz alternating method. III. A variant for nonoverlapping subdomains, in Third International Symposium on Domain Decomposition Methods for Partial Differential Equations, T. Chan, R. Glowinski, J. Periaux, and O. Widlund, eds., SIAM, Philadelphia, 1990, pp. 202-223.

[20] Y. MADAY and F. MAGOULES, Non-overlapping additive Schwarz methods tuned to highly heterogeneous media, C. R. Math. Acad. Sci. Paris, 341 (2005), pp. 701-705.

[21] O. S. Madsen, A realistic model of the wind-induced Ekman boundary layer, J. Phys. Oceanogr., 7 (1977), pp. 248-255.

[22] G. Manfredi and M. Ottaviani, Finite-difference schemes for the diffusion equation, in Dynamical Systems, Plasmas and Gravitation, P. Leach, S. Bouquet, J.-L. Rouet, and E. Fijalkow, eds., Lecture Notes in Physics, 518, Springer, Berlin, 1999, pp. 82-92.

[23] J. C. McWilliams, Irreducible imprecision in atmospheric and oceanic simulations, Proc. Natl. Acad. Sci. USA, 104 (2007), pp. 8709-8713.

FLORIAN LEMARIE ([dagger]), LAURENT DEBREU ([dagger]), AND ERIC BLAYO ([double dagger])

* Received September 29, 2010. Accepted March 11, 2013. Published online on July 5, 2013. Recommended by Martin Gander.

([dagger]) INRIA Grenoble Rhone-Alpes, Montbonnot, 38334 Saint Ismier Cedex, France and Jean Kuntzmann Laboratory, BP 53, 38041 Grenoble Cedex 9, France, ({florian.lemarie, laurent.debreu}@inria.fr).

([double dagger]) University of Grenoble, Jean Kuntzmann Laboratory, BP 53, 38041 Grenoble Cedex 9, France (eric.blayo@imag.fr).

Printer friendly Cite/link Email Feedback | |

Author: | Lemarie, Florian; Debreu, Laurent; Blayo, Eric |
---|---|

Publication: | Electronic Transactions on Numerical Analysis |

Article Type: | Report |

Date: | Jan 1, 2013 |

Words: | 11457 |

Previous Article: | A combinatorial approach to nearly uncoupled Markov chains I: reversible Markov chains. |

Next Article: | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients. Part 2: the... |

Topics: |