El Rami

Christian Haddad

Gopalakrishnan Lécot

Venkiteswaran

Haddad

Christian Lécot
email: christian.lecot@univ-savoie.fr

Gopalakrishnan Venkiteswaran

Diffusion in a nonhomogeneous medium: quasi-random walk on a lattice

Keywords: Quasi-Monte Carlo, random walk, low-discrepancy sequences, diffusion equation AMS classification. 11K45, 80M25

We are interested in Monte Carlo (MC) methods for solving the diffusion equation: in the case of a constant diffusion coefficient, the solution is approximated by using particles and in every time step, a constant stepsize is added to or substracted from the coordinates of each particle with equal probability. For a spatially dependent diffusion coefficient, the naive extension of the previous method using a spatially variable stepsize introduces a systematic error: particles migrate in the directions of decreasing diffusivity. A correction of stepsizes and stepping probabilities has recently been proposed and the numerical tests have given satisfactory results. In this paper, we describe a quasi-Monte Carlo (QMC) method for solving the diffusion equation in a spatially nonhomogeneous medium: we replace the random samples in the corrected MC scheme by low-discrepancy point sets. In order to make a proper use of the better uniformity of these point sets, the particles are reordered according to their successive coordinates at each time step. We illustrate the method with numerical examples: in dimensions 1 and 2, we show that the QMC approach leads to improved accuracy when compared with the original MC method using the same number of particles.

Introduction

Random walk techniques are useful computational tools for simulating processes involving the diffusion of substances. The substance is modelled by a set of particles and the particles move according to the dynamics described in the diffusion equation. For diffusion in one dimension and in the case of a constant diffusion coefficient D 0 , there are two main ways of performing random walk. A short time interval ∆t is chosen. With the first method, a particle with coordinate x moves to x + ∆x, where ∆x is a random increment drawn from a centered Gaussian distribution with variance 2D 0 ∆t. This approach is discussed in [START_REF] Chorin | Numerical study of slightly viscous flow[END_REF][START_REF] Ghoniem | Grid-free simulation of diffusion using random walk methods[END_REF][START_REF] Fogelson | Optimal smoothing in function-transport particle methods for diffusion problems[END_REF]. The second method consists of using a constant stepsize ∆x := √ 2D 0 ∆t. A particle with coordinate x jumps to x -∆x or to x + ∆x with equal probability 1/2. This is the approach taken in [START_REF] Kac | Random walk and the theory of Brownian motion[END_REF][START_REF] Feller | An Introduction to Probability Theory and its Applications[END_REF][START_REF] Hald | Convergence of random methods for a reaction-diffusion equation[END_REF]. Both methods are simple to implement and allow to consider complicated geometries. Drawbacks are their slow convergence so that many runs are often needed to obtain reliable results.

When considering a problem with a spatially dependent diffusion coefficient, i.e. D = D(x), the naive extensions of the above methods are to use a spatially variable variance 2D(x)∆t in the first case, and a spatially variable stepsize 2D(x)∆t in the second case. However, both approaches give biased results as they lead to an apparent advection in directions of decreasing diffusivity and a concentration of particles in regions of low diffusivity [START_REF] Hunter | On the use of random walk models with spatially variable diffusivity[END_REF][START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium: correction to the Gaussian steplength[END_REF][START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium: A biased random walk on an asymmetrical lattice[END_REF]. A systematic method for calculating corrections to the Gaussian steplength has been proposed by L. Farnell and W. G. Gibson in [START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium: correction to the Gaussian steplength[END_REF]. The same authors have given in [START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium: A biased random walk on an asymmetrical lattice[END_REF] a method for correcting the fixed step case: both the steplength and the probability of stepping to left and right must be modified. Numerical experiments have shown that both corrected schemes lead to satisfactory results.

One approach to improve the accuracy of MC methods is to change the random numbers used. QMC methods replace pseudo-random numbers with quasi-random numbers. These are points which are very evenly distributed. We refer to [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF] for detailed information on QMC methods. In MC integration, it is not the randomness of the samples that is relevant, but rather that they are spread in a uniform manner on the integration domain, so QMC quadrature often outperforms MC integration . But in particle simulations the independence of the random samples is important. If care is not taken in how the quasi-random points are used, the correlations may destroy the method. An additionnal technique was first introduced in [START_REF] Lécot | A direct simulation Monte Carlo scheme and uniformly distributed sequences for solving the Boltzmann equation[END_REF] and consists of sorting the particles according to position at each time step. This technique was applied to the diffusion equation (in spatially homogeneous media) in [START_REF] Morokoff | A quasi-Monte Carlo approach to particle simulation of the heat equation[END_REF][START_REF] Coulibaly | Simulation of diffusion using quasi-random walk methods[END_REF][START_REF] Lécot | Quasi-Monte Carlo simulation of diffusion[END_REF][START_REF] Venkiteswaran | Quasi-Monte Carlo algorithms for diffusion equations in high dimensions[END_REF]. In a previous communication [START_REF] Haddad | Quasi-Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium[END_REF], we have proposed and tested a QMC algorithm for the simulation of the diffusion equation based on the corrected Gaussian steplength of [START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium: correction to the Gaussian steplength[END_REF]. The aim of the present work is to develop a QMC version of the corrected MC algorithm of [START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium: A biased random walk on an asymmetrical lattice[END_REF].

The rest of the paper is organized as follows. Section 2 presents the classical random walk scheme with constant stepsize and its deterministic version, for the simulation of diffusion in an one-dimensional homogeneous medium. In Section 3, we recall the corrections of the steplength and probability in order to deal with spatially variable diffusion coefficients. Section 4 presents a QMC version of the corrected MC scheme proposed above. In Section 5, we illustrate the new method with numerical examples in dimensions one and two, and we show that the accuracy is improved by the QMC approach, when compared with the standard MC strategy. We draw conclusions in the final section.

Random walk on a lattice

The diffusion equation describes the evolution of the density of a substance due to the displacement of molecules (attributed to Brownian motion) from regions of strong concentration to regions of small concentration. The pure initial value problem for a simple diffusion equation in a one-dimensional domain is:

∂c ∂t (x, t) = ∂ ∂x D ∂c ∂x (x, t), x ∈ R, t > 0, (2.1) c(x, 0) = c 0 (x), x ∈ R, (2.2)
where D := D(x) > 0 is the diffusion coefficient and c 0 (x) ≥ 0 is the initial concentration at position x. Here c(x, t) is the concentration of particles at location x and time t. If we assume

R c 0 (x)dx = 1, (2.3)
then ∀t > 0 R c(x, t)dx = 1. (2.4)
In the spatially homogeneous case, D(x) := D 0 , the random walk technique on a lattice can be described as follows. Let ∆x be a spatial step and ∆t be a time step. The finite difference method replaces Eq. (2.1) with

c(x, t + ∆t) -c(x, t) ∆t ≈ D 0 ∆x 2 (c(x -∆x, t) -2c(x, t) + c(x + ∆x, t)). (2.5)
Hence, if ∆x = √ 2D 0 ∆t, one has c(x, t + ∆t) ≈ 1 2 c(x -∆x, t) + 1 2 c(x + ∆x, t). (2.6)
Given a large integer N , the simulation is conducted by first generating N samples (particles) from the initial probability density function c 0 . Then the particles random walk as follows. During a small time step ∆t, each particle moves by a random distance

L 0 := 2D 0 ∆tB, (2.7)
where B is a Bernoulli random variable such that P(B = -1) = P(B = +1) = 1/2. A QMC version of the previous scheme has been analyzed in [START_REF] Coulibaly | Simulation of diffusion using quasi-random walk methods[END_REF]. The basic idea of QMC methods is to replace the pseudo-random points of MC methods by lowdiscrepancy point sets. We first recall from [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF] some basic concepts. Denote I := [0, 1) and let s be a fixed dimension and λ s be the s-dimensional Lebesgue measure. For a set U = {u 0 , . . . , u N -1 } of points in the s-dimensional unit cube I s , we define the discrepancy by

D N (U ) := sup J 1 N 0≤k<N 1 J (u k ) -λ s (J) ,
where 1 J denotes the indicator function of J, and where the supremum is taken over all subintervals J of I s . We say that a set of N points of I s is a low-discrepancy point set if its discrepancy is of order O((log N ) s-1 /N ), which is the minimum size possible. The theory of (t, m, s)-nets and (t, s)-sequences was developed with a view to providing excellent low-discrepancy point sets for QMC methods. For an integer b ≥ 2, an elementary interval in base b is an interval of the form

s i=1 [a i b -d i , (a i + 1)b -d i ), with integers d i ≥ 0 and 0 ≤ a i < b d i for 1 ≤ i ≤ s. If 0 ≤ t ≤ m are integers, a (t, m, s)-net in base b is a point set consisting of b m points in I s such that every elementary interval J in base b of measure b t-m contains exactly b t points of the set. If b ≥ 2 and t ≥ 0 are integers, a sequence {u 0 , u 1 , . . .} of points in I s is called a (t, s)-sequence in base b if, for all integers n ≥ 0 and m > t, the points u j with nb m ≤ j < (n + 1)b m form a (t, m, s)-net in base b.
We sketch the quasi-random walk method for the simulation of (2.1)-(2.2) in the case of a constant diffusion coefficient D 0 . The first step towards the approximation of the solution entails the discretization of the space variable x only. Let ∆x be a space step and τ ±∆x be the operator defined by

τ ±∆x f (x) := f (x ± ∆x), for any function f on R. The semi-discrete approximation to (2.1)-(2.2) is ∂ c ∂t (x, t) = D 0 ∆x 2 τ -∆x c -2 c + τ ∆x c (x, t), x ∈ R, t > 0, (2.8) c(x, 0) = c 0 (x), x ∈ R. (2.9)
It is easy to show that the conservation law still holds:

∀t > 0 R c(x, t)dx = 1. (2.10)
We first construct a weak formulation of the equation. A mapping σ : R → R is said to be simple if σ(R) is a finite set. Let S + (R) be the set of all nonnegative simple measurable functions on R. Then, for any σ ∈ S + (R),

d dt R σ(x) c(x, t)dx = D 0 ∆x 2 R τ -∆x σ -2σ + τ ∆x σ (x) c(x, t)dx. (2.11)
Let b and m be integers and denote N := b m . We generate N samples (particles)

x 0 k , 0 ≤ k < N from the initial probabillity distribution c 0 (x)dx, i.e. c 0 (x) := 1 N 0≤k<N δ(x -x 0 k ) ≈ c 0 (x)dx,
where δ(x-x * ) denotes the Dirac measure at x * . Next we discretize time into intervals of length ∆t with the following stability requirement

λ := D 0 ∆t ∆x 2 ≤ 1 2 .
(2.12)

For QMC integration, we need a low discrepancy sequence:

U = {u 0 , u 1 , . . .} ⊂ I 2 .
For an integer n, denote U n := {u j : nN ≤ j < (n + 1)N }. If Π 1 is the map defined by Π 1 (u 1 , u 2 ) := u 1 , we assume that for all n ∈ N,

Π 1 U n is a (0, m, 1)-net in base b. (2.13)
Assuming that we have computed an approximation c n of c at time t n := n∆t by a sum of Dirac masses located at positions x n 0 , . . . , x n N -1 :

c n (x) := 1 N 0≤k<N δ(x -x n k ),
we compute an approximation c n+1 at time t n+1 in two steps.

(1) Sorting the particles. The particles are relabeled such that

x n 0 ≤ x n 1 ≤ • • • ≤ x n N -1 . (2.14)
This reordering was first introduced in [START_REF] Lécot | A direct simulation Monte Carlo scheme and uniformly distributed sequences for solving the Boltzmann equation[END_REF] for the simulation of the Boltzmann equation. It guarantees the convergence of the method: since a step of the random walk can be described by a numerical integration (see below), the sorting reverts to minimizing the amplitude of the jumps of the function to be integrated.

(2) Quasi-random walk. A measure γ n+1 (x) which approximates the solution at time t n+1 is obtained by using a forward Euler scheme to discretize Eq. (2.11): for any σ ∈ S + (R),

1 ∆t R σ(x)(γ n+1 (x) -c n (x)) = D 0 ∆x 2 R τ -∆x σ -2σ + τ ∆x σ (x)c n (x). (2.15) Consequently, R σ(x)γ n+1 (x) = 1 N 0≤k<N (λσ(x n k -∆x)+(1-2λ)σ(x n k )+λσ(x n k +∆x)). (2.16)
Denote by 1 k the indicator function of the elementary interval

I k := [k/N, (k +1)/N )
and let χ -1 , χ 0 and χ +1 denote, respectively, the indicator functions of the intervals

I -:= [0, λ), I 0 := [λ, 1 -λ), I + := [1 -λ, 1). To σ ∈ S + (R), we associate Σ n (u) := 0≤k<N 1 k (u 1 ) ε=-1,0,+1 χ ε (u 2 )σ(x n k + ε∆x), u = (u 1 , u 2 ) ∈ I 2 . Then R σ(x)γ n+1 (x) = I 2
Σ n (u)du.

(2.17)

We recover c n+1 by performing a QMC integration:

R σ(x)c n+1 (x) = 1 N nN ≤j<(n+1)N Σ (n) (u j ). (2.18)
This can be reworded as follows. For u ∈ [0, 1), put k(u) := N u , where ξ denotes the greatest integer ≤ ξ. It follows from (2.13) that the mapping

j ∈ {nN, nN + 1, . . . , (n + 1)N -1} → k(u j,1 ) ∈ {0, 1, . . . , N -1} (2.19)
is a bijection. Then, the random walk restatement of (2.18) reads: for nN ≤ j < (n + 1)N ,

x n+1 k(u j,1 ) =        x n k(u j,1 ) -∆x, if 0 ≤ u j,2 < λ, x n k(u j,1 ) , if λ ≤ u j,2 < 1 -λ, x n k(u j,1 ) + ∆x, if 1 -λ ≤ u j,2 < 1.
(2.20)

If λ = 1/2, this reduces to:

x n+1 k(u j,1 ) = x n k(u j,1 ) -∆x, if 0 ≤ u j,2 < 1/2, x n k(u j,1 ) + ∆x, otherwise, (2.21)
which must be compared with (2.7). For every j, the first coordinate u j,1 of the quasirandom point u j is used to select the particle, while the second coordinate u j,2 is used to decide whether the particle moves left or right. This approach may be extended to several dimensions. In [START_REF] Coulibaly | Simulation of diffusion using quasi-random walk methods[END_REF], convergence of the algorithm is established in the multidimensional case; the results of computational experiments indicate that a significant improvement is achieved over standard random walk simulation.

In some applications, the diffusion coefficient D(x) varies with position. In this case, the obvious extension of the previous method using a spatially variable space step 2D(x)∆t leads to the accumulation of particles in regions of low diffusivity. In [START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium: A biased random walk on an asymmetrical lattice[END_REF], L. Farnell and W. G Gibson have derived simple formulas for correcting the standard fixed-step MC method in the case of a spatially variable diffusion coefficient. We recall their MC algorithm in the next section.

A biased random walk in a nonhomogeneous medium

We suppose that the diffusion coefficient D(x) > 0 varies with distance. Let ∆t be the time step and L(x) := 2D(x)∆t denote the uncorrected steplength. During the time

∆ m (x) = 1 2 (L(x -∆ m (x)) + L(x)), (3.1)
∆ p (x) = 1 2 (L(x) + L(x + ∆ p (x))), (3.2)
the random walk stays on a (irregular) lattice: see Fig.

∆ m (x) = L(x) - 1 2 L(x)L (x) + O((∆t) 3/2 ), (3.3)
∆ p (x) = L(x) + 1 2 L(x)L (x) + O((∆t) 3/2 ). (3.4)
Denote by T (x) (resp. T r (x)) the probability that a particle at x moves to the left (resp. to the right) during ∆t: clearly T (x) + T r (x) = 1. If we consider a large number of particles, the average number N (x, t) of particles at point x and time t satisfies:

N (x, t + ∆t) = N (x -∆ m , t)T r (x -∆ m ) + N (x + ∆ p , t)T (x + ∆ p ), (3.5)
where ∆ m := ∆ m (x) and ∆ p := ∆ p (x). The concentration of particles at point x and time t is c(x, t) = N (x, t)/L(x) and the time rate of change of concentration at point x verifies:

∂c ∂t (x, t) = lim ∆t→0 1 L(x) N (x, t + ∆t) -N (x, t)) ∆t . (3.6)
It follows then from (3.5):

∂c ∂t (x, t) = lim ∆t→0 1 L(x) 1 ∆t (L(x + ∆ p )c(x + ∆ p , t)T (x + ∆ p ) +L(x -∆ m )c(x -∆ m , t)T r (x -∆ m ) -L(x)c(x, t)). (3.7)
Using the Taylor expansions

c(x -∆ m , t) = c(x, t) -∆ m ∂c ∂x (x, t) + 1 2 ∆ 2 m ∂ 2 c ∂x 2 (x, t) + O((∆) 3/2 ), (3.8) c(x + ∆ p , t) = c(x, t) + ∆ p ∂c ∂x (x, t) + 1 2 ∆ 2 p ∂ 2 c ∂x 2 (x, t) + O((∆) 3/2 ) (3.9)
in Eq. (3.7), we get

∂c ∂t (x, t) = lim ∆t→0 1 ∆t d 0 c(x, t) + d 1 ∂c ∂x (x, t) + d 2 ∂ 2 c ∂x 2 (x, t) , (3.10)
where These conditions are satisfied if we have

d 0 := R p (x)T (x + ∆ p ) + R m (x)T r (x -∆ m ) -1, d 1 := R p (x)∆ p T (x + ∆ p ) -R m (x)∆ m T r (x -∆ m ), d 2 := 1 2 R p (x)∆ 2 p T (x + ∆ p ) + R m (x)∆ 2 m T r (x -∆ m ) ,
T (x) = 1 2 - 1 4 L (x) + O(∆t) and T r (x) = 1 2 + 1 4 L (x) + O(∆t). (3.12)
Details of the calculation are given in the original paper [START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium: A biased random walk on an asymmetrical lattice[END_REF].

The corrected random walk can be summarized as follows: a particle at position x at time t should at time t + ∆t, with probability T (x) move to position x -∆ m (x) and with probability T r (x) move to position x + ∆ p (x). Here

∆ m (x) := L(x) - 1 2 L(x)L (x), T (x) := 1 2 - 1 4 L (x), (3.13)
∆ p (x) := L(x) + 1 2 L(x)L (x), T r (x) := 1 2 + 1 4 L (x) (3.14)
and L(x) := 2D(x)∆t is the uncorrected steplength.

A quasi-random walk on an asymmetrical lattice

In this section, we develop a QMC approach for the simulation of the diffusion equation in a multi-dimensional nonhomogeneous medium. The algorithm is based on the QMC scheme described in Section 2 and uses the correction method of Section 3. We consider the pure initial value problem:

∂c ∂t (x, t) = ∇ • (D∇c)(x, t), x ∈ R s , t > 0, (4.1) c(x, 0) = c 0 (x), x ∈ R s , (4.2)
where the initial concentration c 0 (x) ≥ 0 satisfies

R s c 0 (x)dx = 1. (4.3) This implies ∀t > 0 R s c(x, t)dx = 1. (4.4)
Let ∆t be a time step and t n := n∆t, for n ∈ N.

Constant number of particles

We first use a fixed number of particles. Let b and d 1 , . . . , d s be integers. Put m :=

d 1 + • • • + d s and N = b m .
We are looking for an approximation of the solution at time t n of the form

c n := 1 N 0≤k<N δ(x -x n k ),
where x n 0 , . . . , x n N -1 are the (positions of the) particles at time t n . The QMC algorithm has several steps:

Initialization

The initial locations {x 0 0 , . . . , x 0 N -1 } are sampled from the probability distribution c 0 (x)dx. This is done by mapping a (0, m, s)-net in base b to R s . The choice of the mapping depends on the initial distribution to be sampled; in our numerical experiments (section 5), the particles are released from a point source at origin, so we take

x 0 0 = • • • = x 0 N -1 := 0.

Displacement of particles

We use a low-discrepancy sequence U = {u 0 , u 1 , . . .} ⊂ I 2s for moving the particles. For an integer n, denote U n := {u j : nN ≤ j < (n + 1)N }. Let Π : I 2s → I s be the projection defined by Π (u 1 , . . . , u 2s ) := (u 1 , . . . , u s ); we assume that for all n ∈ N,

Π U n is a (0, m, s)-net in base b. (4.5)
For u ∈ I 2s , we put u := Π (u). We generalize formulas (3.13)-(3.14) as follows. If L(x) := 2D(x)∆t, we define for 1 ≤ i ≤ s:

∆ i m (x) := L(x) - 1 2 L(x) ∂L ∂x i (x), T i (x) := 1 2 - 1 4
∂L ∂x i (x), (4.6)
∆ i p (x) := L(x) + 1 2 L(x) ∂L ∂x i (x), T i r (x) := 1 2 + 1 4 ∂L ∂x i (x). (4.7)
Suppose that at time t n , the particles are located at positions x n 0 , . . . , x n N -1 : we determine the positions at time t n+1 in two steps.

b = 2, s = 2, d 1 = 2, d 2 = 3)
Relabeling the particles. The particles are labeled x (n) a using a multi-dimensional index a = (a 1 , . . . , a s ) with 0 ≤ a i < b d i for 1 ≤ i ≤ s, such that if a and a are distinct indices, then

a 1 = a 1 , . . . , a i-1 = a i-1 , a i < a i ⇒ x n a,i ≤ x n e a,i . (4.8)
This type of sorting was first used in [START_REF] Lécot | A quasi-Monte Carlo scheme using nets for a linear Boltzmann equation[END_REF] and generalizes the classical one-dimensional sorting of (2.14). This is done in s steps as follows: the values of the first coordinate x 1 of the particles are used to sort them into b d 1 subsets of level 1. In each subset of level 1, the values of the second coordinate x 2 of the particles are used to sort them into b d 2 subsets of level 2, and so on: see Fig. 2 for an example. x n a(u j ),i -∆ i m (x n a(u j ) ), if u j,s+i < T i (x n a(u j ) ), x n a(u j ),i + ∆ i p (x n a(u j ) ), otherwise.

Performing a random walk

(4.10)

The first s coordinates of the quasi-random point u j are used to select the particle, while the last s coordinates are used to update the coordinates of the particle.

Varying number of particles

The QMC algorithm of [START_REF] Coulibaly | Simulation of diffusion using quasi-random walk methods[END_REF] works with a constant number N of numerical particles and N is equal to a power of some prime base b. This is required if we want the mapping (2.19) to be a bijection. In an application of interest to us, physical particles are released from a source so that the number of particles increases over time. We can overcome the difficulty as follows. If N is the number of particles at time t n , we write the digit expansion of N in base b:

N = ∞ h=0 α h (N )b h , (4.11)
where 0 ≤ α h (N ) < b and α h (N ) = 0 for sufficiently large h. The set of N particles is split into α 0 (N ) subsets of size b 0 , α 1 (N ) subsets of size b 1 , and so forth. Each subset is treated separately (for relabeling and displacement). At the end of the time step all the subsets are merged in one set of N particles before adding new particles that are released from the source during the time step.

Numerical experiments

In this section, we present the results of numerical experiments which illustrate the precision of the QMC scheme. We focus on biological applications with Ca 2+ diffusion inside the neuromuscular junction of the crayfish [START_REF] Matveev | New and corrected simulations of synaptic facilitation[END_REF]. Calcium ions are released from a point source as the result of a train of action potentials arriving at that location, with each impulse resulting in release. The ions then diffuse independently through the surrounding medium and act on receptors at varying distances from the release point, chosen as the origin of the axes. The diffusion coefficient varies with position.

The number of particles released from the source is N = 2 14 . We have considered the following cases in dimensions 1 and 2:

(i) instantaneous emission of particles at t = 0;

(ii) emission of particles during a time interval [0, τ e ].

We evaluate the concentration of particles in a region of space as the number of particles in the region divided by the total number of particles. We compute the evolution of concentration in different regions during a time interval [0, T f ]. The time step is ∆t = 1.6 10 -4 . In both cases τ e := 1.31072 = 2 13 ∆t and T f := 10 = 62 500 ∆t.

We compare the results given by the QMC scheme with those produced by the MC scheme of [START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium: A biased random walk on an asymmetrical lattice[END_REF]. Since no analytical solution is available, a reference solution is obtained by using a MC simulation with a very large number of simulation particles.

Non-instantaneous emission

During τ e , two particles are emitted from the source at every time step ∆t. The results of MC and QMC simulations, as plotted on Fig. 6, show the superiority of QMC approach.

Conclusion

We have presented a QMC algorithm for the simulation of diffusion processes in spatially nonhomogeneous media. The scheme is a generalization of the QMC method described in [START_REF] Coulibaly | Simulation of diffusion using quasi-random walk methods[END_REF] in the case of a constant diffusion coefficient. The diffusing substance is approximated by a set of N particles, which are sampled at t = 0 from the initial distribution. Time is discretized, and the particles are reordered according to their coordinates at the beginning of every time step. A mapping between a low-discrepancy sequence and the set of the positions of the particles determines the movement of each of them. The steps of the random walk are corrected according to formulas derived in [START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium: A biased random walk on an asymmetrical lattice[END_REF]. The number N is allowed to change during the simulation. The results of some numerical examples have shown that the new QMC method is more accurate than the standard MC scheme. Directions for future work include the study of theoretical convergence of the method and its application to more general diffusion problems with boundary conditions.

Figure 1 .

 1 Figure 1. The lattice used for the random walk.

  1. Since L(x) = O( √ ∆t), an approximate solution of (3.1)-(3.2) can be found from Taylor expansions:

  together with R m (x) := L(x-∆ m )/L(x) and R p (x) := L(x+∆ p )/L(x). We recover the diffusion equation(2.

Figure 2 .

 2 Figure 2. Sorting 2 2+3 = 32 particles (b = 2, s = 2, d 1 = 2, d 2 = 3)

  We use the quasi-random points u j . To each u = (u 1 , . . . , u s ) ∈ I s , we associate the index a(u) := ( b d 1 u 1 , . . . , b ds u s ). It follows from (4.5) that the mapping j ∈ {nN, nN + 1, . . . , (n + 1)N -1} → a(u j ) (4.9) is a bijection onto {0, . . . , b d 1 -1}ו • •×{0, . . . , b ds -1}.At time t n+1 , the coordinates of the particles are updated as follows: for nN ≤ j < (n + 1)N and 1 ≤ i ≤ s,x n+1 a(u j ),i =   

Figure 3 .

 3 Figure 3. Instantaneous emission in one-dimensional medium. Time evolution of the concentration in the intervals [20, 40] (up) and [60, 80] (down). Comparison of MC (left) and QMC (right) simulations with a reference solution.

Figure 4 .

 4 Figure 4. Non-instantaneous emission in one-dimensional medium. Time evolution of the concentration in the intervals [20, 40] (up) and [60, 80] (down). Comparison of MC (left) and QMC (right) simulations with a reference solution.

Figure 5 .

 5 Figure 5. Instantaneous emission in two-dimensional medium. Time evolution of the concentration in the squares [20, 40] × [10, 30] (up) and [60, 80] × [10, 30] (down). Comparison of MC (left) and QMC (right) simulations with a reference solution.

Figure 6 .

 6 Figure 6. Non-instantaneous emission in two-dimensional medium. Time evolution of the concentration in the squares [20, 40] × [10, 30] (up) and [60, 80] × [10, 30] (down). Comparison of MC (left) and QMC (right) simulations with a reference solution.

One-dimensional diffusion

We consider the spatially varying diffusion coefficient used in [START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium: correction to the Gaussian steplength[END_REF][START_REF] Farnell | Monte Carlo simulation of diffusion in a spatially nonhomogeneous medium: A biased random walk on an asymmetrical lattice[END_REF]: D(x) := D(1 -0.8u(x)), (5.1) where D := 4 000 and

The quasi-random sequence used in 1-D simulations is the (0, 2)-sequence of Faure in base b = 2 (see [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF]). We compare MC and QMC results by calculating the time evolution of the concentration of ions in the space intervals [20, 40] and [60, 80]. To smooth the curves and to make comparison clearer, we average the results over a time interval of amplitude 0.16. The reference solution is obtained by a MC simulation using a sample 2 6 times larger than the original one, i.e., with 2 20 particles.

Instantaneous emission

All the particles are released at time t = 0 from the point source. The concentrations are displayed in Fig. 3. Using quasi-random numbers in place of pseudo-random numbers and reordering the particles clearly reduce scattering of the results.

Non-instantaneous emission

Now the ions are released over the time interval [0, τ e ]: two particles are emitted at each time step ∆t. Figure 4 shows the outputs of the computations. Once again, the quasi-random strategy produces more accurate approximations than the standard random walk method, and the splitting technique used during emission seems to be efficient.

Two-dimensional diffusion

The diffusion coefficient is that of Eq. (5.1), with x replaced by the radial coordinate

Here, the low-discrepancy sequence used is the (3, 4)-sequence of Niederreiter in base b = 2 (see [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF]). We compare MC and QMC results by calculating the time evolution of the concentration of ions in the squares [20, 40] × [10, 30] and [60, 80] × [START_REF] Hunter | On the use of random walk models with spatially variable diffusivity[END_REF]30]. We average the results over a time interval of amplitude 0.24. The reference solution is obtained by a MC simulation using a sample 2 7 times larger than the original one, i.e., with 2 21 particles.

Instantaneous emission

At t = 0, all the particles are released from the point source. The concentrations are displayed in Fig. 5. The improvement due to QMC is clear even in dimension 2.


Lori Badea

Marius Cocou
email: cocou@lma.cnrs-mrs.fr

Approximation results and subspace correction algorithms for implicit variational inequalities

Keywords: Evolution problems, implicit variational inequalities, approximation and existence results, subspace correction, domain decomposition, Schwarz method, multilevel methods

This paper deals with the mathematical analysis and the subspace approximation of a system of variational inequalities representing a unified approach to several quasistatic contact problems in elasticity. Using an implicit time discretization scheme and some estimates, convergence properties of the incremental solutions and existence results are presented for a class of abstract implicit evolution variational inequalities involving a nonlinear operator. To solve the corresponding semi-discrete and the fully discrete problems, some general subspace correction algorithms are proposed, for which global convergence is analyzed and error estimates are established.

Introduction

This work is concerned with the mathematical analysis and the approximation of a system of evolution variational inequalities representing a unified approach to several quasistatic contact problems in elasticity.

The general results presented here constitute a generalization of the cases studied in [START_REF] Capatina | A class of implicit variational inequalities and applications to frictional contact[END_REF] and [START_REF] Cocou | Formulation and approximation of quasistatic frictional contact[END_REF] and can be applied to various quasistatic contact problems, including unilateral contact with nonlocal friction, normal compliance conditions with friction or more complex interaction laws, as, for example, interface laws coupling unilateral contact, adhesion and nonlocal friction between two elastic bodies [START_REF] Raous | A consistent model coupling adhesion, friction, and unilateral contact[END_REF].

Using an implicit time discretization scheme and some estimates of the incremental solutions, approximation and existence results are presented for a class of abstract implicit evolution variational inequalities involving a nonlinear operator.

To solve the general elliptic quasi-variational inequalities of the second kind that are obtained by the previous incremental procedure, some subspace correction algorithms are proposed, for which global convergence is analyzed and error estimates are established.

If the subspaces are the finite element spaces of the fine grid associated with a decomposition of the domain, or with the space corresponding to the coarse mesh, these algorithms are in fact one-and two-level Schwarz methods. In this case, we are able to write the convergence rate depending on the overlapping and mesh parameters. Following this way, we can show that our methods have an optimal convergence rate, i.e. their convergence is the same as in the case of linear equations.

Schwarz methods are widely applied for solving linear problems, because they provide robust and efficient solution methods but their generalization to non-linear problems as, for example, quasi-variational inequalities, is not straightforward. In particular, gaining an estimate for the convergence speed of a two-level or multilevel Schwarz method in the case of non-linear problems is far from trivial.

The methods we deal here generalize the iterative method suggested in [START_REF] Mandel | A multilevel iterative method for symmetric, positive definite linear complementarity problems[END_REF] and [START_REF] Mandel | Etude algébrique d'une méthode multigrille pour quelques problèmes de frontière libre[END_REF] for complementarity problems. For these problems, this projected multilevel relaxation was later developed in [START_REF] Kornhuber | Monotone multigrid methods for elliptic variational inequalities I[END_REF]- [START_REF] Kornhuber | Adaptive Monotone Multigrid Methods for Nonlinear Variational Problems[END_REF] and named as monotone multigrid method. On the other hand, the application of this method to other types of convex sets in general abstract spaces and monotone minimizing functionals have been investigated in [START_REF] Badea | Convergence rate of a multiplicative Schwarz method for strongly nonlinear inequalities[END_REF] and [START_REF] Badea | Convergence rate of a Schwarz multilevel method for the constrained minimization of non-quadratic functionals[END_REF], for instance. Also, the case where the inequality contains extra terms which do not stem from the minimization of a functional has been investigated in [START_REF] Badea | Schwarz methods for inequalities with contraction operators[END_REF]. Additional non-linear terms have also to be considered in the case of quasi-variational, or implicit, inequalities.

Let us finally emphasize that in our opinion, the conditions for a global and optimal convergence rate of the methods with more than two levels (multilevel or multigrid methods) are not yet well understood and that their theoretical understanding will likely require a different approach.

The paper is organized as follows. In Section 2, using an implicit time discretization scheme and some estimates, convergence properties are given and an existence result is established.

In Section 3, general convergence results, based on an semi-discrete internal approximation, are presented.

In Section 4, some subspace correction algorithms are introduced for the solution of the problem discretized in time. These algorithms represent oneand two-level Schwarz methods for the finite element form of the problem. First, three algorithms are introduced in a general framework in a Hilbert space. We introduce here an assumption on the convex set and the correction subspaces, which will be useful in the proof of the convergence of the algorithms. Mainly, this hypothesis refers to the decomposition of the elements in the convex set, and introduces a constant C 0 which will play an important role in the writing of the convergence rate. Another hypothesis is made on the non differentiable term in the inequality. Under these assumptions, we prove that the three subspace algorithms are convergent and give an estimation of the convergence rate. The reminder of the section is devoted to the one-and two-level methods. If we associate the correction subspaces to a domain decomposition, the abstract algorithms are Schwarz methods. We show that the assumptions introduced in the general theory hold and explicitly write the constant C 0 depending on the mesh and domain decomposition parameters. In this way, we get that the convergence rates of the one-and two-level methods for our quasi-variational inequalities are similar with the convergence rates obtained for equations, i.e., we get an optimal convergence. In the case of the two-level methods, the convergence rate is almost independent of the mesh and domain decomposition parameters.

A system of implicit evolution inequalities

Let (V, . , . ), (H, (•, •) H ) be two real Hilbert spaces with the associated norms

• and • H , respectively. Let K be a non empty closed convex cone contained in V with its vertex at 0 and (K(g)) g∈V be a family of nonempty convex subsets of K satisfying the following conditions: 0 ∈ K(0) and

if g n → g in V, v n ∈ K(g n ) and v n v in V then v ∈ K(g). (1)
We consider a functional F : V → R Gateaux differentiable on V and assume that there exist two constants α, β > 0 for which

α v -u 2 ≤ F (v) -F (u), v -u (2)
and

F (v) -F (u) V ≤ β v -u (3)
for all u, v ∈ V , where F is the Gateaux derivative of F .

Using the relations

F (v) -F (u) = 1 0 F (u + r(v -u)), v -u dr = F (u), v -u + 1 0 F (u + r(v -u)) -F (u), v -u dr
and ( 2), (3), it is easily seen that for all u, v ∈ V we have

F (u), v -u + α 2 v -u 2 ≤ F (v) -F (u) ≤ F (u), v -u + β 2 v -u 2 . ( 4
)
Since F satisfies (4), it follows that F is strictly convex, sequentially weakly lower semicontinuous and differentiable on V . We assume that for all g ∈ V there exists an operator γ(g,

•) : K(g) → H such that γ(0, 0) = 0, if g n → g in V, v n ∈ K(g n ) and v n v in V then γ(g n , v n ) → γ(g, v) in H (5)
and for all

g i ∈ V , v i ∈ K(g i ), i = 1, 2, γ(g 1 , v 1 ) -γ(g 2 , v 2 ) H ≤ k 1 ( g 1 -g 2 + v 1 -v 2 ), (6)
where k 1 is a positive constant. For all g ∈ V , let j(g, •, •) : K(g) × V → R be a functional satisfying the following conditions:

j(g, v, •) is sequentially weakly continuous on V ∀ g ∈ V, v ∈ K(g), (7) j(g, v, •) is sub-additive for all g ∈ V, v ∈ K(g), that is (8) j(g, v, w 1 + w 2 ) ≤ j(g, v, w 1 ) + j(g, v, w 2 ) ∀ g, w 1,2 ∈ V, v ∈ K(g), j(g, v, •) is positively homogeneous for all g ∈ V, v ∈ K(g), (9)
that is j(g, v, θw) = θj(g, v, w) ∀ g, w ∈ V, v ∈ K(g), θ ≥ 0,

j(0, 0, w) = 0 ∀ w ∈ V, (10)
and there exists k 2 > 0 such that

|j(g 1 , v 1 , w 2 ) + j(g 2 , v 2 , w 1 ) -j(g 1 , v 1 , w 1 ) -j(g 2 , v 2 , w 2 )| ≤k 2 ( g 1 -g 2 + γ(g 1 , v 1 ) -γ(g 2 , v 2 ) H ) w 1 -w 2 ∀ g i , w i ∈ V, v i ∈ K(g i ), i = 1, 2. (11)
We assume that k 1 and k 2 satisfy the following property:

k 1 k 2 < α. (12)
For all g ∈ V , we consider a functional b(g,

•, •) : K(g) × V → R such that ∀ g ∈ V, v ∈ K(g), b(g, v,
•) is linear and continuous on V

and

|b(g 1 , v 1 , w) -b(g 2 , v 2 , w)| ≤ k b ( g 1 -g 2 + v 1 -v 2 ) w ∀ g i , w ∈ V, v i ∈ K(g i ), i = 1, 2, (14)
where k b is a positive constant. From the above properties of F, j and K and by a classical argument, it follows that for all g ∈ V, d ∈ K, w ∈ K(g) the elliptic variational inequality

u ∈ K F (u), v -u + j(g, w, v -d) -j(g, w, u -d) ≥ 0 ∀ v ∈ K
has a unique solution, so that we can define the mapping S g,d : K(g) → K by S g,d (w) = u. We assume that for all g ∈ V, d ∈ K K(g) is stable under S g,d i.e. S g,d (K(g)) ⊂ K(g).

For all g ∈ V, d ∈ K, we consider the following problems:

( P )

     u ∈ K(g) F (u), v -u + j(g, u, v -d) -j(g, u, u -d) ≥ b(g, u, v -u) ∀ v ∈ V, b(g, u, z -u) ≥ 0 ∀ z ∈ K, ( Q) u ∈ K(g) F (u), v -u + j(g, u, v -d) -j(g, u, u -d) ≥ 0 ∀ v ∈ K,
and we assume that if u is a solution of ( Q), then u is a solution of ( P ). ( 16)

Remark 1. If u satisfies ( P ), then u obviously satisfies ( Q).

Let f ∈ W 1,2 (0, T ; V ) be given. Using the hypotheses (2), ( 3), ( 15), ( 11) and [START_REF] Capatina | Internal approximation of quasi-variational inequalities[END_REF], it follows that S f (0),0 : K(f (0)) → K(f (0)) is a contraction if the condition (12) holds. Thus, the following implicit elliptic variational inequality has a unique solution u 0 ∈ K(f (0)) (see, e.g. [START_REF] Capatina | Internal approximation of quasi-variational inequalities[END_REF]):

F (u 0 ), w -u 0 + j(f (0), u 0 , w) -j(f (0), u 0 , u 0 ) ≥ 0 ∀ w ∈ K. ( 17
)
We consider the following evolution problems, involving implicit variational inequalities.

Problem P: Find u ∈ W 1,2 (0, T ; V ) such that (P )

         u(0) = u 0 , u(t) ∈ K(f (t)) ∀ t ∈]0, T [, F (u(t)), v -u(t) + j(f (t), u(t), v) -j(f (t), u(t), u(t)) ≥ b(f (t), u(t), v -u(t)) ∀ v ∈ V a.e. on ]0, T [, b(f (t), u(t), w -u(t)) ≥ 0 ∀ w ∈ K, ∀ t ∈]0, T [, Problem Q: Find u ∈ W 1,2 (0, T ; V ) such that (Q)      u(0) = u 0 , u(t) ∈ K(f (t)) ∀ t ∈]0, T [, F (u(t)), w -u(t) + j(f (t), u(t), w -u(t) + u(t)) -j(f (t), u(t), u(t)) ≥ 0 ∀ w ∈ K a.e. on ]0, T [, and
Problem Q: Find u ∈ W 1,2 (0, T ; V ) such that ( Q)        u(0) = u 0 , u(t) ∈ K(f (t)) ∀ t ∈]0, T [, F (w) -F (u(t)) + j(f (t), u(t), w -u(t) + u(t)) -j(f (t), u(t), u(t)) ≥ α 2 w -u(t) 2 ∀ w ∈ K a.e. on ]0, T [.
Remark 2. i) By ( 16) and Remark 1 with d = uu, the problems P and Q are equivalent. ii) Since F satisfies (4) and using the definition of the Gateaux derivative of F with the convexity of j(f, u, •), it follows that the problems Q and Q are equivalent.

We shall prove the existence of a solution to problem P by using an implicit time discretization scheme and its convergence properties.

For ν ∈ N * , we set ∆t := T /ν, t ι := ι ∆t and

K ι := K(f (t ι )), ι = 0, 1, ..., ν. If θ is a continuous function of t ∈ [0, T ]
valued in some vector space, we use the notations θ ι := θ(t ι ) unless θ = u, and if ζ ι , ∀ ι ∈ {0, 1, ..., ν}, are elements of some vector space, then we set

∂ζ ι := ζ ι+1 -ζ ι ∆t ∀ ι ∈ {0, 1, ..., ν -1}.
We denote u 0 := u 0 and we approximate (P ) using the following sequence of incremental problems (P ι ν ) ι=0,1,...,ν-1 . Problem P ι ν : Find u ι+1 ∈ K ι+1 such that

(P ι ν )      F (u ι+1 ), v -∂u ι + j(f ι+1 , u ι+1 , v) -j(f ι+1 , u ι+1 , ∂u ι ) ≥ b(f ι+1 , u ι+1 , v -∂u ι ) ∀ v ∈ V, b(f ι+1 , u ι+1 , w -u ι+1 ) ≥ 0 ∀ w ∈ K.
By [START_REF] Mandel | Etude algébrique d'une méthode multigrille pour quelques problèmes de frontière libre[END_REF] and Remark 1 for g = f ι+1 , d = u ι , and using similar arguments as in Remark 2 ii), respectively, it is easily seen that for all ι ∈ {0, 1, ..., ν -1} the problem P ι ν is equivalent to each of the following variational inequalities:

find u ι+1 ∈ K ι+1 such that (Q ι ν ) F (u ι+1 ), w -u ι+1 + j(f ι+1 , u ι+1 , w -u ι ) -j(f ι+1 , u ι+1 , u ι+1 -u ι ) ≥ 0 ∀ w ∈ K, find u ι+1 ∈ K ι+1 such that ( Qι ν )    F (w) -F (u ι+1 ) + j(f ι+1 , u ι+1 , w -u ι ) -j(f ι+1 , u ι+1 , u ι+1 -u ι ) ≥ α 2 w -u ι+1 2 ∀ w ∈ K.
From the hypotheses (2), ( 12), ( 15), ( 11) and ( 6), it follows that S f ι+1 ,u ι : K ι+1 → K ι+1 is a contraction. Therefore (Q ι ν ) has a unique solution which is equally the unique solution of (P ι ν ) and of ( Qι ν ), for all ι ∈ {0, 1, ..., ν -1}. Let us define the following functions:

               u ν (0) = ûν (0) = u 0 , f ν (0) = f 0 and ∀ ι ∈ {0, 1, ..., ν -1}, ∀ t ∈]t ι , t ι+1 ], u ν (t) = u ι+1 , ûν (t) = u ι + (t -t ι )∂u ι , f ν (t) = f ι+1 .
Then for all ν ∈ N * the sequence of inequalities (P ι ν ) ι=0,1,...,ν-1 is equivalent to the following incremental formulation: for almost every t ∈ [0, T ]

(P ν )            u ν (t) ∈ K(f ν (t)), F (u ν (t)), v - d dt ûν (t) + j(f ν (t), u ν (t), v) -j(f ν (t), u ν (t), d dt ûν (t)) ≥ b(f ν (t), u ν (t), v - d dt ûν (t)) ∀ v ∈ V, b(f ν (t), u ν (t), w -u ν (t)) ≥ 0 ∀ w ∈ K.
The existence of a solution of the problem P will be proved by using the following results and their proofs will be presented in a forthcoming paper.

Lemma 2.1. There exist a subsequence of (u ν , ûν ) ν , denoted by (u νp , ûνp ) p , and an element u ∈ W 1,2 (0, T ; V ) such that

u νp (t) u(t) in V ∀ t ∈ [0, T ], (18)
ûνp u in W 1,2 (0, T ; V ), (19)
d dt ûν u in L 2 (0, T ; V ). ( 20
)
Also, for all s ∈ [ 0, T ], we have u(s) ∈ K(f (s)) and

lim inf p→∞ s 0 j(f νp (t), u νp (t), d dt ûνp (t)) dt ≥ s 0 j(f (t), u(t), u(t)) dt. (21)
We have the following strong convergence and existence result.

Theorem 2.2. Under the assumptions ( 1)-( 3), ( 5)-( 16) every convergent subsequence of (u ν , ûν ) ν , still denoted by (u ν , ûν ) ν , and its limit u ∈ W 1,2 (0, T ; V ), given by lemma 2.1, satisfy the following properties:

u ν (t) → u(t) in V ∀ t ∈ [0, T ], ( 22
) ûν → u in L 2 (0, T ; V ), ( 23
)
and u is a solution of problem P .

Internal approximation and convergence results

In this section we shall study the approximation of problem P by using a convergence result for a method based on an internal approximation and a backward difference scheme. The full proofs will be presented in a forthcoming paper. First, we consider a semi-discrete approximation of (P ), which extends the classical internal approximations as presented in, e.g. [START_REF] Glowinski | Analyse numérique des inéquations variationnelles[END_REF], [START_REF] Glowinski | Numerical Methods for Nonlinear Variational Problems[END_REF]. For a positive parameter h converging to 0, let (V h ) h be an internal approximation of V, that is a family of finite-dimensional subspaces of V which satisfies:

there exist U ⊂ V such that U = V and ∀ v ∈ U, ∃ v h ∈ V h for each h, such that v h → v in V. (24)
Let (K h ) h be a family of convex cones with their vertices at 0 such that

K h ⊂ V h for all h and (K h ) h is an internal approximation of K, i.e. if v h ∈ K h for all h and v h v then v ∈ K, (25)
∀ v ∈ K, ∃ v h ∈ K h for each h, such that v h → v in V. (26)
Let (K h (g)) g∈V be a family of nonempty convex subsets of K h such that 0 ∈ K h (0) for all h, satisfying the following conditions:

if g n → g in V, v hn ∈ K h (g n ) and v hn → v h in V h then v h ∈ K h (g), ( 27
) if v h ∈ K h (g) for all h and v h v then v ∈ K(g) ∀ g ∈ V. ( 28
)
We assume that for all g ∈ V there exists an operator γ h (g,

•) : K h (g) → H such that γ h (0, 0) = 0 and for all g i ∈ V , v hi ∈ K h (g i ), i = 1, 2, γ h (g 1 , v h1 ) -γ h (g 2 , v h2 ) H ≤ k 1 ( g 1 -g 2 + v h1 -v h2 ). ( 29
) For all g ∈ V , let j h (g, •, •) : K h (g) × V h → R be a functional satisfying the following conditions for all g ∈ V : if v h ∈ K h (g) for all h, v h v in V and w h w in V then lim h→0 j h (g, v h , w h ) = j(g, v, w), (30)
for all h and

v h ∈ K h (g) j h (g, v h , •) is sub-additive, (31)
for all h and v h ∈ K h (g) j h (g, v h , •) is positively homogeneous, (32)

j h (0, 0, w h ) = 0 ∀ w h ∈ V h , (33)
and

if v h (t) ∈ K h (g(t)) for all h and t ∈ [0, T ], v h v in W 1,2 (0, T ; V ) then lim inf h→0 T 0 j h (g(t), v h (t), vh (t)) dt ≥ T 0 j(g(t), v(t), v(t)) dt (34) for all g ∈ C([0, T ]; V ), |j h (g 1 , v h1 , w h2 ) + j h (g 2 , v h2 , w h1 ) -j h (g 1 , v h1 , w h1 ) -j h (g 2 , v h2 , w h2 )| ≤ k 2 ( g 1 -g 2 + γ h (g 1 , v h1 ) -γ h (g 2 , v h2 ) H ) w h1 -w h2 (35) ∀ g i ∈ V, v hi ∈ K h (g i ), w hi ∈ V h , i = 1, 2.
From the properties of F, j h and K h , it follows that for all g ∈ V, d h ∈ K h , w h ∈ K h (g), the elliptic variational inequality:

u h ∈ K h F (u h ), v h -u h + j h (g, w h , v h -d h ) -j h (g, w h , u h -d h ) ≥ 0 ∀ v h ∈ K h
has a unique solution. Hence we can define a mapping S h g,d h :

K h (g) → K h by S h g,d h (w h ) = u h . We suppose that for all g ∈ V, d h ∈ K h S h g,d h (K h (g)) ⊂ K h (g). ( 36
)
For all g ∈ V, d h ∈ K h , we consider the following problems:

( Ph )

     u h ∈ K h (g), F (u h ), v h -u h + j h (g, u h , v h -d h ) -j h (g, u h , u h -d h ) ≥ b(g, u h , v h -u h ) ∀ v h ∈ V h , b(g, u h , z h -u h ) ≥ 0 ∀ z h ∈ K h ,

and

( Qh )

u h ∈ K h (g), F (u h ), v h -u h + j h (g, u h , v h -d h ) -j h (g, u h , u h -d h ) ≥ 0 ∀ v h ∈ K h .
We assume that if u h is a solution of ( Qh ), then u h is a solution of ( Ph ).

(37)

Remark 3. If u h satisfies ( Ph ), then u h obviously satisfies ( Qh ). Now, we consider the following equivalent semi-discrete problems.

Problem P h : Find u h ∈ W 1,2 (0, T ; V h ) such that (P h )          u h (0) = u 0h , u h (t) ∈ K h (f (t)) ∀ t ∈]0, T [, F (u h (t)), v h -uh (t) + j h (f (t), u h (t), v h ) -j h (f (t), u h (t), uh (t)) ≥ b(f (t), u h (t), v h -uh (t)) ∀ v h ∈ V h a.e. on ]0, T [, b(f (t), u h (t), z h -u h (t)) ≥ 0 ∀ z h ∈ K h , ∀ t ∈]0, T [, Problem Q h : Find u h ∈ W 1,2 (0, T ; V h ) such that (Q h )      u(0) = u 0h , u h (t) ∈ K h (f (t)) ∀ t ∈]0, T [, F (u h (t)), w h -u h (t) + j h (f (t), u h (t), w h -u h (t) + uh (t)) -j h (f (t), u h (t), uh (t)) ≥ 0 ∀ w h ∈ K h a.e. on ]0, T [, and
Problem Qh : Find u h ∈ W 1,2 (0, T ; V h ) such that ( Qh )            u(0) = u 0h , u h (t) ∈ K h (f (t)) ∀ t ∈]0, T [, F (w h ) -F (u h (t)) + j h (f (t), u h (t), w h -u h (t) + uh (t)) -j h (f (t), u h (t), uh (t)) ≥ α 2 w h -u h (t) 2 ∀ w h ∈ K h a.e. on ]0, T [.
where u 0h ∈ K h (f (0)) is the unique solution of the variational inequality

F (u 0h ), w h -u 0h + j h (f (0), u 0h , w h ) -j h (f (0), u 0h , u 0h ) ≥ 0 ∀ w h ∈ K h . ( 38
)
The full discretization of (P h ) is obtained by using an implicit scheme as in Section 2 for (P ). For u 0 h := u 0h and ι ∈ {0, 1, ..., ν -1}, we define u ι+1 h as the solution of the following problem.

Problem P ι hν : Find u ι+1 h ∈ K ι+1 h such that (P ι hν )      F (u ι+1 h ), v h -∂u ι h + j h (f ι+1 , u ι+1 h , v h ) -j h (f ι+1 , u ι+1 h , ∂u ι h ) ≥ b(f ι+1 , u ι+1 h , v h -∂u ι h ) ∀ v h ∈ V h , b(f ι+1 , u ι+1 h , z h -u ι+1 h ) ≥ 0 ∀ z h ∈ K h ,
where K ι+1 h := K h (f ι+1 ). As in Section 2, it follows that for all ι ∈ {0, 1, ..., ν -1} the problem (P ι hν ) is equivalent to each of the following variational inequalities:

find u ι+1 h ∈ K ι+1 h such that (Q ι hν ) F (u ι+1 h ), w h -u ι+1 h + j h (f ι+1 , u ι+1 h , w h -u ι h ) -j h (f ι+1 , u ι+1 h , u ι+1 h -u ι h ) ≥ 0 ∀ w h ∈ K h , find u ι+1 h ∈ K ι+1 h such that ( Qι hν )    F (w h ) -F (u ι+1 h ) + j h (f ι+1 , u ι+1 h , w h -u ι h ) -j h (f ι+1 , u ι+1 h , u ι+1 h -u ι h ) ≥ α 2 w h -u ι+1 h 2 ∀ w h ∈ K h .
From the relations (2), (29), (35), ( 12) and (36) it follows that the mapping

S h f ι+1 ,u ι h : K ι+1 h → K ι+1 h
is a contraction, so that (Q ι hν ) has a unique solution which is also the unique solution of (P ι hν ) and of ( Qι hν ), for all ι ∈ {0, 1, ..., ν -1}.

If we define the functions

         u hν (0) = ûhν (0) = u 0h and ∀ ι ∈ {0, 1, ..., ν -1}, ∀ t ∈]t ι , t ι+1 ], u hν (t) = u ι+1 h , ûhν (t) = u ι h + (t -t ι )∂u ι h ,
then for all ν ∈ N * the sequence of inequalities (P hι ν ) ι=0,1,...,ν-1 is equivalent to the following incremental formulation: for almost every t ∈ [0, T ]

(P hν )            u hν (t) ∈ K h (f ν (t)), F (u hν (t)), v h - d dt ûhν (t) + j h (f ν (t), u hν (t), v h ) -j h (f ν (t), u hν (t), d dt ûhν (t)) ≥ b(f ν (t), u hν (t), v h - d dt ûhν (t)) ∀ v h ∈ V h , b(f ν (t), u hν (t), w h -u hν (t)) ≥ 0 ∀ w h ∈ K h .
We have the analogue to theorem 2.2 in the finite-dimensional case.

Theorem 3.1. Assume that (2), ( 3), ( 13), ( 14), ( 27), ( 29), ( 31)-( 33), ( 35)-(37) hold. Then there exists a subsequence of (u hν , ûhν ) ν , still denoted by (u hν , ûhν ) ν , such that

u hν (t) → u h (t) in V ∀ t ∈ [0, T ], ( 39
) ûhν → u h in L 2 (0, T ; V ), ( 40
)
where u h is a solution of (P h ).

Using the previous theorem and a priori estimates for the solutions of (P h ) we have the following convergence and existence result. Theorem 3.2. Under the assumptions (2), (3), ( 13), ( 14), ( 24)-(37) there exists a subsequence of (u h ) h such that

u h (t) → u(t) in V ∀ t ∈ [ 0, T ], (41)
u h → u in L 2 (0, T ; V ), ( 42
) uh u in L 2 (0, T ; V ), ( 43
)
where u is a solution of (P ).

Using theorems 3.1 and 3.2, we obtain the following main approximation result.

Theorem 3.3. Under the assumptions of theorem 3.2, there exists a subsequence of (u hν ) hν such that

u hν (t) → u(t) in V ∀ t ∈ [0, T ], (44) uhν u in L 2 (0, T ; V ), (45)
where u ∈ W 1,2 (0, T ; V ) is a solution of (P ). Furthermore any cluster point of (u hν ) hν is a solution of (P ).

Subspace correction approximation

The aim of this section is to give, for a fixed time step ι, one-and two-level Schwarz methods for problem (Q ι ν ), prove their global convergence and estimate the convergence rate. First, we present three correction algorithms in an abstract Hilbert space, prove, under some assumptions, their global convergence and estimate the error. The domain decomposition methods are obtained from these general algorithms by associating to a domain decomposition some subspaces of a Sobolev space. In particular, the one-and two-level methods are obtained using the finite element spaces. In this case, we can estimate the convergence rate depending on the mesh and overlapping parameters.

Subspace correction algorithms

As in the previous section, we consider a Hilbert space V , and let V 1 , • • • , V m be some closed subspaces. We also consider a convex subset K ⊂ V , and assume that it satisfies the following assumption. Assumption 4.1. There exists a constant C 0 such that for any w, v ∈ K and w i ∈ V i with w

+ i j=1 w j ∈ K, i = 1, • • • , m, there exist v i ∈ V i , i = 1, • • • , m, satisfying w + i-1 j=1 w j + v i ∈ K for i = 1, • • • , m, (46)
v -w = m i=1 v i , (47)
and

m i=1 v i ≤ C 0 v -w + m i=1 w i . ( 48
)
This assumption looks complicated enough, but we will see in what follows that it is satisfied for many different convex sets in Sobolev spaces, in particular, for the convex sets K and K(g) in the previous section. In the proofs, v is the exact solution, w is the current iterate, and w i are the corrections on the subspaces V i , i = 1, . . . , m. In the case where the convex set is written as a sum of convex subsets or the problem has no constraint, (46) and (47) are always satisfied with w i = 0. In these cases, the above assumption is much more simple, and for this reason (48) usually is known without the extra terms given by w i . Now, we consider a functional ϕ : K × K → R and we assume that is convex and lower semicontinuous in the second variable, and

|ϕ(v 1 , w 2 ) + ϕ(v 2 , w 1 ) -ϕ(v 1 , w 1 ) -ϕ(v 2 , w 2 )| ≤k 1 k 2 ||v 1 -v 2 ||||w 1 -w 2 ||, ∀v 1 , v 2 , w 1 , w 2 ∈ K. ( 49
)
Also, we suppose that

Assumption 4.2. m i=1 [ϕ(u, w + i-1 j=1 w j + v i ) -ϕ(u, w + i-1 j=1 w j + w i )] ≤ϕ(u, v) -ϕ(u, w + m i=1 w i ) (50)
for any u ∈ K, and for v, w ∈ K and v i , w i ∈ V i , i = 1, . . . , m, in Assumption 4.1.

This assumption has been introduced for proof reasons. As in the previous section, we take a Gateaux differentiable functional F : V → R defined satisfying (2) and (3), we consider the problem of finding u ∈ K, the solution of the following quasi-variational inequality:

F (u), v -u + ϕ(u, v) -ϕ(u, u) ≥ 0 ∀ v ∈ K. ( 51
)
Since ϕ satisfies (49), with similar arguments as for problem (Q ι ν ), we can prove that problem (51) has a solution and it is unique.

Evidently, since ϕ is convex in the second variable, and F is differentiable and satisfies (4), problem (51) is equivalent with the minimization problem

u ∈ K : F (u) + ϕ(u, u) ≤ F (v) + ϕ(u, v) ∀ v ∈ K. ( 52
)
Also, in view of (4) we see that the solution u of (51) satisfies

α 2 v -u 2 ≤ F (v) -F (u) + ϕ(u, v) -ϕ(u, u) ∀ v ∈ K. (53)
A first algorithm corresponding to the subspaces V 1 , • • • , V m and the convex set K is written as Algorithm 4.1. We start the algorithm with an arbitrary u 0 ∈ K. At iteration n + 1, having u n ∈ K, n ≥ 0, we compute sequentially for i

= 1, • • • , m, the local corrections w n+1 i ∈ V i , u n+ i-1 m + w n+1 i ∈ K satisfying F (u n+ i-1 m + w n+1 i ), v i -w n+1 i + ϕ(u n+ i-1 m + w n+1 i , u n+ i-1 m + v i ) -ϕ(u n+ i-1 m + w n+1 i , u n+ i-1 m + w n+1 i ) ≥ 0, ∀ v i ∈ V i , u n+ i-1 m + v i ∈ K, (54)
and then we update

u n+ i m = u n+ i-1 m + w n+1 i .
A simplified variant of Algorithm 4.1 can be written as Algorithm 4.2. We start the algorithm with an arbitrary u 0 ∈ K. At iteration n + 1, having

u n ∈ K, n ≥ 0, we compute sequentially for i = 1, • • • , m, the local corrections w n+1 i ∈ V i , u n+ i-1 m + w n+1 i ∈ K satisfying F (u n+ i-1 m + w n+1 i ), v i -w n+1 i + ϕ(u n+ i-1 m , u n+ i-1 m + v i ) -ϕ(u n+ i-1 m , u n+ i-1 m + w n+1 i ) ≥ 0, ∀ v i ∈ V i , u n+ i-1 m + v i ∈ K, (55)
and then we update

u n+ i m = u n+ i-1 m + w n+1 i .
We can simplify Algorithm 4.1 even more as Algorithm 4.3. We start the algorithm with an arbitrary u 0 ∈ K. At iteration n + 1, having u n ∈ K, n ≥ 0, we compute sequentially for i

= 1, • • • , m, the local corrections w n+1 i ∈ V i , u n+ i-1 m + w n+1 i ∈ K satisfying F (u n+ i-1 m + w n+1 i ), v i -w n+1 i + ϕ(u n , u n+ i-1 m + v i ) -ϕ(u n , u n+ i-1 m + w n+1 i ) ≥ 0, ∀ v i ∈ V i , u n+ i-1 m + v i ∈ K, (56)
and then we update

u n+ i m = u n+ i-1 m + w n+1 i .
As for problem (51), we can prove that problems (54)-( 56) are equivalent with some minimization problems and they have unique solutions.

The following theorem proves that if k 1 k 2 is small enough in comparison with α and β, then Algorithms 4.1, 4.2 and 4.3 are convergent. Then, if u is the solution of problem (51), u n+ i m , n ≥ 0, i = 1, . . . , m, are its approximations obtained from one of Algorithms 4.1, 4.2 or 4.3, and

α 2 ≥ mk 1 k 2 + 2m(25C 0 + 8)βk 1 k 2 , (57)
then we have the following error estimations

F (u n ) + ϕ(u, u n ) -F (u) -ϕ(u, u) ≤ C 1 C 1 + 1 n F (u 0 ) + ϕ(u, u 0 ) -F (u) -ϕ(u, u) , (58)
u n -u 2 ≤ 2 α C 1 C 1 + 1 n F (u 0 ) + ϕ(u, u 0 ) -F (u) -ϕ(u, u) , (59)
where the constant C 1 > 0 depends on α, β, k 1 , k 2 , the number of subspaces m, and on the constant C 0 introduced in Assumption 4.1.

In a slightly more general form, this theorem has been proved in [START_REF] Badea | One-and two-level Schwarz methods for variational inequalities of the second kind and their application to frictional contact[END_REF] (see also [START_REF] Badea | One-and two-level multiplicative Schwarz methods for variational and quasi-variational inequalities of the second kind: Part I -general convergence results[END_REF]). In the case of Algorithm 4.1, the constant C 1 can be written as,

C 1 = C 2 /C 3 C 2 = βm(1 + 2C 0 + C 0 ε 1 ) + k 1 k 2 m(1 + 2C 0 + 1 + 3C 0 ε 2 ) C 3 = α 2 -k 1 k 2 (1 + ε 3 )m (60)
where

ε 1 = ε 2 = 2k 1 k 2 m α 2 -k 1 k 2 m , ε 3 = α 2 -k 1 k 2 m 2k 1 k 2 m ,
For Algorithms 4.2 and 4.3, the constant C 1 has similar expressions.

Convergence rates for the one-and two-level methods

Algorithms 4.1-4.3 can be viewed as multiplicative Schwarz methods in a subspace correction variant if the solution space is a Sobolev space and associate the subspaces to the subsets in a domain decomposition. Let us assume that the domain Ω is decomposed as

Ω = m i=1 Ω i (61)
If for instance V = H 1 0 (Ω) and we take V i = H 1 0 (Ω i ), i = 1, . . . , m, then the above algorithms are Schwarz methods.

The convergence rates given in Theorem 4.1 depend on the functionals F and ϕ, the number m of the subspaces and the constant C 0 introduced in Assumption 4.1. The number of subspaces can be associated with the number of colors needed to mark the subdomains such that the subdomains with the same color do not intersect with each other. Since this number of colors depends in general on the dimension of the Euclidean space where the domain lies, we can conclude that our convergence rate essentially depends on the constant C 0 .

We make now some remarks on the application of Algorithms 4.1-4.3 for the solution of problem (Q ι ν ). First, we see that we can take K ι+1 in the place of K in problem (Q ι ν ), i.e. this problem is equivalent with: find

u ι+1 ∈ K ι+1 such that (R ι ν ) F (u ι+1 ), w -u ι+1 + j(f ι+1 , u ι+1 , w -u ι ) -j(f ι+1 , u ι+1 , u ι+1 -u ι ) ≥ 0 ∀ w ∈ K ι+1 .
Indeed, since K ι+1 ⊂ K, it follows that the solution of (Q ι ν ) is also the solution of (R ι ν ). But, with similar arguments as for (Q ι ν ), we can prove that problem (R ι ν ) has a unique solution. It is proved in [START_REF] Badea | Convergence rate of a multiplicative Schwarz method for strongly nonlinear inequalities[END_REF] that if the convex set K has the property Property 4.1. If v, w ∈ K, and if θ ∈ C 0 (Ω), θ ∈ C 1 (Ω i ), i = 1, . . . , m, with 0 ≤ θ ≤ 1, then θv + (1θ)w ∈ K then Assumption 4.1 is satisfied with a C 0 depending on 1/δ, where δ is the overlapping parameter of the domain decomposition (61). In general, the convex sets K which are defined only by the values of the functions, but not by their derivatives, have the above property. The convex sets of one-obstacle type, like in contact problems (see [START_REF] Capatina | A class of implicit variational inequalities and applications to frictional contact[END_REF], for instance), or those of two-obstacle type have Property 4.1. Our problems (Q ι ν ), or (R ι ν ), are defined in an abstract Hilert space V , and we shall suppose in the following that their convex set K ι+1 has this property.

Since f ι+1 and u ι are fixed in problem (R ι ν ), writing

ψ(u, v) = j(f ι+1 , u, v -u ι ) ( 62
)
this functional has the properties asked of ϕ in problem (51), i.e. it is lower semicontinuous and convex in the second variable, and satisfies (49). But, such a ψ does not satisfy Assumption 4.2. In general, j is given by an integral, and we can overcome this difficulty by considering an approximation ϕ of ψ obtained by a numerical quadrature.

To conclude, writing u := u ι+1 and K := K ι+1 , problem (R ι ν ) (or (Q ι ν )) could be approximated by a problem (51) which satisfies the conditions of Theorem 4.1 if ψ in (62) could be approximated by a ϕ, which is lower semicontinuous and convex in the second variable, and satisfies (49) and Assumption 4.2. In this case, the Schwarz methods given by Algorithms 4.1-4.3 are globally convergent with a convergence rate depending on 1/δ.

We focus in the following on the application of Algorithms 4.1-4.3 for the solution of the finite element form of problem (Q ι ν ). Evidently, Theorem 4.1 can be applied for finite element spaces, too. More precisely, we study the convergence of the one-and two-level methods for this problem. As above, these methods are directly obtained from Algorithms 4.1-4.3. We shall see that Assumption 4.1 holds for closed convex sets K h satisfying a similar property with that given in Property 4.1, in particular for the discretized form of K ι+1 . Moreover, we are able to explicitly write the dependence of C 0 on the overlapping and mesh parameters. Also, we give some numerical approximations ϕ of the functional j for which Assumption 4.2 holds. Therefore, from Theorem 4.1, we can conclude that these methods globally converge for the discretized form of (Q ι ν ) if conditions ( 2) and (3) on F , and condition (11) on j hold. Moreover, from the dependence of C 0 on the mesh and domain decomposition parameters, we will conclude that the convergence rate is optimal, i.e. it is similar with that in the case of linear equations, for instance. The convergence rate of the two-level method depends very weakly on the mesh and domain decomposition parameters, and, for some particular choices, it is even independent of them.

One-level method

We consider a simplicial regular mesh partition T h (see [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF], p. 124, for instance) of mesh size h over the domain Ω ⊂ R d . The domain Ω is decomposed as in (61) with the overlapping parameter δ, and we assume that T h supplies a mesh partition for each subdomain Ω i , i = 1, . . . , m.

We associate to the decomposition (61), some functions

θ i ∈ C 0 ( Ω), θ i | τ ∈ P 1 (τ ) for any τ ∈ T h , i = 1, • • • , m, such that 0 ≤ θ i ≤ 1 on Ω, θ i = 0 on ∪ m j=i+1 Ω j \Ω i , θ i = 1 on Ω i \ ∪ m j=i+1 Ω j , (63)
and

|∂ x k θ i | ≤ C/δ, a.e. in Ω, ∀ k = 1, . . . , d (64)
Such functions θ i with the above properties exist (see [START_REF] Toselli | Domains Decomposition Methods -Algorithms and Theory[END_REF], p. 59, for instance). As in (64), we denote in the following by C a generic constant which does not depend on either the mesh or the domain decomposition. Now, writing the discretized variant of (R ι ν ) as a problem of the form (49), we can conclude from the following proposition that the error estimations in Theorem 4.1 hold for the one-level multiplicative Schwarz methods given by Algorithms 4.1-4.3 if F satisfies (2) and (3), and ϕ satisfies (49). In general, ϕ satisfies this condition if j satisfies [START_REF] Glowinski | Analyse numérique des inéquations variationnelles[END_REF]. Proposition 1. Assumption 4.1 holds for the piecewise linear finite element spaces, V = V h and V i = V i h , i = 1, . . . , m, and any convex set K h ⊂ V h having Property 4.2. Also, Assumption 4.2 on the functional ϕ of the form (67) is satisfied. The constant in (48) can be written as

C 0 = Cm 1 + m -1 δ , ( 69
)
where C is independent of the mesh and domain decomposition parameters.

The proof of this proposition can be found in [START_REF] Badea | One-and two-level Schwarz methods for variational inequalities of the second kind and their application to frictional contact[END_REF] or [START_REF] Badea | One-and two-level multiplicative Schwarz methods for variational and quasi-variational inequalities of the second kind: Part I -general convergence results[END_REF]. It is proved there that Assumptions 4.1 and 4.2 hold if we consider

v i = L h θ i (v -w - i-1 j=1 v j ) + (1 -θ i )w i , i = 1, . . . , m, (70)
where the functions θ i satisfy (63) and (64).

Two-level method

We consider two simplicial regular mesh partitions T h and T H on the domain Ω ⊂ R d of mesh sizes h and H, respectively. As in the previous section, we consider an overlapping decomposition (61), the mesh partition T h of Ω supplying a mesh partition for each Ω i , 1 ≤ i ≤ m. Also, we assume that the overlapping size is δ. In addition, we suppose that there exists a constant C, independent of both meshes, such that the diameter of the connected components of each Ω i are less than CH. The domain Ω may be different from

Ω 0 = τ ∈T H τ, (71)
but we assume that the mesh T h is a refinement of T H on Ω 0 , and if a node of T H lies on ∂Ω 0 then it also lies on ∂Ω. Finally, we suppose that there exists a constant C, independent of both meshes, such that dist(x, Ω 0 ) ≤ CH (72)

for any node x of T h . Now, besides the spaces V h and V i h , i = 1, . . . , m, defined in (65) and (66), we introduce the continuous, piecewise linear finite element space corresponding to the H-level,

V 0 H = v ∈ C 0 ( Ω0 ) : v| τ ∈ P 1 (τ ), τ ∈ T H , v = 0 on ∂Ω 0 , (73)
where the functions v are extended with zero in Ω\Ω 0 . The convex set K h ⊂ V h is defined as a subset of V h having Property 4.2.

The two-level Schwarz methods are also obtained from Algorithms 4.1-4.3 in which we take V = V h , K = K h , and the subspaces

V 0 = V 0 H , V 1 = V 1 h , V 2 = V 2 h , . . . , V m = V m h . As in the previous section, the spaces V h , V 0 H , V 1 h , V 2
h , . . . , V m h , are considered as subspaces of H 1 . We note that, this time, the decomposition of the domain Ω contains m overlapping subdomains, but we use m + 1 subspaces of V , V 0 , V 1 , . . . , V m , in Algorithms 4.1-4.3. Assuming that F satisfies (2) and ( 3) and ϕ satisfies (49), if we prove that Assumption 4.1, written for m + 1 subspaces, is satisfied for the previous choice of the convex set K and the subspaces V 0 , V 1 , . . . , V m of V , and the functional ϕ of the form (67) satisfies Assumption 4.2, we can conclude that these Algorithms 4.1-4.3 converge. To this end, we consider the operator I H : V h → V 0 H , which has been introduced in [START_REF] Badea | Convergence rate of a Schwarz multilevel method for the constrained minimization of non-quadratic functionals[END_REF] and has the following properties (see Lemma 4.2 in [START_REF] Badea | Convergence rate of a Schwarz multilevel method for the constrained minimization of non-quadratic functionals[END_REF]) for any v ∈ V h :

I H v -v 0 ≤ CHC d (H, h)|v| 1 (74)
and

I H v 0 ≤ C v 0 and |I H v| 1 ≤ CC d (H, h)|v| 1 , (75)
where

C d (H, h) =      1 if d = 1 ln H h + 1 1 2 if d = 2 H h 1 2 if d = 3. (76)
Moreover, for any x ∈ Ω, we have

0 ≤ I H v(x) ≤ v(x) if v(x) > 0, 0 ≥ I H v(x) ≥ v(x) if v(x) < 0,
and I H v = 0 on τ ∈ T H if there exists a x ∈ τ such that v(x) = 0 for any v ∈ V h . Consequently, writing

θ v (x) = I H v(x) v(x) if v(x) = 0 0 if v(x) = 0, (77)
then θ v ∈ C 0 ( Ω), θ v | τ ∈ C 1 (τ ) for any τ ∈ T h , 0 ≤ θ v ≤ 1, and

I H v = θ v v
for any v ∈ V h . Now, we can prove the following proposition which shows that the constant C 0 in Assumption 4.1 is independent of the mesh and domain decomposition parameters if H/δ and H/h are constant. Proposition 2. Assumption 4.1 is satisfied for the piecewise linear finite element spaces V = V h and V 0 = V 0 H , V i = V i h , i = 1, . . . , m, defined in (65), (66), and (73), respectively, and any convex set K = K h with Property 4.2. Also, Assumption 4.2 for the functional ϕ of the form (67) is satisfied. The constant in (48) of Assumption 4.1 can be taken of the form

C 0 = Cm 1 + (m -1) H δ C d (H, h), ( 78
)
where C is independent of the mesh and domain decomposition parameters, and C d (H, h) is given in (76).

A detailed proof of this proposition can be found in [START_REF] Badea | One-and two-level Schwarz methods for variational inequalities of the second kind and their application to frictional contact[END_REF] or [START_REF] Badea | One-and two-level multiplicative Schwarz methods for variational and quasi-variational inequalities of the second kind: Part I -general convergence results[END_REF]. By means of I H and the functions θ i , i = 1, . . . , m, with properties (63) and (64), we define v 0 = w 0 + I H (vww 0 ).

and

v i = L h (θ i (v -w -v 0 - i-1 j=1 v j ) + (1 -θ i )w i ), (80)
for i = 1, . . . , m. Using properties (74) and (75) of the operator I H , it is proved that v 0 , v 1 , . . . , v m , defined in (79) and (80), satisfy Assumption 4.1 with the constant C 0 given in (78). Also, using (68) in which θ is replaced by θ i or by functions of the type θ v , defined in (77), it is proved that Assumption 4.2 hold, too.

Remark 4. In this section, we have assumed that the functional ϕ is of the form (67). With similar proofs, we can show that Propositions 1 and 2 also hold if we replace the functional ϕ(u, v) of the form (67) with

ϕ(u, v) = κ∈N h I κ (φ(u(x κ ), v(x κ ))) = k∈N h φ κ (u, v), (81)
where φ : R×R → R, I κ : R → R, κ ∈ N h , are continuous and I κ (φ(r, •)) : R → R, κ ∈ N h , are convex functions for any r ∈ R. For the writing simplicity, we have written φ κ (u, v) = I(φ(u(x κ ), v(x κ ))), κ ∈ N h . In general, (67) or (81) represent numerical approximations of some integrals. Concerning to condition (49) imposed on ϕ of the form (67) or (81), we have to check it for each particular problem we solve. The results of this subsection have referred to problems in H 1 with Dirichlet boundary conditions. We point out that similar results can be obtained for problems in (H 1 ) d or problems with mixed boundary conditions.

Theorem 4 . 1 .

 41 Let us assume that Assumptions 4.1 and 4.2 are satisfied.

Acknowledgments

The support of this work by "Laboratoire Européen Associé CNRS Franco-Roumain de Mathématiques et Modélisation" LEA Math-Mode is gratefully acknowledged.

We consider the piecewise linear finite element space

and also, for i = 1, . . . , m, we take

as some subspaces of V h corresponding to the domain decomposition Ω 1 , . . . , Ω m . The spaces V h and V i h , i = 1, . . . , m, are considered as subspaces of H 1 . We denote by • 0 the norm in L 2 , and by • 1 and | • | 1 the norm and seminorm in H 1 , respectively.

We have assumed that the convex set K ι+1 of problem (R ι ν ) is a particular case of a convex set K with Property 4.1. We consider that the discretized form of K = K ι+1 is defined as a subset K h ⊂ V h which satisfies a similar property,

Above, we have denoted by L h the P 1 -Lagrangian interpolation operator which uses the function values at the nodes of the mesh T h . This property is satisfied, in general by the convex sets which are defined by the values of the functions at the mesh nodes. Like Property 4.1, Property 4.2 is satisfied for convex sets of the one-or two-obstacle type.

Let us consider some continuous functionals I κ : L 2 (Ω) → R, and assume that the functional ϕ, the finite element approximation of ψ in (62), is of the form ϕ(u, v) =

where N h is the set of nodes of the mesh partition T h , φ :

is continuous, and we assume that for any u ∈ K h

are convex functions. For the ease of notation, for all κ ∈ N h we have written φ κ (u, v) = I κ (φ(u, v(x κ ))). We see that φ κ can be viewed as functionals

for any u, v, w ∈ K h , and any function θ : Ω → R with the properties θ ∈ C 0 ( Ω), θ| τ ∈ C 1 (τ ) for any τ ∈ T h , and 0 ≤ θ ≤ 1.