Luis A Salomón

Jean-Claude Fort
email: fort43@gmail.com

Li-Vang Lozada

Average Competitive Learning Vector Quantization

Keywords: Optimal quantization, Stochastic algorithms, Asymptotic convergence, Kushner-Clark method. 2000 MSC: 94A29, 62L20, 60F15, 90C51, 60F05, 65D30

We propose a new algorithm for vector quantization:Average Competitive Learning Vector Quantization(ACLVQ). It is a rather simple modification of the classical Competitive Learning Vector Quantization(CLVQ). This new formulation gives us similar results for the quantization error to those obtained by the CLVQ and reduce considerably the computation time to achieve the optimal quantizer. We establish the convergence of the method via the Kushner-Clark approach, and compare the two algorithms via the central limit Theorem. A simulation study is carried out showing the good performance of our proposal.

Introduction

The quantization problem comes from the theory of signal processing in electrical engineering in the late 1940's. More precisely the quantization of probability distributions is related to the best approximation of a ddimensional probability distribution P by a discrete probability distribution with finite support. The mathematical aspects of quantization are treated extensively in ? ]. This theory has been applied in several areas such as cluster analysis, pattern recognition, finance and numerical probability. For a detailed exposition on the applications of vector quantization we refer the reader to ? ], ? ], ? ], ? ] and references therein.

Searching an optimal quantizer is one of the main interest on this framework. Only for a few cases the close form is known. It is posible to compute the optimal quantizer at level n in one dimension as the solution of a stationarity equation either by a zero search method (Newton-Raphson gradient descent) or a fixed point procedure (like the Lloyd I procedure). In higher dimensions, deterministic gradient descent methods become intractable and one uses Lloyd I procedure and/or stochastic procedures to compute optimal quantizers (see for instance ? ]).

The basic stochastic approximation algorithms was introduced in the 1950s by the works of Robbins and Monro. The stochastic gradient method (i.e. a method of minima searching using the gradient) is based on the integral representation of the criterion to be optimized and can be seen as a particular case of the Robbins-Monro algorithm (see ? ? ] for instance).

The Competitive Learning Vector Quantization (CLVQ) is an on-line algorithm that can be seen as a particular case of the "stochastic gradient method with decreasing step". A good reference here is ? ]. We also refer the reader to ? ] where the authors give a detailed exposition of the application of the CLVQ. In ? ], ? ], ? ] and references therein we can found several variations around the CLVQ method.

The classical formulation of the CLVQ method has two phases: competitive and learning. The first one is the most time consuming, because it uses at each step a closest neighbor search to find the nearest point to the random vector generated at each iteration of the method. The second phase only actualizes the value of the quantizer. The performance of a vector quantization system depends on the composition of the codebook, the criteria to find an optimal quantizer and the calculations performed in order to find the best solution.

We propose a close algorithm to the CLVQ: Average Competitive Learning Vector Quantization. The advantage of using our proposal lies in the fact that we reduce the computation time for the obtention of optimal or stationary quantizers. We emphasize that our proposal attains similar results regarding the quantization error to those obtained by the CLVQ. Our method introduce a slight Lloyd modification in competitive phase to find the "winning index" which allows to reduce considerably the time and amount of calculations.

The rest of article is organized as follows. Next section describes the CLVQ and ACLVQ algorithms. There we recall briefly the approach of Kushner-Clark in the framework of stochastic algorithms. Section 3 is devoted to the proof of the convergence of ACLVQ algorithm. Section 4 shows the performance of our proposal by numerical experiments. Finally in Section 5 the conclusions are presented.

2. The Average Competitive Learning Vector Quantization.

Quadratic vector quantization.

Let us denote by P the target probability distribution and O = {x ∈ (R d ) n , ∀i = j, x i = x j }. We assume that P is diffuse, with a moment of order 2 + δ, δ > 0. This assumption ensures that P(x ∈ O) = 1. Throughout this work we will assume that our algorithms lives in O.

Let X be a square integrable R d -valued random vector with distribution P defined on a probability space (Ω, F, P). For n ∈ N, let Q n be the set of all Borel measurable maps f :

R d → α = {x 1 , • • • , x n } ⊂ R d with |f (R d )| ≤ n (where |A| denote the cardinal of the set A). The elements of Q n are called n-quantizers.
Among all Borel functions f ∈ Q n the one with the best properties is that obtained by the rule of closest neighbour:

f (X) = n i=1 x i 1 C i (X) (X), where {C i (X)} 1≤i≤n is a Voronoï partition for R d satisfying C i (X) ⊂ {x ∈ R d : ||x -x i || = min 1≤j≤n ||x -x j ||}.
Every n-quantizer induce an error.

Definition 2.1. The n-th quadratic quantization error for X is defined by

e n (X) = inf (E min a∈α ||X -a|| 2 ) 1/2 : α ⊂ R d , |α| ≤ n .
We said that α is the optimal n-quantizer for X if

e n (X) = E min 1≤j≤n ||X -x j || 2 1/2 .
Finding an optimal quantizer could be a hard task for some distributions.

Definition 2.2. Stationary quantizers. A n-quantizer X α is called stationary if any of the two following conditions are satisfied:

1. If X α is a nearest neighbour projection that satisfies

X α = E(X| X α ).
(1)

2. If X α is a critical point of the distortion, i.e.:

∇e n ( X α ) = 0.

Obviously an optimal n-quantizer is always n-stationary.

The classical Competitive Learning Vector Quantization

In applications, one generally prefers recursive algorithms, owing to their relative computational simplicity. Each successive estimate is obtained as a simple function of the last estimate and the current observation. Competitive Learning Vector Quantization is a classical example of these recursive methods. The general formulation of the CLVQ algorithm can be described briefly as follows.

CLVQ scheme.

1. Generate the initial n-quantizers x 1 0 , • • • , x n 0 . 2. Generate (independently) a random vector ξ with distribution P. 3. Find the winning index (competitive phase)

i 0 (k + 1) = arg min i x i k -ξ 2 .
4. Update the quantizers (learning phase)

   x i k+1 = x i k -γ k+1 (x i k -ξ), i = i 0 (k + 1) x i k+1 = x i k , i = i 0 (k + 1)
.

5. Repeat steps 2-4 until to satisfy some convergence criteria.

As we mentioned before this method is a "stochastic gradient method with decreasing step". Thus the previous scheme could be written under that theory.

Definition 2.3. Stochastic gradient method. Let g be a twice differentiable function from E to R such that dg has an integral representation on E with respect to P, i.e.:

dg(x) = R d dG(x, ξ)P(dξ), with dG : E × R d → R, dG(x, •) ∈ L 1 (P) for every x ∈ E.
The stochastic gradient method is defined by a triplet of sequences ((X k ) k≥0 , (ξ k ) k≥1 , (γ k ) k≥1 ) with values respectively in E, R d and R + satisfying for every k ≥ 1 that

X k+1 = X k -γ k+1 dG(X k , ξ k+1 ),
where (ξ k ) k≥1 are i.i.d random variables with distribution P and (γ k ) k≥1 is a decreasing positive deterministic sequence going to zero such that γ k = +∞.

The theorem which ensures the convergence of such algorithms reads as follows.

Theorem 2.1. Convergence a.s. Let g : E → R + be a continuously differentiable function whose differential dg admits an integral representation on E with respect to P

dg(x) = R d dG(x, ξ)P(dξ).
Assume that dg and dG satisfy

lim |x|→+∞ g(x) = +∞ and dg is Lipschitz continuous, dg(x) = O(g(x)), when |x| → +∞. Let ((X k ) n≥0 , (ξ k ) k≥1 , (γ k ) k≥1
) be a stochastic gradient method with a positive gain parameter sequence satisfying

k≥1 γ k = +∞, k≥1 γ 2 k < +∞.
Then g(X k ) a.s. converges to some nonnegative random variable g ∞ and X k a.s. converges toward some random connected component

χ * of {dg = 0} ∩ {g = g ∞ }.
In particular, if {dg = 0} is a finite set, with probability 1 it exists an x * ∈ {dg = 0} such that:

X k → x * when k → +∞.
In general we can not ensure that the algorithm converges to a global minimum. In higher dimension, uniqueness of stationary quantizers clearly often fails (see ? ]), so that the CLVQ reaches, a priori, only locally optimal quantizers. Even if using a procedure of Lloyd the numerical precision can be improved, there is no reason to suppose that the obtained quantizers are the optimal ones.

One crucial fact for this kind of algorithms is the choice of the gain sequence. Usually some heuristic methods are made attending to the distribution P to quantize (see ? ]). For instance, a gain sequence (γ k ) k≥1 usually used is

γ k = γ 0 a a + (γ 0 bk) c ,
for some constants γ 0 , a, b and c. That is the gain sequence used in the numerical experiments presented in Section 4.

When we know the a.s. convergence of the algorithm toward a minimum, we may look for the speed of convergence. It is given by the following central limit theorem (see e.g [? ]):

Theorem 2.2. Central Limit Theorem The assumptions are those of the theorem(2.1), plus:

• g is C 2 and S = d 2 g(x ) is definite positive, with smallest eigenvalue λ 1 .

• Denoting ε k+1 = dG(X k , ξ k+1 ) -dg(X k ), we assume that it exists a limiting variance for ε k :

lim k→∞ E(ε k ε t k |F k ) = V,
where V is definite positive. Here F k is the σ-algebra of the past, generated by (X 0 , (ξ l ) 1≤l≤k ).

• for some δ > 0, sup k ε k 2+δ < ∞.
Then on the event

X k converges toward x * , √ γ k (X k -x * )
weakly converges toward a N (0, Σ), a normal distribution with mean 0 and variance Σ. The matrix Σ is the solution of the following equation:

if γ k ∼ α k+β with 2α > 1 λ 1 , (S - 1 2α I)Σ + Σ(S - 1 2α I) t = V, (2)
if γ k = o( 1 k ) SΣ + ΣS t = V.
(3)

The Average Competitive Learning Vector Quantization

The CLVQ only changes one quantizer at each iteration, and thus needs a large number of simulations to achieve good results. Our proposal is a modification of the competitive phase of the CLVQ method. This new procedure introduce a sort of small Lloyd in this phase. The ACLVQ generates a set of N random vectors ξ instead of one as the CLVQ method do. Our proposal requires less time to achieve quantization errors of the same order despite it use the same number of simulations needed by the CLVQ. The underlying idea could be described as follows ACLVQ scheme.

1. Generate the initial n-quantizers x 1 0 , • • • , x n 0 . 2. Assign x i = x i 0 for all i. 3. Generate a set of i.i.d. random vectors ξ = (ξ 1 , • • • , ξ N ) (
with the same distribution P). 4. Calculate the distance matrix between the n-quantizer and the random vectors (competitive phase). 5. Identify the "Voronoï region" for each x i in the step k using ξ, i.e.:

I i k = {j : ξ j -x i k < ξ j -x m k , ∀ m = i}, i = 1, • • • , n, j = 1, • • • , N.
6. Calculate the mass centroid of the "Voronoï region" associated to each

I i k ξ i = i∈I i k ξ i N i k 1 I i k =∅ , with N i k = |I i k |,
where |A| denotes the cardinal of the set A. 7. Update the quantizers (learning phase)

   x i k+1 = x i k -γ k+1 (x k i -ξ i ), I i k = ∅ x i k+1 = x i k , I i k = ∅
where γ k > 0 is the gain sequence as usual. 8. Repeat steps 2-7 until reach some convergence criteria.

It is worth pointing out that this new method is not a stochastic gradient method with decreasing step as the CLVQ, even if it is very close to. The convergence theorem stated for this kind of algorithms can not be used directly. For that reason we need to write our proposal in the Kushner-Clark settings to obtain the asymptotic behavior of the ACLVQ solution.

The Kushner-Clark theorem has been used to describe the behavior of stochastic algorithms. It is also called method of ordinary differential equation and was introduced in ? ]. The most important results on a.s. convergence of stochastic algorithms is due to the mentioned theorem. The main idea comes from the link between the behavior of the solution of the ODE ( ẋ = -h(x)) and the convergence of each sample path of the stochastic algorithm.

Let us consider the general recursive algorithm:

   X 0 ∈ R d , X k+1 = X k -γ k+1 f (X k , ξ k+1 ), (4)
where k ∈ N, (ξ k ) k≥1 , (X k ) k≥0 are two sequences on R d , γ k is the gain of the method (i.e. γ k 0) and H some general function, defined as f (X k , ξ k+1 ) = H(X k , ξ k+1 ) -η k+1 where η k is a noisy term such as η k → 0 when k → +∞.

Let us define the mean function h as h(X k ) = E(H(X k , ξ k+1 )). In the Kushner-Clark approach the convergence of (X k ) k≥1 is "conditional a.s." and related to the average ODE, which describes the mean behavior of the stochastic algorithm. For instance the limiting points x * of the algorithm in (4) lies in {h = 0}.

Let us define X = (X i ) 1≤i≤n , where X i ∈ R d for all i = 1, • • • , n and γ k is the gain sequence as before. The general model under the Kushner-Clark settings can be written as:

X i k+1 = X i k -γ k+1 H i (X k , ω k+1 ), (5)
where

H i (X k , ω k+1 ) is a function from R d × R N d to R, (ω k ) k≥1 , is a sequence of random vectors in R N d : ω k = (ω 1 k , • • • , ω N k ), with ω j k i.i.d
∼ P for all k and j = 1, • • • , N , N ∈ N. Therefore ω k ∼ P N = N j=1 P and obviously X k is of the same kind as X.

The theorem which assures the convergence a.s of the ACLVQ to an optimal or stationary quantizer reads as follows

Theorem 2.3. ACLVQ convergence. If the sequence X k = (X 1 k , • • • , X n k
) defined by the stochastic algorithm given in equation ( 5) lives in a compact subset of O, then it converges P-a.s to a limiting point

x = (x ,1 , • • • , x ,n ) (a stationary quantizer) of the function h = (h 1 , • • • , h n ), where h i (X) = E [H i (X k , ω k+1 )|X k = X].
The proof of previous theorem is left to the next section.

ACLVQ convergence

We use the Kushner-Clark theorem to prove the asymptotic convergence of the solution of (5) toward one "optimal" or stationary quantizer of the distribution.

Preliminaries

Let us rewrite the general equation ( 4)

   X 0 ∈ R d , X k+1 = X k -γ k+1 H(X k , ξ k+1 ) + γ k+1 η k+1 . Defining ∆M k+1 = h(X k ) -H(X k , ξ k+1 ) the previous formulation reads as follows    X 0 ∈ R d , X k+1 = X k -γ k+1 h(X k ) + γ k+1 (∆M k+1 + η k+1 ). (6)
We state the following assumption 4) is bounded and the assumption R holds, then

(R) ≡        (i.) +∞ s=1 γ s ∆M s < +∞ (ii.) lim s→∞ η s = 0 . Theorem 3.1. Kushner-Clark. If x * is an equilibrium point of the ODE ẋ = -h(x), h a continuous function, the sequence (X k ) k≥1 in (
X k → x * , as k → +∞, on the event A x * K = {X k ∈ K ⊂ G x * infinitely often}, where K is a compact set.
The Kushner-Clark theorem needs the boundedness of the solutions (X k ) k≥1 and the continuity of h. In some cases the choice of h leads to the boundedness of the solution. The continuity requirement for h could be relaxed using Lyapunov functional (see for instance ? ]).

Proof of the Theorem 2.3

Let us begin with the following lemma. We will denote p i (x) = P(C i (x) and

q i (x) = (1 -(1 -P(C i (x))) N ). Lemma 3.1. The function h = (h 1 , • • • , h n ) (h(x) is defined for x = (x 1 , • • • , x n ), x i ∈ R d for all i = 1, • • • , n) is P-a.
s. continuous and its zeros are exactly stationary quantizers for the distribution P, since h writes:

h i (x) = q i (x) p i (x) dg i (x). ( 7
)
here dg i is the derivatives of g with respect to x i .

Proof. Firstly we prove the second statement of the lemma. We have that

H i (X k , ω k+1 ) =      X i k - N j=1 1 C i (X k ) (ω j k+1 )ω j k+1 N i k+1 if N i k+1 > 0 X i k if N i k+1 = 0,
where

N i k+1 = N j=1 1 C i (X k ) (ω j k+1 ), N i k+1
is the number of times that a random vector lies in the Voronoï region of X i k . Furthermore

P(N i k+1 = 0) = P(ω j k+1 / ∈ C i (X k )), ∀j = 1, 2, • • • , N = (1 -P(C i (X k ))) N , thus P(N i k+1 > 0) = 1 -(1 -P(C i (X k ))) N .
We use ω instead ω k in order to simplify the notations.

h i (X k ) = X i k P(N i k+1 > 0) - {N i k+1 >0} 1 N i k+1 N j=1 1 C i (X k ) (ω j k+1 )ω j k+1 P N (dω).
Denoting

Ω r = {ω 1 k+1 ∈ C i (X k ), • • • , ω r k+1 ∈ C i (X k ), ω r+1 k+1 / ∈ C i (X k ) • • • , ω N k+1 / ∈ C i (X k )},
and knowing that ω j k i.i.d ∼ P it follows

h i (X k ) = X i k P(N i k+1 > 0) - N r=1 1 r N r r j=1 Ω r 1 C i (X k ) (ω j k+1 )ω j k+1 P N (dω) = X i k - E1 C i (X k ) (ω 1 )ω 1 P(C i (X k )) P(N i k+1 > 0). Set U i (X k ) = E1 C i (X k ) (ω)ω P(C i (X k )) = E [X|X ∈ C i (X k )] , ∀i = 1, 2, • • • , n,
where ω ∼ P. We have

h i (X k ) = E [H i (X, ω k+1 )|X = X k ] = P(N i k+1 > 0)(X i k -U i (X k )).
This gives formula ( 7), and it is quite obvious that h i (X k ) = 0 if and only if Proof of Theorem 2.3. By Lemma 3.1 the proof is a direct consequence of the Kushner-Clark theorem.

X i k = U i (X k ). Which means that the zeros of h = (h 1 , • • • , h n ) are

About the asymptotics of the CLVQ and ACLVQ

The Central Limit Theorem (theorem 2.2) straightforwardly applied to CLV Q. For ACLV Q we need again a theorem a bit more general since it is only close to a stochastic gradient algorithm. Nevertheless, the assumptions and result are identical to that of theorem 2.2, replacing dg by h and d 2 g(x * ) by the derivatives of h (see [? ]). The Hessian d 2 g has been computed by Pollard and Fort-Pagès (see ([? ]),([? ]).

By noticing that dg(x * ) = 0 we easily obtain dh(x * ) :

dh(x * ) = diag( q i (x * ) p i (x * ) )d 2 g(x * ).
Now in order to compare the asymptotics of the two algorithms, we need a rescaling. Indeeed the ACLVQ use N times more samples at each iteration. So we will compare CLVQ at step k with ACLVQ at step [k/N ], where [u] is the entire part of u ∈ R. Of course the gains γ are to be rescaled too: if the gain is γ k at iteration k for CLVQ, it will be γ kN at iteration k for ACLVQ. Denoting for convenience the ACLVQ X N we thus obtain the following asymptotics.

For CLVQ :

√ γ k (X k -x * ) weakly converges toward a N (0, Σ). For ACLVQ γ [ k N ]N (X N [ k N ] -x * ) weakly converges toward a N (0, Σ N ). As γ [ k N ]N ∼ γ k
, it remains to compare the two covariance matrices Σ and Σ N . Let us denote V the asymptotic variance matrix of the noise for CLVQ (notice that V is a diagonal matrix), and N i , 1 ≤ i ≤ N a multinomial random variable with parameter (p i (x * ) 1≤i≤N ). As x * is a stationary quantizer, one can verify that the asymptotic variance of the noise for ACLVQ is:

V N = diag(E 1 N i > 0 p i (x * )N i )V.
The equations ( 2) and (3) remain valid, but dh(x * ) is no more symetric. Then, for instance in the case of equation ( 3), the matrix Σ N satisfies:

diag( q i (x * ) p i (x * ) )d 2 g(x * )Σ N +Σ N [diag( q i (x * ) p i (x * ) )d 2 g(x * )] t = diag( 1 p i (x * ) E 1 N i > 0 N i )V. Now we notice that E 1 N i >0 N i ≤ P(N i > 0) = q i (x * ), and setting D = diag( q i (x * ) p i (x * )
) we obtain:

Dd 2 g(x * )Σ N + Σ N d 2 g(x * ) t D ≤ DV,
as semi definite matrices. Introducing D 1 2 , using equation (3) for Σ we deduce that the matrix:

M = D 1 2 (d 2 g(x * )Σ N -d 2 g(x * )Σ)D -1 2 + (D 1 2 d 2 g(x * )Σ N -d 2 g(x * )ΣD -1 2 ) t ,
is semidefinite negative. From this we conclude that d 2 g(x * )Σ N -d 2 g(x * )Σ has negative eigenvalues, which yields Σ N ≤ Σ as semi definite matrices.

We now conclude that ACLVQ with the same number of samples is asymptotically more accurate than CLVQ.

This result as to be softened by the fact that when N is large, the asymptotic of ACLVQ is reached much longer after CLVQ. So choosing N too large would cause ACLVQ to fail with respect to CLVQ.

Numerical results

In this section, numerical comparisons between the classical CLVQ and our proposal are carried out in order to investigate the accuracy and computation time of each method. Our main interest is to illustrate the advantage of the ACLVQ regarding the computation time.

We carry out numerical experiments with distributions in dimension one and two. In dimension one it is not practical using the ACLVQ because determinist algorithms like the Newton one perform very well (see ? ]). However, these runs let us to evaluate the performance of the ACLVQ and the classical CLVQ for the standard Gaussian distribution. In dimension 2 we make the major portion of the numerical simulations. In this setting we consider 4 different distributions. Only for the uniform distribution in the unit cube for any dimension and particular sizes of the quantizer it is known the exact value of the distortion error. Therefore we consider these distributions in both dimensions concerned. More precisely, for n = k 2 , in ? ] is proved that

e 2 n (U ([0, 1] d )) = n -2/d d 12 .
In Table 1 we show the choice of the parameters of the probability distributions concerned in the study in dimension two.

Distribution Parameters

N (0, I 2 ) N (0, Σ 2 ) Σ 2 = (σ ij ) 1≤i,j≤2 : σ 11 = 7.1, σ 22 = 1.7 and σ ij = 0 if i = j U ([0, 1] 2 ) Exp(Ψ 2 ) Ψ 2 = (µ ij ) 1≤i,j≤2 : µ ii = 1 and µ ij = 0 if i = j
Table 1: Probability distributions in dimension 2 used in the study.

In the competitive phase we use five different values of N > 1, the number of random vector generated in each iteration of the method. The value N = 1 corresponds to the CLVQ. For each, distribution, n and N we run 100 repetitions.

To make a fair comparison between both methods we assure generate similar number of random vectors for each case. More precisely, we compare the runs of the CLVQ having total number of iteration K with those of the ACLVQ of N simulations of the random vector by generation and total number of generations K/N . Hence, both methods use the same "total information": K copies of the random vector. We set the same initial vector x 0 for each n and all N . For the gain sequence

γ k = γ 0 a a + (γ 0 N bk) c
we set a = 1, b = 0.1 and γ 0 = 0.2. In Table 2 we present the parameter values chosen in the study. All computations were made with MATLAB 7.11.0.584 R2010b in a computer with the following characteristics: CPU Pentium(R) Dual Core T4500 2.30 GHz and 3Gb RAM.

Quantization error

First of all, we point out that the computation of the quantization error in all the cases was made by Monte Carlo simulation using a dataset of 10 6 random vector copies.

In ? ] the authors reported quantization errors for the n-quantizer of the standard Gaussian distribution for n = 2, 3, • • • , 8. Our algorithm reaches results quite similar to those. The simulation study shows that our method exhibits a small standard deviation (of order of 10 -4 ) over the 100 simulation for each distribution and each value of N . Similar behavior was observed for the one dimensional uniforme distribution on [0, 1]. These results support our assumption that the ACLVQ has a good accuracy.

In Figure 1 and Figure 2 we display the boxplots of the quantization errors after 100 repetitions of the experiment for each of values of N chosen and n = 10 and n = 50 respectively. In Table 3 we shows the rest of the results obtained. As can be observed the results obtained by the CLVQ (N = 1) and by the ACLVQ are close. The standard deviation is small. Moreover, the quantizers obtained for each case had the same spatial configuration (graphs not shown).

Computation time

Now, we turn our attention to the computation time. The major advantage of our proposal is the reduction of the computation time to obtain optimal quantizers. The results regarding the computation time can be seen in Table 4. As we can see the ACLVQ reduces considerably the computation time when N takes large values. For instance, the ACLVQ for n = 10 and N = 200 achieves the optimal quantizer 30 times faster than the CLVQ.

One could think that increasing N is all what one need to do in order to improve the performance of algorithm. However, in experiments we do not show here it was observed a decline of the accuracy of the algorithm. This is mainly due to the fact that the ACLVQ can be trapped more easily in local minima. Similar behavior is observed in the Lloyd's method I when the initial set is far from the global minimum [? ]. Actually, the choice of N is a tradeoff between accuracy and time consuming.

n = 10 n = 50 N D 1 D 2 D 3 D 4 N D 1 D 2 D 3 D 4 1 66,
The gain in time is due to two factors: the first one and the most important is related with the number of computations that take place on the learning phase: the ACLVQ is able to update more quantizers at each iteration than the CLVQ and, at the same time, in the overall process the ACLVQ update less times. The number of iterations needed by the ACLVQ is smaller. This fact relies on the Lloyd modification of the ACLVQ. The computation time needed for this modification is negligible for the whole process. It is important to that the number of random vector simulations generation for both algorithms is similar. The second aspect concerns the procedure used for the computation of the distance matrix in the competitive phase (see Appendix A) which exploits very efficiently the concept of vectorizing algorithms of MATLAB.

Some considerations in conclusion

The Average Competitive Learning Vector Quantization attains quite similar results to those obtained by the CLVQ in the searching of the optimal quantizers in R d regarding the distortion error. The generation of set of random vectors in the competitive phase instead of one allows to the algorithm updating more than one quantizer at each iteration. This improvement leads to an important dismiss of consuming of computation time.

Our proposal gains in interest into problems in where the time consumption is a crucial constraint. Another important advantage of the ACLVQ is that the competitive phase is parallelizable what is not possible in the CLVQ because the strong recursive nature of this algorithm.

There are natural few ways to improve our algorithm. In a forthcoming paper we study some of them. They concern, principally, an adaptive scheme to select the number of simulations N at each competitive phase of the ACLVQ. We focus in the choice of N as a tradeoff between accuracy and computation time.

  exactly stationary quantizers. The first assertion of the lemma is clear since P is supported by O (O ⊂ R d is an open set in R d ): the functions C i (X k ) P(dω), C i (X k ) ωP(dω) are continuous on O and the first one is positive on O.

Figure 1 :

 1 Figure 1: Boxplots of quantization error for n = 10 obtained after 100 runs: N (0, I 2 ) (left) and N (0, Σ 2 ) (right).

Figure 2 :

 2 Figure 2: Boxplots of quantization error for n = 50 obtained after 100 runs: N (0, I 2 ) (left) and N (0, Σ 2 ) (right).P n N e 2 n,N s n,N 10 1 0.0170 0.3523×10 -4 50 0.0170 0.3505×10 -4 U ([0, 1] 2 )200 0.0170 0.3682×10 -4 50 1 0.0033 0.1471×10 -4 250 0.0033 0.0953×10 -4 1000 0.0033 0.0938×10 -4 10 1 0.2492 0.1134×10 -2 50 0.2494 0.1706×10 -2 Exp(Ψ 2 ) 200 0.2494 0.1577×10 -2 50 1 0.0531 0.3213×10 -3 250 0.0517 0.2856×10 -3 1000 0.0517 0.2117×10 -3

Table 2 :

 2 Parameter values used in the simulation study in dimension 2.

Table 3 :

 3 Mean and standard deviation of the quantization after 100 runs.

Table 4 :

 4 Mean of computation times (in seconds).

Appendix A. Euclidean distance matrix: a MATLAB procedure

We present a method to compute the Euclidean distance matrix of two sets of vectors in dimension greater than 2. This method is based in product of matrices and componentwise matrix multiplication which permits speed up computing in MATLAB.

Let us define Q = (Q i, j ) 1≤i≤d, 1≤j≤n the matrix of quantizers and ξ = (ξ i, k ) 1≤i≤d, 1≤k≤N the matrix of random vector with N ∈ N.

To calculate the norm of n • N vectors: D j,k := Q •, j -ξ •, k for all j and all k, we propose to decompose

The vectors

are easily computed using elementwise self-multiplication of matrix Q and ξ respectively, and summing by rows. Set D = D 2 j,k j,k it follows

The code is available at http://profesores.matcom.uh.cu/~salomon.


Christine Bernardi
email: bernardi@ann.jussieu.fr

Adel Blouza
email: adel.blouza@univ-rouen.fr

Linda El Alaoui
email: elalaoui@math.univ-paris13.fr

The rain on underground porous media Part I. Analysis of a Richards model

Richards equation models the water flow in a partially saturated underground porous medium under the surface. When it rains on the surface, boundary conditions of Signorini type must be considered on this part of the boundary. We first study this problem which results into a variational inequality. We propose a discretization by an implicit Euler's scheme in time and finite elements in space. The convergence of this discretization leads to the well-posedness of the problem.

Résumé: L'équation de Richards modélise l'écoulement d'eau dans un milieu poreux partiellement saturé souterrain. Lorsqu'il pleut, des conditions aux limites de type Signorini doivent être imposées sur la surface. Nous étudions tout d'abord ce problème qui se traduit par une inéquation variationnelle. Nous écrivons une discrétisation par schéma d'Euler implicite en temps et éléments finis en espace. La convergence de cette discrétisation permet de vérifier que le problème est bien posé.

Introduction.

The following equation

∂ t Θ(ψ) -∇ • K w Θ(ψ) ∇(ψ + z) = 0, (1.1)
models the flow of a wetting fluid, mainly water, in the underground surface, hence in an unsaturated medium, see L.A. Richards [START_REF]Richards -Capillary conduction of liquids through porous mediums[END_REF] for the introduction of this type of models.

In opposite to Darcy's or Brinkman's systems (see [START_REF]On a hierarchy of approximate models for flows of incompressible fluids through porous solid[END_REF] for all these models), this equation which is derived by combining Darcy's generalized equation with the mass conservation law is highly nonlinear: This follows from the fact that, due to the presence of air above the surface, the porous medium is only partially saturated with water. The unknown ψ is the difference between the pressure of water and the atmospherical pressure.

This equation is usually provided with Dirichlet or Neuman type boundary conditions. Indeed, Neumann boundary conditions on the underground part of the boundary are linked to the draining of water outside of the domain and Dirichlet boundary conditions on the surface are introduced to take into account the rain. However, when the porous media can no longer absorb the rainwater that falls, the upper surface ot the domain allows to exfiltration and infiltration. In other words, the upper surface is divided into a saturated and unsaturated zone. We assume that the re-infiltration process is negligible. This leads to variational inequalities of the following type :

-ψ ≥ 0, v(ψ) • n ≥ v r • n, ψ v(ψ) • n -v r • n = 0, (1.2)
where v(ψ) is the flux v(ψ) = -K w Θ(ψ) ∇(ψ + z).

(1.3)

Here, n stands for the unit normal vector to the surface and v r for a given rain fall rate. We refer to the thesis of H. Berninger [START_REF] Berninger | Domain decomposition methods for elliptic problems with jumping nonlinearites and application to the Richards equation[END_REF] for the full derivation of this model from hydrology laws and more specifically to [START_REF] Berninger | Domain decomposition methods for elliptic problems with jumping nonlinearites and application to the Richards equation[END_REF]Section 1.5] for the derivation of the boundary inequalities (1.2).

It is not so easy to give a mathematical sense to the system (1.1) -(1.2). As standard, the key argument for the analysis of problem (1.1) is to use Kirchhoff's change of unknowns. Indeed, after this transformation, the new equation fits the general framework proposed in [START_REF] Alt | Luckhaus -Quasilinear elliptic-parabolic differential equations[END_REF] (see also [START_REF] Fabrié | Gallouët -Modelling wells in porous media flows[END_REF] for the analysis of a different model). Thus, the existence and uniqueness of a solution to this equation when provided with appropriate linear initial and boundary conditions can be derived from standard arguments. In order to handle the inequality in (1.2), we again use a variational formulation. We refer to [START_REF] Alt | Visintin -On nonstationary flow through porous media[END_REF] for a first analysis of very similar systems, see also [START_REF] Brezzi | Raviart -Error estimates for the finite element solution of variational inequalities. II. Mixed methods[END_REF]. In this paper, we prove that problem (1.1)-(1.2) is well-posed when the coefficients and the data are smooth enough but without any other restriction.

The discretization of problem (1.1) has been proposed and/or studied in many papers when provided with standard boundary conditions, see [START_REF] Gabbouhy | Analyse mathématique et simulation numérique des phénomènes d'écoulement et de transport en milieux poreux non saturés. Application à la région du Gharb[END_REF], [START_REF] Radu | Knabner -Order of convergence estimates for an Euler implicit, mixed finite element discretization of Richards' equation[END_REF], [START_REF] Schneid | Radu -A priori error estimates for a mixed finite element discretization of the Richards' equation[END_REF], [START_REF] Sochala | Piperno -Mass conservative BDF-discontinuous Galerkin/explicit finite volume schemes for coupling subsurface and overland flows[END_REF] and [START_REF] Woodward | Analysis of expanded mixed finite element methods for a nonlinear parabolic equation modeling flow into variably saturated porous media[END_REF] and also [START_REF] Sochala | Ern -Numerical methods for subsurface flows and coupling with runoff[END_REF] for a more general equation, however it seems that it has not been treated when provided with the boundary inequality (1.2). We propose here a discretization of system (1.1) -(1.2), in two steps: (i) We first use the Euler's implicit scheme to build a time semi-discrete problem, where one of the nonlinear terms is treated in an explicit way for simplicity; (ii) We then construct a fully discrete problem that relies on the Galerkin method and finite elements in the spatial domain. In both cases, we prove that the corresponding variational problem is well-posed.

To conclude, we prove that the solution of this discrete problem converges to a solution of the continuous one when the discretization parameters tend to zero. This makes complete our existence result, since no restrictive condition is needed here.

The outline of the paper is as follows.

• In Section 2, we present the variational formulation of the full system and investigate its wellposedness in appropriate Sobolev spaces.

• Section 3 is devoted to the description of the time semi-discrete problem and of the fully discrete problem. We check their well-posedness.

• In Section 4, we investigate the convergence of the solution of the discrete problem towards a solution of the continuous one.

The continuous problem and its well-posedness.

Let Ω be a bounded connected open set in R d , d = 2 or 3, with a Lipschitz-continuous boundary ∂Ω, and let n denote the unit outward normal vector to Ω on ∂Ω. We assume that ∂Ω admits a partition without overlap into three parts Γ B , Γ F , and Γ G (these indices mean "bottom", "flux" and "ground", respectively) and that Γ B has a positive measure. Let also T be a positive real number.

In order to perform the Kirchhoff's change of unknowns in problem (1.1), we observe that, since the conductivity coefficient K w is positive, the mapping:

x → K(x) = x 0 K w Θ(ξ) dξ,
is one-to-one from R into itself. Thus, by setting

u = K(ψ), b(u) = Θ • K -1 (u), k • b(u) = K w • Θ • K -1 (u), ,
and thanks to an appropriate choice of the function Θ, we derive the equation (more details are given in [3, Remark 2.1] for instance)

α ∂ t u + ∂ t b(u) -∇ • ∇u + k • b(u)e z = 0 in Ω×]0, T [,
where -e z stands for the unit vector in the direction of gravity. Moreover, the Kirchhoff's change of unknowns has the further property of preseving the positivity: u is positive if and only if ψ is positive, negative if and only if ψ is negative, So, writing inequation (1.2) in terms of the unknown u is easy.

As a consequence, from now on, we work with the following system

                               α ∂ t u + ∂ t b(u) -∇ • ∇u + k • b(u)e z = 0 in Ω×]0, T [, u = u B on Γ B ×]0, T [, -∇u + k • b(u)e z • n = f F on Γ F ×]0, T [, u ≤ 0, -∇u + k • b(u)e z • n ≥ q r • n, u ∇u + k • b(u)e z + q r • n = 0 on Γ G ×]0, T [, u| t=0 = u 0 in Ω. (2.1)
The unknown is now the quantity u. The data are the Dirichlet boundary condition u B on Γ B and the initial condition u 0 on Ω, together with the boundary conditions f F and q r on the normal component of the flux, with f F corresponding to the draining of water and q r corresponding to the rain. Finally, the coefficients b and k are supposed to be known, while α is a positive constant. From now on, we assume that the function b is continuously differentiable on R, with a bounded and Lipschitz-continuous first derivative b , and that the function k • b is continuous, bounded, and uniformly Lipschitz-continuous on R.

In what follows, we use the whole scale of Sobolev spaces W m,p (Ω), with m ≥ 0 and 1 ≤ p ≤ +∞, equipped with the norm • W m,p (Ω) and seminorm | • | W m,p (Ω) , with the usual notation H m (Ω) when p = 2. As standard, the range of H 1 (Ω) by the trace operator on any part Γ of ∂Ω is denoted by H 1 2 (Γ). For any separable Banach space E equipped with the norm • E , we denote by C 0 (0, T ; E) the space of continuous functions from [0, T ] with values in E. For each integer m ≥ 0, we also introduce the space H m (0, T ; E) as the space of measurable functions on ]0, T [ with values in E such that the mappings:

v → ∂ t v E , 0 ≤ ≤ m, are square-integrable on ]0, T [.
To write a variational formulation for the problem, we introduce the time-dependent subset

V(t) = v ∈ H 1 (Ω); v| Γ B = u B (•, t) and v| Γ G ≤ 0 . (2.2)
It is readily checked that each V(t) is closed and convex, see [4, Prop. 1.5.5], when u B belongs to C 0 (0, T ; H 1 2 (Γ B )). Thus, we are led to consider the following variational problem (with obvious notation for L 2 (0, T ; V) )

Find u in L 2 (0, T ; V) such that u| t=0 = u 0 , (2.3)
and that, for a.e. t in ]0, T [,

∀v ∈ V(t), α Ω (∂ t u)(x, t)(v -u)(x, t) dx + Ω ∂ t b(u) (x, t)(v -u)(x, t) dx + Ω ∇u + k • b(u)e z (x, t) • ∇(v -u) (x, t) dx ≥ - Γ F f F (τ , t)(v -u)(τ , t) dτ - Γ G (q r • n)(τ , t)(v -u)(τ , t) dτ , (2.4)
where τ denotes the tangential coordinates on ∂Ω. The reason for this follows.

(∂ t u)(x, t)(v -u)(x, t) dx + Ω ∂ t b(u) (x, t)(v -u)(x, t) dx + Ω ∇u + k • b(u)e z (x, t) • ∇(v -u) (x, t) dx = Γ F ∪Γ G ∇u + k • b(u)e z • n(τ )(v -u)(τ , t) dτ .
To conclude, we observe that, on Γ G , either u is zero and ∇u + k • b(u)e z is smaller than -q τ • n or u is not zero and ∇u + k • b(u)e z is equal to -q r • n. All this yields (2.4).

2) Conversely, let u be any solution of (2.3) -(2.4).

• By noting that, for any function w in D(Ω), (u + w)(•, t) belongs to V(t) and taking v equal to u ± w in (2.4), we obtain the first line of (2.1) in the distribution sense.

• The second line in (2.1) follows from the definition of V(t).

• By taking v equal to u ± w for any w in D(Ω ∪ Γ F ), we also derive the third line in (2.1).

• The fact that u is nonpositive on Γ G comes from the definition of V(t). On the other hand, the previous equations imply that, for any v in V(t),

Γ G ∇u + k • b(u)e z • n(τ )(v -u)(τ , t) dτ ≥ - Γ G (q r • n)(τ , t)(v -u)(τ , t) dτ .
Taking v equal to u + w where w vanishes on Γ B and is nonpositive on Γ G yields that -∇u + k • b(u)e z • n is larger than q r • n. Finally taking v equal to zero on Γ G yields that

Γ G ∇u + k • b(u)e z + q r • n(τ )u(τ , t) dτ ≤ 0,
and since the two quantities u and ∇u + k • b(u)e z + q r are nonpositive on Γ G , their product is zero.

• And finally the last line of (2.1) is written in (2.3).

Proving that problem (2.3) -(2.4) is well-posed is not at all obvious. We begin with the simpler result, i.e., the uniqueness of the solution. For brevity, we set:

X = L 2 (0, T ; V) ∩ H 1 (0, T ; L 2 (Ω)).
(2.5)

We also refer to [START_REF] Lions | Magenes -Problèmes aux limites non homogènes et applications[END_REF]Chap. 1,Th. 11.7] for the definition of the space H

1 2 00 (Γ B ).
Proposition 2.2. For any data u B , f F , q r and u 0 satisfying

u B ∈ H 1 (0, T ; H 1 2 00 (Γ B )), f F ∈ L 2 (0, T ; L 2 (Γ F )), q r ∈ L 2 (0, T ; L 2 (Γ G ) d ), u 0 ∈ H 1 (Ω), (2.6
)

problem (2.3) -(2.
4) has at most a solution in X.

Proof: Let u 1 and u 2 be two solutions of problem (2.3) -(2.4). Thus, the function u = u 1 -u 2 vanishes on Γ B and at t = 0. Taking v equal to u 2 in the problem satisfied by u 1 and equal to u 1 in the problem satisfied by u 2 and subtracting the second problem from the first one, we obtain

α Ω (∂ t u)(x, t)u(x, t) dx + Ω ∂ t b(u 1 ) -∂ t b(u 2 ) (x, t)u(x, t) dx + Ω (∇u) 2 (x, t) dx + Ω k • b(u 1 ) -k • b(u 2 ) (x, t)e z • ∇u) (x, t) dx ≤ 0.
(2.7)

We integrate this equation with respect to t, use the decomposition

Ω ∂ t b(u 1 ) -∂ t b(u 2 ) (x, t)u(x, t) dx = Ω b (u 1 )(x, t) ∂ t u) (x, t)u(x, t) dx + Ω b (u 1 ) -b (u 2 ) (x, t)∂ t u 2 (x, t)u(x, t) dx.
Next, the nonnegativity and Lipschitz-continuity of b , combined with the interpolation inequality (see [11, Chap. 1, Prop. 2.3]) and the Poincaré-Friedrichs inequality

u L 4 (Ω) ≤ u 1-d 4 L 2 (Ω) (c|u| H 1 (Ω) ) d 4 ≤ c (1 - d 4 ) u L 2 (Ω) + d 4 |u| H 1 (Ω) ,
lead to handle the second term. The Lipschitz-continuity of k • b together with a Young's inequality is needed to handle the fourth term. All this gives

α u(•, t) 2 L 2 (Ω) + 1 2 t 0 |u(•, s)| 2 H 1 (Ω) ds ≤ c(u 2 ) t 0 u(•, s) 2 L 2 (Ω) ds.
Thus, applying Grönwall's lemma yields that u is zero, whence the uniqueness result.

Proving the existence is much more complex. We begin with a basic result.

Lemma 2.3. If the function u B belongs to C 0 (0, T ; H 1 2 00 (Γ B )), for all t in [0, T ], the convex set V(t) is not empty.
Proof: Denoting by u B (•, t) the extension by zero of u B (•, t) to ∂Ω, we observe that any lifting of u B (•, t) in H 1 (Ω) belongs to V(t), whence the desired result.

In a first step, we consider the linear problem, for any datum F in L 2 (0, T ; L 2 (Ω)):

Find u in L 2 (0, T ; V) satisfying (2.
3) and such that, for a.e.

t in ]0, T [, ∀v ∈ V(t), α Ω (∂ t u)(x, t)(v -u)(x, t) dx + Ω ∇u (x, t) • ∇(v -u) (x, t) dx ≥ - Ω F (x, t)(v -u)(x, t) dx - Γ F f F (τ , t)(v -u)(τ , t) dτ - Γ G (q r • n)(τ , t)(v -u)(τ , t) dτ .
(2.8) However a weaker formulation of this problem can be derived by integrating in time.

It reads

Find u in L 2 (0, T ; V) satisfying (2.3) and such that ∀v ∈ X, α T 0 Ω (∂ t u)(x, t)(v -u)(x, t) dxdt + T 0 Ω ∇u (x, t) • ∇(v -u) (x, t) dxdt ≥ - T 0 Ω F (x, t)(v -u)(x, t) dxdt - T 0 Γ F f F (τ , t)(v -u)(τ , t) dτ dt - T 0 Γ G (q r • n)(τ , t)(v -u)(τ , t) dτ dt.
(2.9)

We recall in the next lemma the properties of this problem which are standard.

Lemma 2.4. Assume that the data u B , f F , q r and u 0 satisfy (2.6). Then, for any F in L 2 (0, T ; L 2 (Ω)), problem (2.3) -(2.9) has a unique solution u in L 2 (0, T ; V).

Proof: It follows from Lemma 2.3 and the further assumption on u B that X is a non empty closed convex set. We also consider a lifting u B of the extension by zero of u B to ∂Ω in H 1 (0, T ; H 1 (Ω)). Then it is readily checked that u -u B is the solution of a problem which satisfies all the assumptions in [10, Chap. 6, Th. 2.2], whence the existence and uniqueness result.

Any solution of (2.3) -(2.8) is a solution of (2.3) -(2.9), but the converse property is not obvious in the general case, see [START_REF] Glowinski | Trémolières -Analyse numérique des inéquations variationnelles. 2. Applications aux phénomènes stationnaires et d'évolution[END_REF]Chap. 6]. However, in our specific case, it is readily checked by a density argument that inequation (2.9) is satisfied for any v in L 2 (0, T ; V), so that problems (2.3) -(2.8) and (2.3) -(2.9) are fully equivalent.

To go further, we assume that the following compatibility condition holds: for a.e. x ∈ Γ B , u 0 (x) = u B (x, 0) and for a.e. x ∈ Γ G , u 0 (x) ≤ 0, (2.10) and we introduce a lifting u * B of an extension of u B to ∂Ω which belongs to H 1 (0, T ; V) and satisfies for a.e. x ∈ Ω, u * B (x, 0) = u 0 (x), (2.11) together with the stability property

u * B H 1 (0,T ;H 1 (Ω)) ≤ c u B H 1 (0,T ;H 1 2 00 (Γ B ))
.

(2.12)

Then it is readily checked that u is a solution of problem (2.3) -(2.4) if and only if the function

u * = u -u * B is a solution of Find u * in L 2 (0, T ; V 0 ) with ∂ t u in L 2 (0, T ; L 2 (Ω)) such that u * | t=0 = 0, (2.13)
and that, for a.e. t in ]0, T [,

∀v ∈ V 0 , α Ω (∂ t u * )(x, t)(v -u * )(x, t) dx + Ω ∂ t b * (u * ) (x, t)(v -u * )(x, t) dx + Ω ∇u * + k • b * (u * )e z (x, t) • ∇(v -u * ) (x, t) dx ≥ - Ω F B (x, t)(v -u * )(x, t) dx - Γ F f F (τ , t)(v -u * )(τ , t) dτ - Γ G (q r • n)(τ , t)(v -u * )(τ , t) dτ , (2.14)
with obvious definition of the subset V 0 :

V 0 = v ∈ H 1 (Ω); v| Γ B = 0 and v| Γ G ≤ 0 , (2.15)
and where the new application b * is defined by b

* (u * ) = b(u * + u * B ). The datum F B is defined by, for a.e. t in ]0, T [, Ω F B (x, t)v(x) dx = α Ω (∂ t u * B )(x, t)v(x) dx + Ω (∇u * B )(x, t) • (∇v)(x) dx, (2.16)
and clearly belongs to L 2 (0, T ; W ), where W is the smallest linear space containing V 0 , namely

W = v ∈ H 1 (Ω); v| Γ B = 0 . (2.17)
It can be noted that the existence result stated in Lemma 2.4 is still valid for any F in L 2 (0, T ; W ).

We denote by T the operator which associates with any pair (F, D), with F in L 2 (0, T ; W ) and data D = (0, f F , q r , 0) satisfying (2.6), the solution

u of problem (2.3) - (2.8). It follows from (2.13) -(2.14) that u * satisfies u * -T (F B + F (u * ), D) = 0, (2.18)
where the quantity F (u) is defined by duality, for a.e. t in [0, T ]:

F (u), v = Ω ∂ t b(u) (x, t)v(x) dx + Ω k • b(u)(x, t)e z • (∇v)(x) dx. (2.19)
We first prove some further properties of the operator T .

Lemma 2.5. The operator T is continuous from

L 2 (0, T ; W ) × L 2 (0, T ; L 2 (Γ F )) × L 2 (0, T ; L 2 (Γ G ) d
) into the space L 2 (0, T ; V 0 ). Moreover the following estimate holds

T 0 |T (F, f F , q r )(•, t)| 2 H 1 (Ω) dt 1 2 ≤ F L 2 (0,T ;W ) + c f F L 2 (0,T :L 2 (Γ F )) + c q r L 2 (0,T ;L 2 (Γ G ) d ) .
(2.20)

Proof: We set: u = T (F, f F , q r ) and we only prove the estimate (indeed, it is readily checked that it implies the continuity property). We take v equal to u 2 in problem (2.8). This obviously gives

α 2 Ω (∂ t u 2 )(x, t) dx + |u(•, t)| 2 H 1 (Ω) ≤ F (•, t) W + c f F (•, t) L 2 (Γ F ) + c q r (•, t) L 2 (Γ G ) d |u(•, t)| H 1 (Ω) ,
where c is the norm of the trace operator. Thus, integrating with respect to t gives estimate (2.20).

Lemma 2.6. The operator T is continuous from

L 2 (0, T ; L 2 (Ω)) × H 1 (0, T ; L 2 (Γ F )) × H 1 (0, T ; L 2 (Γ G ) d
) into the space H 1 (0, T ; L 2 (Ω). Moreover the following estimate holds for any positive ε

α ∂ t T (F, f F , q r ) L 2 (0,T ;L 2 (Ω)) ≤ (1 + ε) F L 2 (0,T ;L 2 (Ω)) + c f F H 1 (0,T :L 2 (Γ F )) + c q r H 1 (0,T ;L 2 (Γ G ) d ) .
(2.21)

Proof: The continuity property of T is proved in [10, Chap. 6, Th. 2.1]. Next, setting: u = T (F, f F , q r ), we take v equal to u -η ∂ t u in (2.8) for a positive η. Indeed,

• since u vanishes on Γ B , so does ∂ t u;

• since u is nonpositive on Γ G and u(x, t -η) which is very near of u(x, t) -η ∂ t u(x, t) is also nonpositive, there exists a η > 0 such that u -η ∂ t u belongs to V 0 . This yields

α ∂ t u 2 L 2 (Ω) + 1 2 ∂ t |u| 2 H 1 (Ω) ≤ F L 2 (Ω) ∂ t u L 2 (Ω) - Γ F f F (τ , t)∂ t u(τ , t) dτ - Γ G (q r • n)(τ , t)∂ t u(τ , t) dτ .
To bound the first term, we use Young's inequality

F L 2 (Ω) ∂ t u L 2 (Ω) ≤ α 2 ∂ t u 2 L 2 (Ω) + 1 2α F 2 L 2 (Ω) .
To handle the last integrals, we integrate by parts with respect to t. For instance, we have, for any ε > 0,

t 0 Γ F f F (τ , s)∂ t u(τ , s) dτ ds = Γ F f F (τ , t)u(τ , t) dτ ds - t 0 Γ F ∂ t f F (τ , s) u(τ , s) dτ ds ≤ 1 4 |u(•, t)| 2 H 1 (Ω) + c f F (•, t) 2 L 2 (Γ F ) + c ∂ t f F 2 L 2 (0,t;L 2 (Γ )) + ε u 2 L 2 (0,t;H 1 (Ω)) .
Thus, the desired estimate follows by combining all this and using (2.20).

We are thus in a position to prove a first existence result.

Theorem 2.7. Assume that the coefficient α satisfies

1 α b L ∞ (R) < 1. ( 2

.22)

For any data u B , f F , q r and u 0 satisfying

u B ∈ H 1 (0, T ; H 1 2 00 (Γ B )), f F ∈ H 1 (0, T ; L 2 (Γ F )), q r ∈ H 1 (0, T ; L 2 (Γ G ) d ), u 0 ∈ H 1 (Ω).
(2.23) and (2.10), problem (2.3) -(2.4) has at least a solution in X.

Proof: We proceed in several steps. 1) Let X 0 be the space of functions of X vanishing in t = 0. We provide it with the norm:

v X 0 = ∂ t v L 2 (0,T ;L 2 (Ω)) .
It follows from the previous Lemma 2.6 that

T (F B + F (u * ), D) X 0 ≤ 1 + ε α F (u * ) L 2 (0,T ;L 2 (Ω)) + c(D),
where the constant c(D) only depends on the data u B , f F and q r . Due to te boundedness of b and k • b (see (2.19) for the definition of F (u * )), we have

T (F B + F (u * ), D) X 0 ≤ 1 + ε α b L ∞ (R) u * X 0 + c (D).
Thus, due to Assumption (2.22), the application: u * → T (F B + F (u * ), D) maps the ball in X 0 with radius R into itself for all R such that, for an appropriate ε,

1 - 1 + ε α b L ∞ (R) R > c (D). (2.24)
2) Since X 0 is separable, there exists an increasing sequence of finite-dimensional spaces X n which is dense in X 0 . If Π n denotes the orthogonal projection operator (for the scalar product associated with the norm of X 0 ) onto X n , the mapping:

u → Π n T (F B + F (u), D)
is continuous from X n into itself. The same arguments as previously yield that it sends the ball of X n with radius R into itself for all R satisfying (2.24). Thus, applying Brouwer's fixed point theorem (see [9, Chap. IV, Th. 1.1] for instance) yields that this mapping admits a fixed point in this same ball, namely that there exists u n in X n satisfying the equation u n = Π n T (F B + F (u n ), D). Moreover, it follows from Lemma 2.5 that this sequence is also bounded in L 2 (0, T ; H 1 (Ω)).

3) The function u n thus satisfies

∀v ∈ X n , α Ω (∂ t u n )(x, t)(v -u n )(x, t) dx + Ω ∂ t b * (u n ) (x, t)(v -u n )(x, t) dx + Ω ∇u n + k • b * (u n )e z (x, t) • ∇(v -u n ) (x, t) dx ≥ - Ω F B (x, t)(v -u n )(x, t) dx - Γ F f F (τ , t)(v -u n )(τ , t) dτ - Γ G (q r • n)(τ , t)(v -u n )(τ , t) dτ .
(2.25)

Moreover, due to the boundedness properties of the sequence (u n ) n , there exists a subsequence still denoted by (u n ) n for simplicity which converges to a function u * of X 0 weakly in X and strongly in L 2 (0, T ; L 2 (Ω)). Next, we observe that, for a fixed v in X n :

• The convergence of all terms in the right-hand side follows from the weak convergence in L 2 (0, T ; W);

• The convergence of the first term is derived by writing the expansion

Ω (∂ t u n )(x, t)(v -u n )(x, t) dx = Ω (∂ t u * )(x, t)(v -u * )(x, t) dx + Ω ∂ t (u n -u * )(x, t)(v -u * )(x, t) dx + Ω (∂ t u n )(x, t)(u * -u n )(x, t) dx,
and checking that the two last terms converge;

• The convergence of the term Ω (∇u n )(x, t) • ∇(v -u n ) (x, t) dx is obtained by using the weak lower semi-continuity of the norm |u n | H 1 (Ω) . Moreover, the convergence of the nonlinear terms follows from the expansions

Ω ∂ t b * (u n ) (x, t)(v -u n )(x, t) dx = Ω ∂ t b * (u * ) (x, t)(v -u * )(x, t) dx + Ω ∂ t b * (u n ) -∂ t b * (u * ) (x, t)(v -u * )(x, t) dx + Ω ∂ t b * (u n ) (x, t)(u * -u n )(x, t) dx,
and Corollary 2.8. Assume that the coefficient b satisfies (2.26), and that the coefficient α is such that

Ω k • b * (u n )(x, t)e z • ∇(v -u n ) (x, t) dx = Ω k • b * (u * )(x, t)e z • ∇(v -u * ) (x, t) dx + Ω k • b * (u * )(x, t)e z • ∇(u * -u n ) (x, t) dx + Ω k • b * (u n ) -k • b * (u * ) (x, t)e z • ∇(v -u n ) (x, t) dx,
1 α + b 0 b -b 0 L ∞ (R) < 1. ( 2

.27)

For any data u B , f F , q r and u 0 satisfying (2.23) and (2. (2.28) Under this condition, problem (2.3) -(2.4) has a solution even for α = 0. We refer to [START_REF] Alt | Visintin -On nonstationary flow through porous media[END_REF] for another proof of this result for a similar problem.

The discrete problems.

We present successively the time semi-discrete problem constructed from the backward Euler's scheme, next a finite element discretization of this problem relying on standard, conforming, finite element spaces.

The time semi-discrete problem.

Since we intend to work with non uniform time steps, we introduce a partition of the interval

[0, T ] into subintervals [t n-1 , t n ], 1 ≤ n ≤ N , such that 0 = t 0 < t 1 < • • • < t N = T .
We denote by τ n the time step t n -t n-1 , by τ the N -tuple (τ 1 , . . . , τ N ) and by |τ | the maximum of the τ n , 1 ≤ n ≤ N .

As already hinted in Section 1, the time discretization mainly relies on a backward Euler's scheme, where the nonlinear term k•b(u) is treated in an explicit way for simplicity. Thus, the semi-discrete problem reads

Find (u n ) 0≤n≤N in N n=0 V(t n ) such that u 0 = u 0 in Ω, (3.1)
and, for 1

≤ n ≤ N , ∀v ∈ V(t n ), α Ω u n -u n-1 τ n (x)(v -u n )(x) dx + Ω b(u n ) -b(u n-1 ) τ n (x)(v -u n )(x) dx + Ω ∇u n + k • b(u n-1 ) (x)e z • ∇(v -u n )x) dx ≥ - Γ F f F (τ , t n )(v -u n )(τ ) dτ - Γ G (q r • n)(τ , t n )(v -u n )(τ ) dτ . (3.2)
It can be noted that this problem makes sense when both f F and q r are continuous in time. Proving its well-posedness relies on rather different arguments as previously.

Theorem 3.1. For any data u B , f F , q r and u 0 satisfying Proof: We proceed by induction on n. Since u 0 is given by (3.1), we assume that u n-1 is known and we consider problem (3.2) for a fixed n, called (3.2) n , that can equivalently be written

u B ∈ H 1 (0, T ; H 1 2 00 (Γ B )), f F ∈ C 0 (0, T ; L 2 (Γ F )), q r ∈ C 0 (0, T ; L 2 (Γ G ) d ), u 0 ∈ H 1 (Ω), ( 3
∀v ∈ V(t n ), Ω α u n + b(u n ) (x)(v -u n )(x) dx + τ n Ω ∇u n (x) • ∇(v -u n )x) dx ≥ Ω α u n-1 + b(u n-1 ) (x)(v -u n )(x) dx -τ n Ω k • b(u n-1 )(x)e z • ∇(v -u n )(x) dx -τ n Γ F f F (τ , t n )(v -u n )(τ ) dτ -τ n Γ G (q r • n)(τ , t n )(v -u n )(τ ) dτ .
Let us now set:

ϕ(z) = z 0 α ζ + b(ζ) dζ, Φ(v) = Ω ϕ(v(x)) dx.
It is readily checked that, since b is nonnegative, both ϕ and Φ are convex and, moreover, that DΦ(u)

• (v -u n ) = Ω α u + b(u) (x)(v -u n )(x) dx.
Thus, when taking

a(u, v) = Ω ∇u(x) • ∇v(x) dx, (v) = Ω α u n-1 + b(u n-1 ) (x)v(x) dx -τ n Ω k • b(u n-1 )(x)e z • ∇v(x) dx -τ n Γ F f F (τ , t n )v(τ ) dτ -τ n Γ G (q r • n)(τ , t n )v(τ ) dτ , problem (3.
2) n can also be written

∀v ∈ V(t n ), DΦ(u n ) • (v -u n ) + a(u n , v -u n ) -(v -u n ) ≥ 0.
We now set:

Ψ(v) = Φ(v) + J(v), with J(v) = 1 2 a(v, v) -(v). Problem (3.2) n can finally be written ∀v ∈ V(t n ), DΨ(u n ) • (v -u n ) ≥ 0, or ∀v ∈ V(t n ), Ψ(u n ) ≤ Ψ(v).
So it is equivalent to the minimization of a convex functional on the convex set V(t n ), hence admits a unique solution. This concludes the proof.

It can be noted that, in contrast with the continuous problem, the existence of a solution for the semi-discrete problem (3.1) -(3.2) does not require any limitation on α.

The fully discrete problem.

From now on, we assume that Ω is a polygon (d = 2) or a polyhedron (d = 3). Let (T h ) h be a regular family of triangulations of Ω (by triangles or tetrahedra), in the sense that, for each h:

• Ω is the union of all elements of T h ; • The intersection of two different elements of T h , if not empty, is a vertex or a whole edge or a whole face of both of them; • The ratio of the diameter h K of any element K of T h to the diameter of its inscribed circle or sphere is smaller than a constant σ independent of h. As usual, h stands for the maximum of the diameters h K , K ∈ T h . We make the further and non restrictive assumption that Γ B , Γ F and Γ G are the union of whole edges (d = 2) or whole faces (d = 3) of elements of T h . From now on, c, c , . . . stand for generic constants that may vary from line to line but are always independent of τ and h.

We now introduce the finite element space

V h = v h ∈ L 2 (Ω); ∀K ∈ T h , v h | K ∈ P 1 (K) , (3.4)
where P 1 (K) is the space of restrictions to K of affine functions on R d . Let I h denote the Lagrange interpolation operator at all the vertices of elements of T h with values in V h and i B h the corresponding interpolation operator on Γ B . Assuming that u B is continuous where needed, we then define for each n, 0 ≤ n ≤ N , the subset of V h :

V h (t n ) = v h ∈ V h ; v h | Γ B = i B h u B (•, t n ) and v h | Γ G ≤ 0 . (3.5)
We are thus in a position to write the discrete problem, constructed from problem (3.1) -(3.2) by the Galerkin method,

Find (u n h ) 0≤n≤N in N n=0 V h (t n ) such that u 0 h = I h u 0 in Ω, (3.6)
and, for 1 ≤ n ≤ N ,

∀v h ∈ V h (t n ), α Ω u n h -u n-1 h τ n (x)(v h -u n h )(x) dx + Ω b(u n h ) -b(u n-1 h ) τ n (x)(v h -u n h )(x) dx + Ω ∇u n h + k • b(u n-1 h ) (x)e z • ∇(v h -u n h )(x) dx ≥ - Γ F f F (τ , t n )(v h -u n h )(τ ) dτ - Γ G (q r • n)(τ , t n )(v h -u n h )(τ ) dτ . (3.7)
The proof of the next theorem is exactly the same as the proof of Theorem 3.1, so we omit it. Theorem 3.2. For any data u B , f F , q r and u 0 satisfying (3.3), (2.10) and

u B ∈ C 0 (Γ B × [0, T ]), u 0 ∈ C 0 (Ω), (3.8)
for any nonnegative coefficient α, problem (3.6) -(3.7) has a unique solution.

There also, the existence result is unconditional.

A convergence result.

The aim of this section is to prove a convergence result for the soutions (u n h ) 0≤n≤N of problem (3.6) -(3.7) when |τ | and h tend to zero. In order to do that and as in Section 2, we use the lifting u * B of u B which satisfies (2.11) and (2.12) and we assume moreover that it is continuous on Ω

× [0, T ]. Indeed, if (u n h ) 0≤n≤N is a solution of (3.6) -(3.7), the family (u * n h ) 0≤n≤N , with u * n h = u n h -I h u * B (t n ), is a solution of Find (u * n h ) 0≤n≤N in V N +1 h0 such that u * 0 h = 0 in Ω, (4.1)
and, for 1 ≤ n ≤ N ,

∀v h ∈ V h0 , α Ω u * n h -u * n-1 h τ n (x)(v h -u * n h )(x) dx + Ω b * n (u * n h ) -b * n-1 (u * n-1 h ) τ n (x)(v h -u * n h )(x) dx + Ω ∇u * n h + k • b * n-1 (u * n-1 h ) (x)e z • ∇(v h -u * n h )(x) dx ≥ - Ω F Bh (x, t n )(v h -u * n h ) dx - Γ F f F (τ , t n )(v h -u * n h )(τ ) dτ - Γ G (q r • n)(τ , t n )(v h -u * n h )(τ ) dτ , (4.2)
where the convex set V h0 and the function F Bh are defined, in analogy with (2.15) and (2.16), by

V h0 = V h ∩ V 0 , (4.3)
and

Ω F Bh (x, t)v(x) dx = α Ω (∂ t I h u * B )(x, t)v(x) dx + Ω (∇I h u * B )(x, t) • (∇v)(x) dx, (4.4)
while each function b * n is given by b

* n (ξ) = b(ξ + I h u * B (•, t n ))
. We now investigate the boundedness of the sequence (u * n h ) 0≤n≤N in appropriate norms, we need a preliminary lemma for that. Lemma 4.1. For each part Γ of ∂Ω which is the union of whole edges (d = 2) or whole faces (d = 3) of elements of T h , the following inequality holds for all functions w h in V h ,

w h H -1 2 (Γ) ≤ c w h L 2 (Ω) . (4.5)
Proof: It relies on standard arguments. We have

w h H -1 2 (Γ) = sup z∈H 1 2 (Γ) Γ z(τ )w h (τ ) dτ z H 1 2 (Γ)
.

Let e be any edge or face of an element K of T h which is contained in Γ. Denoting by K the reference triangle or tetrahedron, we have, with obvious notation for ê, ŵ, ẑ,

e z(τ )w h (τ ) dτ ≤ c h d-1 e ê ẑ(τ ) ŵh (τ ) dτ ≤ c h d-1 K ẑ L 2 (ê) ŵh L 2 (ê) .
By using the equivalence of norms on P 1 ( K) and an appropriate stable lifting π which maps traces on ê into functions of K vanishing at the vertex of K which does not belong to Γ, we derive

e z(τ )w h (τ ) dτ ≤ c h d-1 K |πẑ| H 1 ( K) ŵh L 2 ( K) ≤ c h d-1 K h 1-d 2 K |πz| H 1 (K) h -d 2 K w h L 2 (K) ,
there also with an obvious definition of π. We conclude by summing this last inequality on the e and using a Cauchy-Schwarz inequality and the stability of π:

Γ z(τ )w h (τ ) dτ ≤ c z H 1 2 (Γ) w h L 2 (Ω) ,
whence the desired result.

Lemma 4.2. For any data u B , f F , q r and u 0 satisfying

u B ∈ H 1 (0, T ; H 1 2 00 (Γ B )), f F ∈ C 0 (0, T ; H 1 2 (Γ F )), q r ∈ C 0 (0, T ; H 1 2 (Γ G ) d ), u 0 ∈ H 1 (Ω), (4.6)
and (2.10), the sequence (u * n h ) 0≤n≤N satisfies the following bound, for 1 ≤ n ≤ N ,

α n m=1 τ m u * m h -u * m-1 h τ m 2 L 2 (Ω) + |u * n h | 2 H 1 (Ω) ≤ c 1 + I h u * B 2 H 1 (0,T ;H 1 (Ω)) + f F 2 C 0 (0,T ;H 1 2 (Γ F )) + q r 2 C 0 (0,T ;H 1 2 (Γ G ) d )
.

(4.7)

Proof: Taking v equal to u * n-1 h in (4.2) yields α τ n u * n h -u * n-1 h τ n 2 L 2 (Ω) + Ω ∇u * n h (x) • ∇(u * n h -u * n-1 h )(x) dx ≤ - Ω b * n (u * n h ) -b * n-1 (u * n-1 h ) τ n (x)(u * n h -u * n-1 h )(x) dx - Ω k • b * n-1 (u * n-1 h )(x)e z • ∇(u * n h -u * n-1 h )(x) dx + G, u * n h -u * n-1 h ,
where the data depending quantity G is defined by

G, v = - Ω F Bh (x, t n )v(x) dx - Γ F f F (τ , t n )v(τ ) dτ - Γ G (q r • n)(τ , t n )v(τ ) dτ .
To handle the second term, we use the identity

Ω ∇u * n h • ∇(u * n h -u * n-1 h )(x) dx = 1 2 |u * n h | 2 H 1 (Ω) + |u * n h -u * n-1 h | 2 H 1 (Ω) -|u * n-1 h | 2 H 1 (Ω) .
To handle the third term, we write the expansion

Ω b * n (u * n h ) -b * n-1 (u * n-1 h ) τ n (x)(u * n h -u * n-1 h )(x) dx = Ω b(u * n h + I h u * B (t n )) -b(u * n-1 h + I h u * B (t n )) τ n (x)(u * n h -u * n-1 h )(x) dx + Ω b(u * n-1 h + I h u * B (t n )) -b(u * n-1 h + I h u * B (t n-1 )) τ n (x)(u * n h -u * n-1 h )(x) dx,
and, by using the nonnegativity of b together with the Lipschitz-continuity of b, we derive

Ω b * n (u * n h ) -b * n-1 (u * n-1 h ) τ n (x)(u * n h -u * n-1 h )(x) dx ≤ α 4 τ n u * n h -u * n-1 h τ n 2 L 2 (Ω) + 1 α τ n I h u * B (t n ) -I h u * B (t n-1 ) τ n 2 L 2 (Ω) .
Finally, evaluating the last term is an easy consequence of Lemma 4.1:

G, u * n h -u * n-1 h ≤ α 4 τ n u * n h -u * n-1 h τ n 2 L 2 (Ω) + cτ n F Bh (•, t n ) 2 L 2 (Ω) + f F (•, t n ) 2 H 1 2 (Γ F ) + q r (•, t n ) 2 H 1 2 (Γ G ) d
. By combining all this, we obtain

α 2 τ n u * n h -u * n-1 h τ n 2 L 2 (Ω) + 1 2 |u * n h | 2 H 1 (Ω) ≤ 1 2 |u * n-1 h | 2 H 1 (Ω) + c τ n F Bh (•, t n ) 2 L 2 (Ω) + f F (•, t n ) 2 H 1 2 (Γ F ) + q r (•, t n ) 2 H 1 2 (Γ G ) d - Ω k • b * n-1 (u * n-1 h )(x)e z • ∇(u * n h -u * n-1 h )(x) dx.
We sum this inequality on n. To handle the last term, we observe that

- n m=1 Ω k • b * m-1 (u * m-1 h )(x)e z • ∇(u * m h -u * m-1 h )(x) dx = - Ω k • b * n-1 (u * n-1 h )(x)e z • ∇u * n h (x) dx + n-1 m=1 Ω k • b * m (u * m h ) -k • b * m-1 (u * m-1 h ) (x)e z • ∇u * m h (x) dx,
whence, thanks to the boundedness of k and Lipschitz continuity of k • b,

- n m=1 Ω k • b * m-1 (u m-1 * h )(x)e z • ∇(u * m h -u * m-1 h )(x) dx = c + 1 4 |u * n h | 2 H 1 (Ω) + α 4 n-1 m=1 τ m u * m h -u * m-1 h τ m 2 L 2 (Ω) + c n-1 m=1 τ m I h u * B (t m ) -I h u * B (t m -1) τ m 2 L 2 (Ω) + c n-1 m=1 τ m |u * m | 2 H 1 (Ω) .
We conclude by using the discrete Grönwall's lemma, see [8, Chap. V, Lemma 2.4] for instance.

Let us now introduce the function u * hτ which is affine on each interval [t n-1 , t n ], 1 ≤ n ≤ N , and equal to u * n h at time t n , 0 ≤ n ≤ N . When the data u B , f F , q r and u 0 satisfy

u B ∈ H 1 (0, T ; H s (Γ B )), f F ∈ C 0 (0, T ; H 1 2 (Γ F )), q r ∈ C 0 (0, T ; H 1 2 (Γ G ) d ), u 0 ∈ H s+ 1 2 (Ω), (4.8)
for some s > d-1 2 (in order to ensure the stability the operator I h ), it follows from Lemma 4.2 that this function belongs to the set X 0 = L 2 (0, T ; V 0 ) ∩ H 1 (0, T ; L 2 (Ω)), see (2.5) and (2.14); more precisely it satisfies

u * hτ L 2 (0,T ;H 1 (Ω))∩H 1 (0,T ;L 2 (Ω)) ≤ c(u B , f F , q r ), (4.9)
where the constant c(u B , f F , q r ) only depends on the data. Thus, we are in a position to derive the next result.

Theorem 4.3. For any data u B , f F , q r and u 0 satisfying (4.8) and (2.10), and for any positive coefficient α, problem (2.3) -(2.4) has at least a solution in X.

Proof: Thanks to (4.9), the family of functions u * hτ is bounded in X 0 independently of h and τ . Thus, there exist a sequence (T hk ) k of triangulations T h and a sequence (τ k ) k of pameters τ such that, with obvious notation, the sequence (u * k ) k converges to a function u * of X 0 weakly in L 2 (0, T ; H 1 (Ω) ∩ H 1 (0, T ; L 2 (Ω)) and strongly in L 2 (0, T ; L 2 (Ω)). We now intend to prove that u * is a solution of problem (2.13) -(2.14). Since it obviously satisfies (2.13), we now investigate the convergence of all terms in (4.2). For clarity, we keep the notation u * n h for u * k (t n ). 1) The convergence of the first term follows from the expansion

α Ω u * n h -u * n-1 h τ n (x)(v h -u * n h )(x) dx = α Ω ∂ t u * (x, t n )(v h -u * )(x, t n ) dx + α Ω ∂ t (u * k -u * ) (x, t n )(v h -u * )(x, t n ) dx + α Ω ∂ t u * k (x, t n )(u * -u * n h )(x, t n ) dx. 20
2) To prove the convergence of the term

Ω b * n (u * n h ) -b * n-1 (u * n-1 h ) τ n (x)(v h -u * n h )(x) dx,
we use a rather complex expansion that we skip for brevity, combined with the dominated convergence theorem of Lebesgue: Indeed, since (u * k ) k converges to a function u * in L 2 (0, T ; L 2 (Ω)), it converges almost every where in Ω×]0, T [, so that (b (u * k )) k also converges a.e. to b (u * ); thus, since b is bounded, (b (u * k )) k also converges to b (u * ) in L 2 (0, T ; L 2 (Ω)).

3) The convergence of the term Ω ∇u * n h (x, t n )e z • ∇(v h -u * n h )(x, t n )) dx is a consequence of the weak lower semi-continuity of the norm. 4) The convergence of the term

Ω k • b * n-1 (u n-1 * h )(x)e z • ∇(v h -u * n h )(x) dx is easily derived from the expansion Ω k • b * n-1 (u * n-1 h )(x)e z • ∇(v h -u * n h )(x) dx = Ω k • b * (u * )(x, t n )e z • ∇(v h -u * )(x, t n ) dx + Ω (k • b * n-1 -k • b * )(u * )(x, t n )e z • ∇(v h -u * )(x, t n ) dx + Ω k • b * n-1 (u * )(x, t n )e z • ∇(u * -u * n h )(x) dx + Ω k • b * n-1 (u * n-1 h ) -k • b * n-1 (u * ) (x)e z • ∇(v h -u * n h )(x, t n ) dx,
and there also from the dominated convergence theorem of Lebesgue.

5) The convergence of all terms in the right-hand side of (4.2) is obviously derived from the weak convergence of the sequence (u * k ) k . Finally, using the density of the union of the V h0 in V 0 , we derive that u * is a solution of problem (2.13) -(2.14). Thus, the function u = u * + u * B is a solution of problem (2.3) -(2.4).

Even if this requires a slightly different regularity of the data, Theorem 4.3 combined with Proposition 2.2 yields that, for any positive coefficient α, problem (2.3) -(2.4) is well-posed in X. This of course is a great improvement of the results in Section 2 and leads to think that the discretization proposed in Section 3 is rather efficient. We shall check this in the second part of this work.

  combined with the Lipschitz-continuity of b and k • b. Finally, using the density of the sequence of X n in X 0 yields that u * is a solution of problem (2.13)-(2.14). Thus, u is a solution of problem (2.3) -(2.4). Condition (2.22) is rather restrictive since, in pratical situations, α is small. However, this condition can be relaxed when the coefficient b satisfies, for a positive constant b 0 , ∀ξ ∈ R, b (ξ) ≥ b 0 . (2.26) Indeed all the previous arguments are still valid when replacing α by α + b 0 and the coefficient b(ξ) by b(ξ) -b 0 ξ.

  10), problem (2.3) -(2.4) has at least a solution in X. Assume that the coefficient b satisfies min ξ∈R b (ξ) > 0, max ξ∈R b (ξ) < 2 min ξ∈R b (ξ).

. 3 )

 3 and (2.10), for any nonnegative coefficient α, problem (3.1) -(3.2) has a unique solution in N n=0 V(t n ).