Table of Contents Author Guidelines Submit a Manuscript
Advances in Numerical Analysis
Volume 2011 (2011), Article ID 164581, 24 pages
Research Article

Analysis and Finite Element Approximation of a Nonlinear Stationary Stokes Problem Arising in Glaciology

1Department of Mathematics and Computer Science, Free University of Berlin, 14195 Berlin, Germany
2Chair of Numerical Analysis and Simulation, EPFL, 1015 Lausanne, Switzerland

Received 13 September 2011; Accepted 25 October 2011

Academic Editor: Michele Benzi

Copyright © 2011 Guillaume Jouvet and Jacques Rappaz. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


The aim of this paper is to study a nonlinear stationary Stokes problem with mixed boundary conditions that describes the ice velocity and pressure fields of grounded glaciers under Glen's flow law. Using convex analysis arguments, we prove the existence and the uniqueness of a weak solution. A finite element method is applied with approximation spaces that satisfy the inf-sup condition, and a priori error estimates are established by using a quasinorm technique. Several algorithms (including Newton's method) are proposed to solve the nonlinearity of the Stokes problem and are proved to be convergent. Our results are supported by numerical convergence studies.

1. Introduction

In this paper we consider a model problem that is commonly used by glaciologists to compute the motion of glaciers. Ice is assumed to be an incompressible non-Newtonian fluid governed by Glen's law [1]. Glen's law and the mass momentum equation lead to a nonlinear stationary Stokes problem with a strain-dependent viscosity.

Glacier models based on Glen's law have already been studied by several authors. However, all of them have considered a simplified model, called first-order approximation [2]. This model is obtained by rewriting the Stokes equations into a dimensionless form and by dropping all terms of order 𝒪(𝜖2), where 𝜖 is the typical aspect ratio of glaciers. This simplification results into a nonlinear elliptic problem for the horizontal velocity field, the vertical component, and the pressure field being determined a posteriori. Colinge and Rappaz first demonstrated the well-posedness of this problem and proved the convergence of the finite element approximation with piecewise linear continuous functions in [3]. Inspired by the work of Baranger/Najib [4] and Barrett/Liu [5] on non-Newtonian problems, a priori and a posteriori error estimates were obtained later in [68].

Unlike the first-order approximation, the original Stokes model which is considered in this paper is a saddle point problem for the velocity and the pressure fields. We prove the existence and the uniqueness of a weak solution using an equivalent minimisation problem and an inf-sup stability condition. Next, we establish a priori estimates for a finite element approximation using a quasinorm technique [5]. Eventually, we investigate several successive approximation algorithms to solve the system nonlinearity. In particular, we upgrade by using Newton's method the fixed point algorithm, given in [3] and proved to be convergent in [6, 9].

Boundary conditions describe the basal sliding phenomena that can significantly influence the glacier ice flows. In [3, 68], the first order approximation model was coupled to a Dirichlet condition. However, this approach requires the basal velocity distribution which is unknown. To overcome this difficulty, several sliding laws—including a Coulomb-type law—were considered in [10]. In our model, we use a sliding law that results in a nonlinear Dirichlet-Robin boundary condition.

This paper is organised as follows: the physical model is presented in Section 2. We prove the well-posedness of the weak problem in Section 3. In Section 4, we apply a finite element method and establish a priori error estimates. Successive approximation algorithms to solve the system nonlinearity are proposed and proved to be convergent in Section 5. In Section 6, convergence studies are performed to support the results of Sections 4 and 5.

2. The Model

Let us suppose that ice occupies the domain Ω𝑑, with 𝑑=2 or 3. Ice can be considered as an incompressible non-Newtonian fluid with negligible inertial effects [11]. It follows that the velocity 𝐮 and the pressure 𝑝 of ice solve the stationary nonlinear Stokes problem in Ω: 2div(𝜇𝜀(𝐮))+𝑝=𝐟,div(𝐮)=0,(2.1) where 𝜀(𝐮)=(1/2)(𝐮+𝐮𝑇) denotes the rate of strain tensor, 𝜇 the viscosity of ice, and 𝐟 the gravity force. Here above, the viscosity 𝜇 depends on |𝜀(𝐮)|=𝜀(𝐮)𝜀(𝐮) and is defined by the regularised Glen's flow law [11]. More precisely, for a given velocity field 𝐮, the viscosity 𝜇 satisfies the following nonlinear equation: 12𝜇=𝐴𝜏𝑛10+2𝜇𝑠𝑛1,(2.2) where 𝑠=|𝜀(𝐮)|, 𝐴 is a positive parameter, 𝑛1 is Glen's exponent, and 𝜏0>0 is a small regularization parameter which prevents infinite viscosity for zero strain (𝜏0=0 in the original Glen's law [1]). When 𝑛=1, then the viscosity 𝜇 is constant and (2.1) correspond to the classical linear Stokes problem related to a Newtonian fluid. In the framework of glaciology, 𝑛 is often taken equal to 3; see [12].

Let us set the boundary conditions for the system of (2.1). Three mechanical circumstances may occur at the boundary of a glacier: (i) no force applies on the ice-air interface; (ii) ice slides on the bedrock-ice interface; (iii) ice is stuck to the bedrock-ice interface. The boundary of Ω is thus split into three parts: Γ𝑁, Γ𝑅, and Γ𝐷, referring to circumstances (i), (ii), and (iii), respectively. We assume throughout that Ω is bounded, its boundaries Γ𝑁 and Γ𝑅, are 𝒞1 and Γ𝐷. We consider the free surface condition: 2𝜇𝜀(𝐮)𝐧𝑝𝐧=𝟎,onΓ𝑁,(2.3) where 𝐧 is the unit outward normal vector along the boundary of the domain Ω. We apply the nonlinear sliding condition [10, 13, 14]: 𝐮𝐧=0,(2𝜇𝜀(𝐮)𝐧)𝐭𝑖=𝛼𝐮𝐭𝑖,𝑖=1,𝑑1onΓ𝑅,(2.4) where {𝐭𝑖}𝑖=1,𝑑1 are the orthogonal vectors tangent to the boundary Γ𝑅, that is, 𝐭1 when 𝑑=2 and 𝐭1,𝐭2 when 𝑑=3. Here above, 𝛼=𝛼(|𝐮|) is the sliding coefficient that is given by 𝛼(𝑡)=𝑐𝑡+𝑡01/𝑛1,(2.5) where 𝑡=|𝐮| is the Euclidean norm of 𝐮, 𝑛 is Glen's exponent, 𝑐 is a positive parameter, and 𝑡0>0 is a small parameter which prevents infinite 𝛼 for zero velocity. The no-sliding condition writes 𝐮=𝟎,onΓ𝐷.(2.6)

Note that the conditions applied on boundaries Γ𝑁, Γ𝑅, and Γ𝐷 are Neumann, Robin-Dirichlet, and Dirichlet conditions, respectively. When 𝑛=1 (Newtonian flow) and Γ𝑅=, the problem (2.1) with boundary conditions (2.3), (2.6) has already been widely studied; see, for instance, [1517].

3. Existence and Uniqueness

In this section, we prove that there exists a unique weak solution to problem (2.1) with mixed boundary conditions (2.3), (2.4), and (2.6). Pressure is first eliminated from the system by restricting the velocity space to divergence-free fields. Afterwards, the reduced problem is transformed into a minimisation problem. Following [3, 8], its well-posedness is proved by using convex analysis arguments. The existence and the uniqueness of the pressure field are ensured by an inf-sup condition. We now state in the next lemma several properties of the function 𝜇 that will often be used in Sections 3, 4 and 5.

Lemma 3.1. For all 𝑠+, there exists a unique 𝜇=𝜇(𝑠)+ satisfying (2.2). The function 𝑠𝜇(𝑠) is 𝒞(0,+) and decreasing. There exist 𝐷1,𝐷2,𝐷3,𝐷4>0 such that: 𝐷1(1+𝑠+𝑡)11/𝑛(𝑠𝑡)𝑠𝜇(𝑠)𝑡𝜇(𝑡)𝐷2(1+𝑠+𝑡)1(1/𝑛)(𝑠𝑡),𝑠𝑡0,(3.1)𝐷1(1+𝑠)11/𝑛𝜇(𝑠)𝐷2(1+𝑠)11/𝑛,𝑠0,(3.2)1𝑛𝜇(𝑠)𝜇(𝑠)+𝑠𝜇(𝑠),𝑠>0,(3.3)𝜇||𝜉||𝜉𝜇||𝜂||𝜂(𝜉𝜂)𝐷31+||𝜉||+||𝜉𝜂||(1/𝑛1)||𝜉𝜂||2,𝜉,𝜂𝑑×𝑑,(3.4)||𝜇||𝜉||𝜉𝜇||𝜂||𝜂||𝐷41+||𝜉||+||𝜉𝜂||(1/𝑛1)||𝜉𝜂||,𝜉,𝜂𝑑×𝑑.(3.5)

Proof. The properties of 𝜇 and inequalities (3.1) and (3.2) can be easily deduced from Lemmas  1 and 2 of [7]. Inequality (3.3) is obtained by differentiating (2.2) with respect to 𝑠. Inequalities (3.4) and (3.5) result from inequality (3.1), Lemma  2.1 in [5] and inequality (1/2)(|𝜉|+|𝜂|)|𝜉|+|𝜉𝜂|2(|𝜉|+|𝜂|). Details are given in [12].

Let us notice that property (3.1) was introduced by Barrett and Liu (see [5]) in order to obtain a priori error estimates of a similar problem to the one treated in this paper. Define the Banach spaces: 𝑉=𝐯𝑊1,𝑟(Ω)𝑑,𝐯=𝟎onΓ𝐷,𝐯𝐧=0onΓ𝑅,𝑄=𝐿𝑟(Ω),(3.6) where 𝑟=1+1𝑛,𝑟=𝑛+1(3.7) are conjugate exponents and 𝑛 is Glen's exponent. By using (3.2), we have 𝜇(𝑠)𝑠𝐶𝑠𝑟1 for all 𝑠>0. Then, if 𝐮𝑉, we have 𝜇(|𝜀(𝐮)|)𝜀(𝐮)[𝐿𝑟(Ω)]𝑑×𝑑. By using the trace inequality 𝐯𝐿𝑟(Γ𝑅)𝐯𝑊11/𝑟,𝑟(Γ𝑅)𝐶𝐯𝑊1,𝑟(Ω) for all 𝐯[𝑊1,𝑟(Ω)]𝑑, see [19 page 197], we obtain (𝐮𝐭𝑖)𝐿𝑟(Γ𝑅),𝑖=1,𝑑1. Similarly, we can show 𝛼(|𝐮|)(𝐮𝐭𝑖)𝐿𝑟(Γ𝑅),𝑖=1,𝑑1. Owing to Hölder's inequality, the mixed formulation of problem (2.1) with boundary conditions (2.3), (2.4), and (2.6) that consists of finding (𝐮,𝑝)𝑉×𝑄 such that 2Ω𝜇||𝜀(𝐮)||𝜀(𝐮)𝜀(𝐯)𝑑𝑉+𝑖=1,𝑑1Γ𝑅𝛼(|𝐮|)𝐮𝐭𝑖𝐯𝐭𝑖𝑑𝑆,Ω𝑝div(𝐯)𝑑𝑉+Ω𝑞div(𝐮)𝑑𝑉=Ω𝐠𝐯𝑑𝑉,(𝐯,𝑞)𝑉×𝑄(3.8) is meaningful.

Remark 3.2. If Γ𝑁=, pressure 𝑝 in (3.8) is defined up to a constant. In that case, 𝑄=𝐿𝑟 is replaced by 𝑄=𝐿𝑟0={𝑞𝐿𝑟,Ω𝑞𝑑𝑉=0}. Moreover, if 𝑛=1, then 𝑟=2 and 𝜇 is constant and if Γ𝑅=, then the (linear) problem (3.8) is well posed; see, for instance, [1517].

The next lemma states the equivalence of norms |𝜀()|𝐿𝑟 and 𝑊1,𝑟 on space 𝑉.

Lemma 3.3 (Korn's inequality). If Γ𝐷 and if 1<𝛾<, then there exists a constant 𝐶>0 such that 𝐯𝑊1,𝛾𝐶||𝜀(𝐯)||𝐿𝛾,(3.9) for all 𝐯𝑊1,𝛾(Ω) such that 𝐯=𝟎 on Γ𝐷.

Proof. We apply Corollary  4.1 in [18] (𝐹 being the identity matrix) and Lemma  3.1 page 40 in [16].

We consider the divergence-free velocity space: 𝑉div=𝐯𝑊1,𝑟(Ω)𝑑,div(𝐯)=0,𝐯=𝟎onΓ𝐷,𝐯𝐧=0onΓ𝑅.(3.10)

In 𝑉div, the pressure field 𝑝 vanishes of the variational formulation (3.8). The reduced formulation consists then of finding 𝐮𝑉div such that 2Ω𝜇||𝜀(𝐮)||𝜀(𝐮)𝜀(𝐯)𝑑𝑉+𝑖=1,𝑑1Γ𝑅𝛼(|𝐮|)𝐮𝐭𝑖𝐯𝐭𝑖𝑑𝑆=Ω𝐠𝐯𝑑𝑉,𝐯𝑉div.(3.11) To transform problem (3.11) into a minimisation problem, we introduce the functional 𝐽(𝐮)=Ω𝑀||𝜀(𝐮)||𝑑𝑉+12Γ𝑅𝑁(|𝐮|)𝑑𝑆Ω𝐮𝐟𝑑𝑉,(3.12) where 𝑀(𝑥)=𝑥0𝑠𝜇(𝑠)𝑑𝑠,𝑁(𝑥)=𝑥0𝑡𝛼(𝑡)𝑑𝑡.(3.13) The functional 𝐽 is Gâteaux differentiable, and its first derivative 𝐷𝐽, at point 𝐮𝑉div, in direction 𝐯𝑉div, is given by 𝐷𝐽(𝐮),𝐯=2Ω𝜇||𝜀(𝐮)||𝜀(𝐮)𝜀(𝐯)𝑑𝑉+𝑖=1,𝑑1Γ𝑅𝛼(|𝐮|)𝐮𝐭𝑖𝐯𝐭𝑖𝑑𝑆Ω𝐠𝐯𝑑𝑉.(3.14) Clearly, any minimiser of 𝐽 in 𝑉div satisfies (3.11). We now establish several lemmas that allow us to prove the existence and the uniqueness of this minimiser in Theorem 3.8. We show the continuity of 𝐽 in Lemma 3.5, the strict convexity of 𝐽 in Lemma 3.6, and the coercivity (in the sense of (3.18)) of 𝐽 in Lemma 3.7. The continuity of 𝐽 requires the following result (Lemma  4 in [3])

Lemma 3.4. Let 𝑂 be a measurable set of 𝑑 and 𝑓,𝑔𝐿𝑟(𝑂), then one has the following inequality: 𝑂||||𝑓||𝑟||𝑔||𝑟||𝑑𝑉𝑟||𝑓||+||𝑔||𝑟1𝐿𝑟(𝑂)𝑓𝑔𝐿𝑟(𝑂).(3.15)

Lemma 3.5. The functional 𝐽 is 𝑊1,𝑟-continuous.

Proof. By using (3.2), (2.5), and 11/𝑛=2𝑟, we have, for all 𝐮,𝐯𝑉div||𝑀||𝜀(𝐮)||𝑀||𝜀(𝐯)||||=|||||𝜀(𝐮)|||𝜀(𝐯)||𝑠𝜇(𝑠)𝑑𝑠||||||||𝐷1|𝜀(𝐮)|||𝜀(𝐯)||𝑠𝑟1𝑑𝑠||||=𝐷1||||||𝜀(𝐮)||𝑟||𝜀(𝐯)||𝑟𝑟||||,||𝑁(|𝐮|)𝑁(|𝐯|)||=|||||𝐮||𝐯|𝑡𝛼(𝑡)𝑑𝑡||||||||𝑐|𝐮||𝐯|𝑡𝑟1𝑑𝑡||||=𝑐|||||𝐮|𝑟|𝐯|𝑟𝑟||||.(3.16) These two inequalities together with Lemma 3.4 imply the 𝑊1,𝑟-continuity of 𝐽.

Lemma 3.6. The functional 𝐽 is strictly convex on 𝑉.

Proof. Clearly, 𝑀(𝑠)=𝑠𝜇(𝑠) and 𝑀(𝑠)=𝑠𝜇(𝑠)+𝜇(𝑠). From (3.3), we have 𝑀(𝑠)=𝑠𝜇(𝑠)+𝜇(𝑠)(1/𝑛)𝜇(𝑠)>0 if 𝑠>0, and then 𝑀 is strictly convex. Since 𝑀 is an increasing function, 𝑀(||) is strictly convex. In the same way, we can show that 𝑁(||) is strictly convex by using (2.5). Let 𝐮,𝐯𝑉div satisfying 𝐮𝐯 and 𝜃(0,1). From Korn's inequality (Lemma 3.3), we have 𝜀(𝐮)𝜀(𝐯) in 𝐿𝑟. As a consequence, Ω𝑀||𝜃𝜀(𝐮)+(1𝜃)𝜀(𝐯)||𝑑𝑉<𝜃Ω𝑀||𝜀(𝐮)||𝑑𝑉+(1𝜃)Ω𝑀||𝜀(𝐯)||𝑑𝑉.(3.17) The strict convexity of 𝐽 follows from the previous inequality and the convexity of 𝑁(||).

Since 𝐽 is convex, 𝐮𝑉div satisfies (3.11) if and only if 𝐽(𝐮)𝐽(𝐯),𝐯𝑉div.

Lemma 3.7. There exist two constants 𝐷1,𝐷2>0 such that, for all 𝐮𝑉, 𝐽(𝐮)𝐷1𝐮𝑟𝑊1,𝑟𝐷2.(3.18)

Proof. Let 𝐮𝑉. From (3.2) and 11/𝑛=2𝑟, there exists 𝐶0>0 such that 𝑀||𝜀(𝐮)|||𝜀(𝐮)|0𝐶0𝑠(1+𝑠)2𝑟𝑑𝑠=|𝜀(𝐮)|2/20𝐶01+2𝑡2𝑟𝑑𝑡||𝜀(𝐮)||22𝐶01+||𝜀(𝐮)||𝑟2.(3.19) As a consequence, there exist two constants 𝐶1,𝐶2>0 such that 𝑀(|𝜀(𝐮)|)𝐶1(1+|𝜀(𝐮)|)𝑟𝐶2. By using Korn's inequality (Lemma 3.3), there exists 𝐶3>0 such that Ω𝑀||𝜀(𝐮)||𝑑𝑉𝐶3𝐮𝑟𝑊1,𝑟𝐶2Ω𝑑𝑉.(3.20) From Young's inequality, we have, for all 𝛿>0, Ω||𝐮𝐟||𝑑𝑉Ω1𝑟𝛿𝑟||𝐟||𝑟+𝛿𝑟𝑟|𝐮|𝑟𝑑𝑉=𝐶4𝛿𝑟Ω𝑑𝑉+𝐶5𝛿𝑟𝐮𝑟𝐿𝑟,(3.21) where 𝐶4,𝐶5>0. We set 𝛿 small enough such that 𝐶3𝛿𝑟𝐶5>0. From inequalities (3.20), (3.21), and 𝑁0, we obtain 𝐽(𝐮)=Ω𝑀||𝜀(𝐮)||𝐮𝐟𝑑𝑉+Γ𝑅12𝑁(|𝐮|)𝑑𝑆𝐶3𝐮𝑟𝑊1,𝑟𝐶2Ω𝑑𝑉𝐶4𝛿𝑟Ω𝑑𝑉𝐶5𝛿𝑟𝐮𝑟𝑊1,𝑟,(3.22) which is exactly (3.18) with 𝐷1=𝐶3𝐶5𝛿𝑟 and 𝐷2=(𝐶2+𝐶4/𝛿𝑟)Ω𝑑𝑉.

Theorem 3.8. There exists a unique 𝐮𝑉div such that 𝐽(𝐮)=inf{𝐽(𝐯);𝐯𝑉div}. Moreover, 𝐮 is the unique solution of (3.11).

Proof. Clearly, there exists 𝐮𝑉div such that 𝐽(𝐮)<+. Lemma 3.7 ensures the existence of 𝑚=inf{𝐽(𝐯);𝐯𝑉div}. Let {𝐮𝜈} be a sequence of 𝑉div such that lim𝜈𝐽(𝐮𝜈)=𝑚. There exists an integer 𝐾 such that, for all 𝜈>𝐾, we have 𝑚+1>𝐽(𝐮𝜈). Owing to Lemma 3.7, the sequence {𝐮𝜈} is bounded in 𝑉div. Since 𝑉div is a closed subspace of 𝑉, 𝑉div is reflexive. Consequently, there exist 𝐮𝑉div and a subsequence of {𝐮𝜈} (still denoted {𝐮𝜈}) that converges weakly to 𝐮 in 𝑉div. By Lemmas 3.5 and 3.6, 𝐽 is weakly lower semicontinuous; see, for instance, Corollary III.8 in [19] page 38. Then, we have 𝑚=liminf𝜈+𝐽𝐮𝜈𝐽(𝐮)𝑚,(3.23) and 𝐽 possesses at least one minimum 𝐮𝑉div. Since 𝐽 is strictly convex (Lemma 3.6), this minimum is unique. Moreover, 𝐮 is the unique solution of (3.11).

Spaces 𝑉 and 𝑄 are required to satisfy the inf-sup condition, see [5, 20], to ensure the existence and the uniqueness of 𝑝𝑄 such that (𝐮,𝑝) satisfies the mixed formulation (3.8). The inf-sup condition is proved in [15, 21] when Γ𝐷=𝜕Ω (or, equivalently, Γ𝑁=Γ𝑅=). By following the proof of Proposition  5.3.2 in [22], we can easily generalise this result when Γ𝑅Γ𝑁; see details in [12].

Lemma 3.9. Spaces 𝑉 and 𝑄 satisfy the inf-sup condition; that is, there exists 𝐶>0 such that 𝐶<inf𝑞𝑄sup𝐯𝑉Ω𝑞div(𝐯)𝑑𝑉𝑞𝐿𝑟𝐯𝑊1,𝑟.(3.24)

Theorem 3.10. There exists a unique couple (𝐮,𝑝)(𝑉,𝑄) satisfying (3.8).

Proof. Although the result is a straightforward application of Theorem  2.1 in [20] together with Theorem 3.8 and Lemma 3.9, we give all the arguments of the proof. Let 𝐴𝑉𝑉 and 𝐵𝑉𝑄 be the operators defined by 𝐴𝐮,𝐯=2Ω𝜇||𝜀(𝐮)||𝜀(𝐮)𝜀(𝐯)𝑑𝑉+𝑖=1,𝑑1Γ𝑅𝛼(|𝐮|)𝐮𝐭𝑖𝐯𝐭𝑖𝑑𝑆,𝑣𝑉,𝐵𝐮,𝑞=Ω𝑞div𝐮𝑑𝑉,𝑞𝑄,(3.25) where 𝑉 and 𝑄 are dual to 𝑉 and 𝑄, respectively. From Theorem 3.8, there exists a unique 𝐮ker𝐵 such that 𝐴𝐮𝐟,𝐯=0 for all 𝐯ker𝐵, which means that 𝐴𝐮𝐟(ker𝐵). Owing to the inf-sup condition (4.1), the operator 𝐵𝑉𝑄 is surjective, ker𝐵𝑇=, and (𝐵𝑇) is closed; see Lemma  A.40 in [15]. As a consequence, 𝐴𝐮𝐟(ker𝐵)=(𝐵𝑇)=(𝐵𝑇) and there exists 𝑝𝑄 such that 𝐴𝐮𝐟=𝐵𝑇𝑝. Since ker𝐵𝑇=, the pressure 𝑝 is necessarily unique. Eventually, there exists a unique couple (𝐮,𝑝)𝑉×𝑄 satisfying 𝐴𝐮𝐵𝑇𝑝=𝐟,𝐵𝐮=0,(3.26) or equivalently (3.8).

4. Finite Element Approximation and A Priori Estimates

We assume that Ω is a convex polygonal or polyhedral domain and 𝒯 is a regular mesh of Ω parametrized by , the highest diameter of the elements of 𝒯. We say that 𝑉𝑉 and 𝑄𝑄, some finite-dimensional approximation spaces on 𝒯 of 𝑉 and 𝑄, satisfy the inf-sup condition if, for all 𝜅(1,), there exists a constant 𝐶>0 such that 𝐶<inf𝑞𝑄sup𝐯𝑉Ω𝑞div𝐯𝑑𝑉𝑞𝐿𝜅𝐯𝑊1,𝜅.(4.1) The discrete problem is obtained by replacing the spaces 𝑉 and 𝑄 by 𝑉 and 𝑄, respectively. It consists of finding (𝐮,𝑝)(𝑉,𝑄) such that 2Ω𝜇||𝜀𝐮||𝜀𝐮𝜀𝐯𝑑𝑉+𝑖=1,𝑑1Γ𝑅𝛼||𝐮||𝐮𝐭𝑖𝐯𝐭𝑖𝑑𝑆,Ω𝑝div𝐯𝑑𝑉+Ω𝑞div𝐮𝑑𝑉=Ω𝐟𝐯𝑑𝑉,𝐯,𝑞𝑉,𝑄.(4.2) The discrete similar space to 𝑉div is 𝑉div,=𝐯𝑉;Ω𝑞div𝐯𝑑𝑉=0;𝑞𝑄.(4.3) Note that 𝑉div, is not necessarily included in 𝑉div. The discrete reduced problem consists of finding 𝐮𝑉div, such that 2Ω𝜇||𝜀𝐮||𝜀𝐮𝜀𝐯𝑑𝑉+𝑖=1,𝑑1Γ𝑅𝛼||𝐮||𝐮𝐭𝑖𝐯𝐭𝑖𝑑𝑆=Ω𝐠𝐯𝑑𝑉,𝐯𝑉div,.(4.4) Since 𝑉div, is a closed subspace of 𝑉, Theorem 3.8 and the proof can be rewritten by replacing 𝑉div by 𝑉div,.

Theorem 4.1. There exists a unique 𝐮𝑉div, such that 𝐽(𝐮)=inf{𝐽(𝐯);𝐯𝑉div,}. Moreover, 𝐮 is the unique solution of (4.4).

Remark 4.2. By setting 𝐯=𝐮 in (4.4) and by using inequality (3.2), (2.5), and Korn's inequality (Lemma 3.3), we can show that the solution 𝐮 of problem (4.4) satisfies 𝐮𝑊1,𝑟𝐶𝐟𝐿𝑟,(4.5) where 𝐶>0 does not depend on 𝐮.

From Theorem 4.1 and the inf-sup condition (4.1), we can rewrite Theorem 3.10 and its proof for the discrete mixed problem.

Theorem 4.3. If 𝑉 and 𝑄 satisfy the inf-sup condition (4.1), then there exists a unique couple (𝐮,𝑝)(𝑉,𝑄) satisfying (4.2).

Remark 4.4. The spaces [1/Bulle]𝑑1 and [2]𝑑1 are two examples that satisfy the inf-sup condition (4.1) while 11 does not satisfy (4.1); see [15].

The error analysis that follows is partly inspired from [5, 7]. We give a priori estimates for the numerical approximation of the stationary Stokes problem in Theorem 4.9. For the sake of simplicity, we suppose Γ𝑅=; that is, the boundary Robin-Dirichlet condition is not considered; see also Remark 4.11. The nonlinearity of problem (3.8) is treated by introducing (in Lemma 4.5) a quasi-norm that depends on the solution; see [5]. The orthogonality of the error (Lemma 4.6) together with properties (3.4) and (3.5) of the function 𝜇 allow quasi-norm estimates to be established in Theorem 4.7. The properties of the quasi-norm given in Lemma 4.5 allow estimates with standard norms to be proved in Theorem 4.8. Eventually, these estimates together with interpolation inequalities yield to the main Theorem 4.9.

Lemma 4.5. Let (𝐮,𝑝) be the solution of (3.8); the application 𝐯|𝐯|=Ω||𝜀(𝐯)||21+||𝜀(𝐮)||+||𝜀(𝐯)||2𝑟𝑑𝑉(4.6) is a quasi-norm of 𝑉; that is, it satisfies all properties of norms, except homogeneity. Moreover, there exists 𝐷1>0 such that, for all 𝐯𝑊1,𝑟(Ω), one has 𝐯2𝑊1,𝑟𝐷11+𝐮𝑊1,𝑟+𝐯𝑊1,𝑟2𝑟|𝐯|2,(4.7) and there exists 𝐷2>0 such that, for all 𝜅[𝑟,2] and for all 𝐯𝑊1,𝜅(Ω), one has |𝐯|2𝐷2𝐯𝜅𝑊1,𝜅.(4.8)

Proof. The quasi-norm properties are shown in Lemma  3.1 in [5]. Inequalities (4.7) and (4.8) result from Korn and Hölder's inequalities; see details in [12].

By setting 𝐯=𝐯𝑉div, in (3.8) it is easy to prove the next lemma.

Lemma 4.6. Let (𝐮,𝑝)(𝑉,𝑄) be the solution of problem (3.8) and 𝐮𝑉div, the solution of problem (4.4), then Ω2𝜇||𝜀(𝐮)||𝜀(𝐮)𝜇||𝜀𝐮||𝜀𝐮𝜀𝐯𝑑𝑉Ωdiv𝐯𝑝𝑞𝑑𝑉=0(4.9) holds for all (𝐯,𝑞)(𝑉div,,𝑄). Moreover, if the spaces 𝑉 and 𝑄 satisfy the inf-sup condition (4.1), then the solution (𝐮,𝑝) of (4.2) satisfies Ω2𝜇||𝜀(𝐮)||𝜀(𝐮)𝜇||𝜀𝐮||𝜀𝐮𝜀𝐯𝑑𝑉Ωdiv𝐯𝑝𝑝𝑑𝑉=0,(4.10) for all 𝐯𝑉.

Theorem 4.7. Let (𝐮,𝑝)(𝑉,𝑄) be the solution of (3.8) and 𝐮𝑉div, the solution of (4.4). For all (𝐯,𝑞)𝑉div,×𝑄, one has ||𝐮𝐮||𝐷11+𝐮𝑊1,𝑟+𝐯𝑊1,𝑟(2𝑟)/2||𝐮𝐯||+𝑝𝑞𝐿𝑟,(4.11) where 𝐷1>0. Moreover, if the spaces 𝑉 and 𝑄 satisfy the inf-sup condition (4.1), then the solution (𝐮,𝑝) of (4.2) satisfies, for all 𝜅[𝑟,2] and for all 𝑞𝑄, 𝑝𝑝𝐿𝜅𝐷2||𝐮𝐮||2/𝜅+𝑝𝑞𝐿𝜅,(4.12) where 𝐷2>0. The constants 𝐷1,𝐷2 do not depend on 𝐮 and 𝐯; however, 𝐷2 increasingly depends on (𝐶)1.

Proof. By using, respectively, the definition (4.6) of the quasi-norm ||||||, inequality (3.4) with 11/𝑛=2𝑟, and (4.9), there exists 𝐶1>0 such that ||𝐮𝐮||2𝐶1Ω2𝜇||𝜀(𝐮)||𝜀(𝐮)𝜇||𝜀𝐮||𝜀𝐮𝜀𝐮𝐮𝑑𝑉=𝐶1Ω2𝜇||𝜀(𝐮)||𝜀(𝐮)𝜇||𝜀𝐮||𝜀𝐮𝜀𝐮𝐯𝑑𝑉=𝐴1+Ωdiv𝐯𝐮𝑝𝑞𝑑𝑉=𝐴2,(4.13) where (𝐯,𝑞)𝑉div,×𝑄. For the sake of simplicity, 𝐴1 and 𝐴2 are handled separately. By using inequality (3.5) with 11/𝑛=2𝑟, there exists 𝐶2>0 such that ||𝐴1||Ω2|||𝜇||𝜀(𝐮)||𝜀(𝐮)𝜇||𝜀𝐮||𝜀𝐮|||||𝜀𝐮𝐯||𝑑𝑉𝐶2Ω||𝜀𝐮𝐮||||𝜀𝐮𝐯||1+||𝜀(𝐮)||+||𝜀𝐮𝐮||2𝑟𝑑𝑉.(4.14) By using the inequality (see Lemma  2.2 in [5] or  (3.10) in [8]), (1+𝑎+𝑡)𝑒𝑡𝑠𝛼(1+𝑎+𝑡)𝑒𝑡2+𝛼1(1+𝑎+𝑠)𝑒𝑠2,𝑎,𝑡,𝑠0,𝛼(0,1],𝑒(0,1),(4.15) with 𝑎=|𝜀(𝐮)|, 𝑡=|𝜀(𝐮𝐮)|, 𝑠=|𝜀(𝐮𝐯)|, and 𝑒=2𝑟, we obtain, for all 𝛼[0,1], ||𝐴1||𝐶2𝛼||𝐮𝐮||2+𝛼1||𝐮𝐯||2.(4.16) We now use respectively Hölder's inequality, Young's inequality, and (4.7); there exist 𝐶3,𝐶4,𝐶5>0 such that ||𝐴2||div𝐯𝐮𝐿𝑟𝑝𝑞𝐿𝑟12𝐶3𝛽𝐯𝐮2𝑊1,𝑟+12𝛽1𝑝𝑞2𝐿𝑟𝐶4𝛽1+𝐮𝑊1,𝑟+𝐯𝐮𝑊1,𝑟2𝑟||𝐯𝐮||2+12𝛽1𝑝𝑞2𝐿𝑟𝐶5𝛽1+𝐮𝑊1,𝑟+𝐯𝑊1,𝑟2𝑟||𝐮𝐮||2+||𝐮𝐯||2+12𝛽1𝑝𝑞2𝐿𝑟.(4.17) By setting 𝛼=1/(4𝐶2𝐶1) and 𝛽=1/(4𝐶5𝐶1[1+𝐮𝑊1,𝑟+𝐯𝑊1,𝑟]2𝑟), we obtain ||𝐮𝐮||2𝐶1||𝐴1||+||𝐴2||12||𝐮𝐮||2+4𝐶22𝐶21||𝐮𝐯||2+14||𝐮𝐯||2+2𝐶21𝐶51+𝐮𝑊1,𝑟+𝐯𝑊1,𝑟2𝑟𝑝𝑞2𝐿𝑟(Ω).(4.18) By moving (1/2)|||𝐮𝐮|||2 to the left-hand side, we obtain (4.11). From the inf-sup condition (4.1), we have, for all 𝑞𝑄, 𝐶𝑞𝑝𝐿𝜅<sup𝐯𝑉Ω𝑞𝑝div𝐯𝑑𝑉𝐯𝑊1,𝜅.(4.19) From (4.10), we have, for all 𝐯𝑉, Ωdiv𝐯𝑞𝑝𝑑𝑉=Ω2𝜇||𝜀(𝐮)||𝜀(𝐮)𝜇||𝜀𝐮||𝜀𝐮𝜀𝐯𝑑𝑉+Ωdiv𝐯𝑞𝑝𝑑𝑉.(4.20) From (4.19), (4.20), and (3.5) with 11/𝑛=2𝑟, there exist 𝐶6,𝐶7>0 such that 𝐶𝑞𝑝𝐿𝜅<sup𝐯𝑉|||Ω2𝜇||𝜀(𝐮)||𝜀(𝐮)𝜇||𝜀𝐮||𝜀𝐮𝜀𝐯𝑑𝑉|||𝐯𝑊1,𝜅+sup𝐯𝑉||Ωdiv𝐯𝑞𝑝𝑑𝑉||𝐯𝑊1,𝜅𝐶62𝜇||𝜀(𝐮)||𝜀(𝐮)𝜇||𝜀𝐮||𝜀𝐮𝐿𝜅+𝑞𝑝𝐿𝜅𝐶7Ω1+||𝜀(𝐮)||+||𝜀𝐮𝐮||(𝑟2)𝜅||𝜀𝐮𝐮||𝜅𝑑𝑉1/𝜅+𝐶6𝑞𝑝𝐿𝜅𝐶7𝐶8Ω1+||𝜀(𝐮)||+||𝜀𝐮𝐮||(𝑟2)||𝜀𝐮𝐮||2𝑑𝑉1/𝜅+𝐶6𝑞𝑝𝐿𝜅𝐶7||𝐮𝐮||2/𝜅+𝐶6𝑞𝑝𝐿𝜅,(4.21) where 𝐶8=(1+|𝜀(𝐮)|+|𝜀(𝐮𝐮)|)(𝑟2)|𝜀(𝐮𝐮)|2𝜅1/𝜅𝐿<1. Eventually, the previous inequality together with 𝑝𝑝𝐿𝜅𝑞𝑝𝐿𝜅+𝑝𝑞𝐿𝜅 leads to (4.12).

Theorem 4.8. Let (𝐮,𝑝)(𝑉,𝑄) be the solution of (3.8), and 𝐮𝑉div, the solution of (4.4). For all (𝐯,𝑞)𝑉div,×𝑄 and for all 𝜅[𝑟,2], assuming 𝐮𝑊1,𝜅(Ω), one has 𝐮𝐮𝑊1,𝑟𝐷11+𝐮𝑊1,𝑟+𝐯𝑊1,𝑟2𝑟𝐮𝐯𝜅/2𝑊1,𝜅+𝑝𝑞𝐿𝑟,(4.22) where 𝐷1>0. Moreover, if the spaces 𝑉 and 𝑄 satisfy the inf-sup condition (4.1), then the solution (𝐮,𝑝) of (4.2) satisfies for all (𝐯,𝑞)𝑉×𝑄 and for all 𝜅[𝑟,2], assuming 𝐮𝑊1,𝜅(Ω), 𝐮𝐮𝑊1,𝑟𝐷21+𝐮𝑊1,𝑟+𝐯𝑊1,𝑟2𝑟𝐮𝐯𝜅/2𝑊1,𝜅+𝑝𝑞𝐿𝑟,(4.23)𝑝𝑝𝐿𝜅𝐷31+𝐮𝑊1,𝑟+𝐯𝑊1,𝑟(2𝑟)/𝜅×𝐮𝐯𝜅/2𝑊1,𝜅+𝑝𝑞𝐿𝑟2/𝜅+𝑝𝑞𝐿𝜅,(4.24) where 𝐷2,𝐷3>0. The constants 𝐷1,𝐷2,𝐷3 do not depend on 𝐮 and 𝐯; however, 𝐷2 and 𝐷3 increasingly depends on (𝐶)1.

Proof. On one hand, inequality (4.22) follows from inequalities (4.7), (4.11), and (4.8). On the other hand, (4.23) follows from (4.22) and from the following property (see (1.16), page 115 in [16]): for all 𝐯𝑉div and for all 𝐰𝑉, there exists 𝐯𝑉div, such that 𝐯𝐯𝑊1,𝜅𝐶𝐯𝐰𝑊1,𝜅,(4.25) where 𝐶 depends on the inf-sup constant 𝐶. Eventually, (4.24) follows from (4.12), (4.11), and (4.8).

Theorem 4.9. Assume that, for all 𝜅[𝑟,2], there exists a continuous operator 𝜋[𝑊2,𝜅]𝑑𝑉 that satisfies 𝐮𝜋(𝐮)𝑊1,𝜅𝐶𝐮𝑊2,𝜅,𝐮𝑊2,𝜅𝑑,(4.26) and a continuous operator 𝜌𝑊1,𝜅𝑄 that satisfies 𝑝𝜌(𝑝)𝐿𝜅𝐶𝑝𝑊1,𝜅,𝑝𝑊1,𝜅,(4.27) where is the size of the higher diameter of the elements of 𝒯. Assume that 𝑉 and 𝑄 satisfy the inf-sup condition (4.1). Let (𝐮,𝑝) be the solution of problem (3.8) and let (𝐮,𝑝) be the solution of problem (4.2). Assume that (𝐮,𝑝)([𝑊2,𝜅]𝑑,𝑊1,𝜅), where 𝜅[𝑟,2], then one has 𝐮𝐮𝑊1,𝑟+𝑝𝑝𝐿𝜅𝜅/2𝐷𝜅/2,(4.28) where 𝐷=𝐷(𝐮𝑊2,𝜅,𝑝𝑊1,𝑟,(𝐶)1)>0.

Proof. Apply (4.23) and (4.24) with 𝐯=𝜋(𝐮) and 𝑞=𝜌(𝑝). By using the continuity of 𝜋, (4.5), (4.26), and (4.27), there exist 𝐶1,𝐶2,𝐶3,𝐶4>0 such that 𝐮𝐮𝑊1,𝑟𝐶11+𝐟𝐿𝑟+𝐮𝑊1,𝑟2𝑟(𝜅/2)+𝐶2(𝜅/2),𝑝𝑝𝐿𝜅𝐶31+𝐟𝐿𝑟+𝐮𝑊1,𝑟(2𝑟)/𝜅(𝜅/2)+(2/𝜅)+𝐶4(𝜅/𝜅).(4.29) The estimate (4.28) directly follows from (4.29).

Remark 4.10. The combination [1/Bulle]𝑑1 for spaces 𝑉 and 𝑄, introduced in [23], satisfies the assumptions of Theorem 4.9; see Lemma  4.20 page 190 of [15] for the inf-sup condition (4.1) and [15, 16] for the interpolation inequalities (4.26) and (4.27).

Remark 4.11. If Γ𝑅, a similar analysis can be led by replacing the norm defined by (4.6) by𝐯|𝐯|=Ω||𝜀(𝐯)||21+||𝜀(𝐮)||+||𝜀(𝐯)||2𝑟𝑑𝑉+Γ𝑅|𝐯|2(1+|𝐮|+|𝐯|)2𝑟𝑑𝑆.(4.30)

5. Successive Approximations

In this section, several successive approximation algorithms are proposed for solving the nonlinearity of the discrete problem (4.4) when 𝑛>1. For the sake of simplicity, we suppose Γ𝑅= in this section; see Remark 5.8. We present a unified scheme that contains the classical fixed point method together with Newton's method. The mesh 𝒯 is fixed, and we assume that the approximation spaces satisfy 𝑉𝑉[𝑊1,(Ω)]𝑑 and 𝑄𝑄𝐿(Ω). In what follows, denotes an arbitrary norm of 𝑉div,. Since 𝑉div, is a finite-dimensional space, all norms are equivalent. Let 𝛾[0,1]. We define 𝐸𝑉div,𝑉div,,̃𝐮𝐰,(5.1) where 𝐰𝑉div, solves 2Ω𝜇||𝜀̃𝐮||𝜀𝐰𝜀𝐯𝑑𝑉+2𝛾Ω𝜇||𝜀̃𝐮||||𝜀̃𝐮||𝜀̃𝐮𝜀𝐰̃𝐮𝜀̃𝐮𝜀𝐯𝑑𝑉=Ω𝐠𝐯𝑑𝑉,𝐯𝑉div,.(5.2)

The application 𝐸 is well defined. Indeed, by using, respectively, 𝜇<0, inequalities (3.3) and (3.2), there exist 𝐶1,𝐶2>0 such that 2Ω𝜇||𝜀̃𝐮||||𝜀𝐰||2𝑑𝑉+2𝛾Ω𝜇||𝜀̃𝐮||||𝜀̃𝐮||𝜀̃𝐮𝜀𝐰2𝑑𝑉2Ω𝜇||𝜀̃𝐮||||𝜀𝐰||2𝑑𝑉2𝛾(11/𝑛)Ω𝜇||𝜀̃𝐮||||𝜀𝐰||2𝑑𝑉2𝐶1(1𝛾(11/𝑛))1+𝜀̃𝐮𝐿11/𝑛Ω||𝜀𝐰||2𝑑𝑉𝐶2𝐰2.(5.3) As a consequence, the problem (5.2) is coercive. From the Lax-Milgram Theorem, see [15] page 83, there exists a unique solution 𝐰𝑉div, of (5.2).

In what follows, 𝐮 denotes the solution (4.4), which is also the unique fixed point of 𝐸. Assume that 𝐮,0 is given; we define iteratively a sequence 𝐮,𝑘, for all 𝑘1, by 𝐮,𝑘+1=𝐸𝐮,𝑘.(5.4) Our goal is to prove that 𝐮,𝑘 converges to 𝐮 when 𝑘 goes to the infinity. When 𝛾=0, we obtain the classical fixed point method, widely used to solve the nonlinearity of Glen's law; see [6, 9, 14]. When 𝛾=1, we have an additional term in (5.2) which corresponds to Newton's method; see Remark 5.5. The case 𝛾(0,1) corresponds to a hybrid fixed point—Newton's method. The convergence of sequence 𝐮,𝑘 requires several preliminary results. We compute the first derivative of 𝐸 in Lemma 5.1. Lemma 5.2 provides an upper bound of the first derivative. Eventually, Theorem 5.3 states the linear convergence of 𝐮,𝑘 by using the Banach fixed point theorem. Theorem 5.7 states the second-order convergence when 𝛾=1. By differentiating formally (5.2) at point ̃𝐮 in direction 𝐮, with 𝐰=𝐸(̃𝐮), we obtain the following lemma.

Lemma 5.1. Let (̃𝐮,𝐰) satisfy 𝐸(̃𝐮)=𝐰. The application 𝐸 is Gâteaux differentiable at point ̃𝐮 and its derivative is given by 𝐷𝐸̃𝐮𝑉div,𝑉div,,𝐮𝐰,(5.5) where 𝐰 solves 2Ω𝜇||𝜀̃𝐮||𝜀𝐰𝜀𝐯𝑑𝑉+2Ω𝜇||𝜀̃𝐮||||𝜀̃𝐮||𝜀̃𝐮𝜀𝐮𝜀𝐰𝜀𝐯𝑑𝑉+𝛾2Ω𝜇||𝜀̃𝐮||||𝜀̃𝐮||𝜇||𝜀̃𝐮||||𝜀̃𝐮||3𝜀̃𝐮𝜀𝐮×𝜀̃𝐮𝜀𝐰̃𝐮𝜀̃𝐮𝜀𝐯𝑑𝑉+𝛾2Ω𝜇||𝜀̃𝐮||||𝜀̃𝐮||𝜀𝐮𝜀𝐰̃𝐮𝜀̃𝐮𝜀𝐯𝑑𝑉+𝛾2Ω𝜇||𝜀̃𝐮||||𝜀̃𝐮||𝜀̃𝐮𝜀𝐰𝐮𝜀̃𝐮𝜀𝐯𝑑𝑉+𝛾2Ω𝜇||𝜀̃𝐮||||𝜀̃𝐮||𝜀̃𝐮𝜀𝐰̃𝐮𝜀𝐮𝜀𝐯𝑑𝑉=0,𝐯𝑉div,.(5.6)

The problems (5.2) and (5.6) have the same coercivity properties to compute 𝐰 (resp., 𝐰) from 𝐮 (resp., ̃𝐮). As a consequence, the problem (5.6) is well-posed by the Lax-Milgram theorem. To prove the convergence of the sequence 𝐮,𝑘, we look for a norm that makes 𝐸 a contraction at point 𝐮.

Lemma 5.2. Let 𝛾[0,1], and let 𝐮 be the fixed point of 𝐸. The application 𝐷𝐸(𝐮) satisfies ||𝐷𝐸𝐮||𝜇(1𝛾)(11/𝑛)1(11/𝑛)𝛾<1,(5.7) where ||||||𝜇 is the subordinated norm to 𝜇=Ω𝜇(|𝜀(𝐮)|)|𝜀()|2𝑑𝑉.

Proof. Since 𝐸(𝐮)=𝐮, then (5.6), with ̃𝐮=𝐮 and 𝐰𝐡=𝐮, is rewriten as Ω𝜇||𝜀𝐮||𝜀𝐰𝜀𝐯𝑑𝑉+(1𝛾)Ω𝜇||𝜀𝐮||||𝜀𝐮||𝜀𝐮𝜀𝐮𝜀𝐮𝜀𝐯𝑑𝑉+𝛾Ω𝜇||𝜀𝐮||||𝜀𝐮||𝜀𝐮𝜀𝐰𝜀𝐮𝜀𝐯𝑑𝑉=0,(5.8) for all 𝐯𝑉div,. From (5.8), 𝜇<0, and (3.3), we have Ω𝜇||𝜀𝐮||𝜀𝐰𝜀𝐯𝑑𝑉(1𝛾)Ω𝜇||𝜀𝐮||||𝜀𝐮||||𝜀𝐮||||𝜀𝐯||𝑑𝑉𝛾Ω𝜇||𝜀𝐮||||𝜀𝐮||||𝜀𝐰||||𝜀𝐯||𝑑𝑉(1𝛾)11𝑛Ω𝜇||𝜀𝐮||||𝜀𝐮||||𝜀𝐯||𝑑𝑉+𝛾11𝑛Ω𝜇||𝜀𝐮||||𝜀𝐰||||𝜀𝐯||𝑑𝑉.(5.9) By setting 𝐯=𝐰 in (5.9), we obtain 𝐰2𝜇(1𝛾)11𝑛Ω𝜇||𝜀𝐮||||𝜀𝐮||||𝜀𝐰||𝑑𝑉+𝛾11𝑛𝐰2𝜇.(5.10) From (5.10) and Cauchy-Schwarz's inequality, we obtain 𝐰2𝜇(1𝛾)(11/𝑛)1(11/𝑛)𝛾𝐮𝜇𝐰𝜇.(5.11) Eventually, (5.7) follows from the definition of norm ||||||𝜇.

Theorem 5.3. Let 𝛾[0,1], let 𝒯 be a given mesh of Ω, and let be a norm of 𝑉div,. There exist 𝛿>0 and 𝐶>0 such that if 𝐮,0𝐮<𝛿, then one has 𝐮,𝑘𝐮𝐶(1𝛾)(11/𝑛)1(11/𝑛)𝛾𝑘𝐮,0𝐮,(5.12) for all 𝑘0, and 𝐮,𝑘 is linearly convergent to 𝐮.

Proof. From Lemma 5.2, the spectral radius of 𝐷𝐸(𝐮) is lower than constant: (1𝛾)(11/𝑛)1(11/𝑛)𝛾,(5.13) which is lower than 1. The theorem is then a direct application of the Banach fixed point theorem in 𝑉div,.

Remark 5.4. In Theorem 5.3, it should be stressed that 𝛿 depends on . When 0 (i.e., if we replace 𝐮 by 𝐮), we cannot ensure Theorem 5.3 to remain true. Nevertheless, in practise, 𝛿 seems to be independent of , see Section 6.

When 𝛾=1, we have 𝐷𝐸(𝐮)=𝟎 from Lemma 5.2. It suggests that the convergence of sequence 𝐮,𝑘 is quadratic. To establish the second-order convergence, we define, for all 𝐯𝑉div,, the application (,𝐯)𝑉div, by ̃𝐮;𝐯=2Ω𝜇||𝜀̃𝐮||𝜀̃𝐮𝜀𝐯𝑑𝑉Ω𝐟𝐯𝑑𝑉.(5.14)

Let ̃𝐮,𝐮,𝐮𝑉div,. We compute formally the first-order derivative of at point ̃𝐮 in direction 𝐮: 𝐷̃𝐮;𝐯,𝐮=2Ω𝜇||𝜀̃𝐮||𝜀𝐮𝜀𝐯𝑑𝑉+2Ω𝜇