Pruebas

Spectral discretizations of the Darcy’s equations with non standard boundary conditions

J. M. Bernard

*IUT d’Evry Val d’Essonne, 22,
allee Jean-Rostand, 91025 Evry Cedex, France., Francia

ACI Avances en Ciencias e Ingenierías

Universidad San Francisco de Quito, Ecuador

Received: 27 February 2017

Accepted: 09 July 2021

Abstract: This paper is devoted to the approximation of a nonstandard Darcy problem, which modelizes the ow in porous media, by spectral methods: the pressure is assigned on a part of the boundary. We propose two variational formulations, as well as three spectral discretizations. The second discretization improves the approximation of the divergence-free condition, but the error estimate on the pressure is not optimal, while the third one leads to optimal error estimate with a divergence-free discrete solution, which is important for some applications. Next, their numerical analysis is performed in detail and we present some numerical experiments which conrm the interest of the third discretization.

Keywords: Spectral methods, Darcy equations, boundary conditions, numerical results.

Introduction

We consider the following Darcy problem:

where Ω is the plane square ] — 1, 1[^{2} and n = (n 1 ,n2) is the exterior unit normal to the boundary Γ = ∂Ω. The boundary Γ is divided into two parts: the horizontal portion Γ^{1} = {(x,y)| — 1 < x < 1 ,y = ± 1} and the vertical portion Γ^{2} = {(x,y)|x = ± 1, — 1 < y < 1}. As we can see, the boundary condition on Γ^{2 }are nonstandard, since we prescribe the value of the pressure on Γ^{2}. On the contrary, we have a classical condition on the portion Γ^{1}. These conditions are described in the following figure.

The
equations of the Darcy problem not only modelize the flow in porous media, but
also appear in the projection techniques for the solution of Navier-Stokes
equations (see [10] and [15]). The nonstandard boundary conditions, where the
pressure is assigned on a part of the boundary, have a physical meaning:
typically the portion Γ^{1} corresponds to rigid
walls, whereas the entry or exit of the fluids takes place through Γ^{2}. The spectral discretization with this type of boundary conditions
was only studied within the framework of the Stokes problem (see [6], [4] and
[5]), while Azaiez, Bernardi and Grundmann proposed in [2] the spectral
discretization of the standard Darcy problem, where the normal velocity is
assigned on the boundary.

This paper is devoted to the spectral discretization of the nonstandard Darcy problem. First, we give two variational formulations. Each one leads us to well-posed problems. Second, we study the regularity of the solution by using a mixed problem of Dirichlet- Neumann for the Laplace operator. Next, from the first variational formulation, we derive a first spectral discretization, which is simple, but, in order to improve the approximation of the divergence-free condition, we study a second spectral discretization, where the inf-sup condition is obtained with more difficulty and where the error estimate on the pressure is not optimal. Finally, the second variational formulation yields a third spectral discretization, which leads us to optimal error estimate and a divergence-free discrete solution.

An outline of this paper is as follows. The two continuous variational problems, as well as the regularity of the solution, are studied in Section 2. Section 3 is devoted to the analysis of three spectral approximations of this problem in the case of homogeneous boundary conditions. In Section 4, we present the algorithms that are used to solve the first and third discretizations, together with some numerical experiments.

Statement of the problem and notation

In order to set this problem into
adequate spaces, recall the definition of the following standard Sobolev spaces
(cf. J. Necas [13]). For any multi-index k = (k_{1} ,k_{2}) with k_{i} > 0, set |k| = k_{1} + k_{2} and denote

Then for any integer m ≥ 0 and any plane domain Ω whose boundary is Lipschitz continuous(cf. Grisvard [12]), we define:

equipped with the seminorm

and norm(for which it is an Hilbert space)

For extensions of this definition to non-integral values of m (see [11,12]), let s a real
number such that s = m + σ with m ∈ N and 0 <σ< 1. We denote by H^{s}(Ω) the space
of all distributions u defined in Ω such that u ∈ H^{m}(Ω) and, ∀|α| = m,

It can be shown that H^{s}(Ω) is a Hilbert space for the scalar product

Let Γ' be an open part of the boundary ∂Ω of class Cm−1,1 and ${T}_{1}^{\mathrm{\Gamma \text{'}}}$ the mapping $\upsilon \to {\upsilon}_{|\Gamma \text{'}}$ defined on H^{m}(Ω). We denote by ${H}^{m-\frac{1}{2}}\left(\mathrm{\Gamma \text{'}}\right)$ (see [7,12]) the space ${T}_{1}^{\mathrm{\Gamma \text{'}}}$ (H^{m}(Ω))
which is equipped with the norm:

In this text, we shall use the spaces H^{1/2}(Γ') and H^{3/2}(Γ') corresponding to m = 1 and 2.
Let us define the space H^{1/2}_{00} (Γ') = {v_{|Γ'}, v ∈ H1(Ω), ∀x ∈ ∂Ω \ Γ ', v_{|∂Ω}(x)=0}. We
shall also be interested in the dual space of H^{1/2}_{00} (Γ'),

We shall use the Hilbert space H(div ; Ω) = {v ∈ L^{2}(Ω)^{2} : div v ∈ L^{2}(Ω)}, equipped
with the norm

For vanishing boundary values, we define:

For Λ =]−1,1[ , we denote the norm in L^{2}(Λ) by${\Vert \upsilon \Vert}_{0,\mathrm{\Lambda}}={\left({\int}_{-1}^{1}{\upsilon \left(x\right)}^{2}\phantom{\rule{0.5em}{1em}}dx\right)}^{\frac{1}{2}}$ , the seminorm in H^{1}(Λ)and the norm in H^{1}(Ω) by ${\Vert \upsilon \Vert}_{1,\mathrm{\Lambda}}={\left({\int}_{-1}^{1}{\left(\upsilon \right(x\left)\right)}^{2}+{(\upsilon \text{'}(x\left)\right)}^{2}\phantom{\rule{0.5em}{1em}}dx\right)}^{2}$ .

We note x = (x,y) the generic point of the square Ω and we call Γ_{I}, Γ_{II }, Γ_{III }and Γ_{IV}_{}the edges of Ω, starting from west and turning counterclockwise. For each J, J = I,II,III,IV, the extremities of the edge Γ_{J} are a_{J}-1 and a_{J}, with the convention a_{o} = a_{IV}, the exterior unit normal vector to Γ_{J}is denoted by n_{J} and the counterclockwise unit tangent vector is ${\tau}_{J}$ . Figure 1.2 below presents this notation.

For any domain Δ in R or R^{2} and for any nonnegative integer n, P_{n}(Δ) stands for the space of all polynomials on Δ with degree < n with respect to each variable. We also use the notation P^{0}_{n}(Δ) for the subspace P_{n}(Δ) ∩ H^{1}_{0}(A). For ∆ =]-1, 1[, the family (Ln)n of Legendre polynomials is a basis of the spaces P($\textcolor{red}{}\mathrm{\Lambda}$) of polynomials on $\textcolor{red}{}\mathrm{\Lambda}$ (we refer to [7, Chap. I] for the properties of the orthogonaI poIynomiaIs). These poIynomiaIs are orthogonal to each other in L^{2}(
$\textcolor{red}{}\mathrm{\Lambda}$
) and are caracterized as follows: for any integer n > 0, the polynomial L_{n} is of degree n and satisfies L_{n}(1) = 1. Let us recall some properties that we need. The family (L_{n})_{n} is given by the recursion relation:

Each polynomial is a solution of the differential equation

and its norm is given by

From (2.6) and integration by parts, we derive

Next, let N ≥ 2 be a fixed integer. We denote by ξj , 0 ≤ j ≤ N, the zeros of the polynomial (1 − ζ2)L N (ζ) in increasing order. We recall (see [7, Chap. I]) that there exist positive weights ρj , 0 ≤ j ≤ N, such that the following equality, called the Gauss-Lobatto formula, holds

Moreover, it follows from the identity (see [7, Chap. III])

that the bilinear form: $\textcolor{red}{}\left(u,\upsilon \right)\to \sum _{j=0}^{N}u\left({\xi}_{j}\right)\upsilon \left({\xi}_{j}\right){\rho}_{j}$ is a scalar product on P_{N} (Λ), since we
have

THE CONTINUOUS PROBLEM

First variational formulation

We define the subspace H ^{1}(Ω; Γ') = {v ∈
H
^{1}(Ω)|v = 0 on Γ'} of H ^{1}(Ω)
and we introduce the space:

In the same way as in [1], we consider the following equivalent variational formulation of the problem (1.1)-(1.4):

Find u in L^{2}(Ω)^{2}and p in H^{1}(Ω) such that p − ϕ belongs to M and that

where the bilinear forms a and b are defined by

Theorem 3.1Let fbe in L^{2}(Ω)^{2}, U_{0} in H^{−1/2}(Γ^{1}) and ϕ in H^{1/2}(Γ^{2}), where H^{−1/2}(Γ^{1})
and H^{1/2}(Γ^{2}) are defined respectively in (2.3) and (2.2). Then problem (3.2), (3.3) has a
unique solution satisfying

Proof. First, let us define φ ̃ belonging to H^{1/2}(Γ) such that
φ ̃|_{Γ}^{2} =
φ and ||
φ ̃||_{1/2,Γ} ≤
c||
φ ||_{1/2,Γ}^{2} . To this end, we must extend ϕ to a function belonging to H^{1/2}(Γ). Let µ be a
function defined in [0, 2] by

We define φ ̃|Γ_{11} by

and φ ̃|_{ΓIV }by

Then we have ∥φ ̃∥_{1/2,Γ} ≤ C∥φ∥_{1/2,Γ}^{2} . Next, let Φ in H^{1}(Ω) such that Φ|_{Γ} = φ ̃. Finally, we
obtain a function Φ verifying

Second, let us extend U_{0} to a function belonging to H^{−1/2}(Γ). We set U ̃_{0|Γ}^{1} = U_{0} and
U ̃_{0|Γ}^{2} = −1/2 < U_{0},1 >_{Γ}^{1} . Then, we have ∥U ̃_{0}∥_{−1/2,Γ} ≤ C∥U_{0}∥_{−1/2,Γ}^{1} and < U ̃_{0},1 >_{Γ}= 0.
Next, we define Neumann’s Problem:

and we set u_{0} = ∇ψ. Applying Proposition 1.2 of [11, page 14], we derive that ψ belongs to H^{1} (Ω), which implies that u_{0} belongs to H (div ; Ω) with

since div u_{o} = ∆ψ = 0. Finally, u_{0} verifies

Now, let us split p as: p = Φ + p ̃ with p ̃ in M and u as: u = u_{0} + u ̃. Then, we can
write the problem (3.2),(3.3) as

Since the right-hand side of (3.9) defines a continuous form on L^{2}(Ω)^{2} and since the properties of continuity and ellipticity are obvious we have only to check the following inf-sup
condition on the form b (see [11, pages 58,59]):

with a positive constant β. This ”inf-sup condition” was introduced independently by Babuska [3] and Brezzi [9]. We can verify this condition by taking v = ∇q . Indeed, we have

and, since q|_{Γ2} = 0, using a generalization of Poincar ́e inequality (see [11, Chap. I, page 40]) yields ∥q∥L^{2}_{(Ω)} ≤ P|q|H^{1}_{(Ω)}, which implies

Thus, the ”inf-sup condition” is verified with the positive constant β = 1/√ (P)2 + 1). Hence, applying Theorem 2.3 [7, pages 116,117], the theorem follows.

Second variational formulation

We introduce the space:

In an analogous way as in [1], we consider the following equivalent variational formulation of the problem (1.1)-(1.4):

Find u in X and p in L^{2}(Ω) such that u − u_{0} belongs to X, where u_{0} is the function
previously constructed that verifies (3.8), and that

where the bilinear form b∗ is defined by

In the same way as previously, we split u as: u = u_{0} + u ̃. Then, we can write the
problem (3.13),(3.14) as

Since the right-hand side of (3.16) defines a continuous form on X and since the properties of continuity and ellipticity are obvious, we have only to check the following inf-sup condition on the form b∗:

with a positive constant β∗. Let us note that belongs to L^{2}_{0}(Ω) and, owing to a classic result (see [11, Chap. I]), there exists v_{0 }in H^{1}_{0}(Ω)^{2}, such that

Then, we set

We can verify that v ̃ belongs to X with

since Then, we have

Hence, we derive the inf-sup condition and we obtain the following result.

Theorem 3.2Let fbe in L^{2}(Ω)^{2,} U_{0} in H^{−1/2}(Γ^{1}) and φ in H^{1/2}(Γ^{2}), where H^{−1/2}(Γ^{1})
and H^{1/2}(Γ^{2}) are defined respectively in (2.3) and (2.2). Then problem (3.13), (3.14) has
a unique solution satisfying

Regularity results

When the data f is in H (div ; Ω), taking the divergence of the first equation of the problem (1.1)-(1.4) and owing to the other equations, we obtain a mixed problem of Dirichlet-Neumann for the Laplace operator:

We suppose that f is in H^{1}(Ω)^{2}, φ is in H^{3/2}(Γ^{2}) and U_{0 }is in H^{1/2}(Γ^{1}). In addition, we
assume matching conditions at the vertices of Γ (see [7, Chap I]):

Theorem 3.3For any data f in H^{1}(Ω)^{2}, φ in H^{3/2}(Γ^{2}) and U_{0} in H^{1/2}(Γ^{1}), where
H^{3/2}(Γ^{2}) and H^{1/2}(Γ^{1}) are defined in (2.2) , verifying the matching conditions (3.21),
the solution (u, p) of the problem (1.1)-(1.4) belongs to H^{1}(Ω)^{2} × H^{2}(Ω).

Proof. Owing to matching conditions (3.21), there exists p_{0} in H^{2}(Ω) such that p_{0|Γ}^{2}_{ = φ} and (∂p0/∂n)|Γ^{1} =f.n−U_{0}. Let us set p ̃=p−p_{0}. The problem (3.20) is equivalent to the
following problem: find p ̃ in H^{1}(Ω; Γ^{2}) such that

Since the boundary between Γ^{1} and Γ^{2} is the set of vertices of Γ, the regularity of the data implies that this homogeneous mixed problem of Dirichlet-Neumann for the Laplace operator has a solution p ̃ in H^{2}(Ω) (see [12]). Hence, we derive the regularity of p and u. ♦

Remark 3.4If (f . n−U_{0}) is Lipschitz-continuous on Γ_{J} or belongs to H^{1}(Γ_{J }), J = II, IV and if φ belongs to C^{1,1} (Γ_{J)} or to H^{2} (Γ_{J}), J = I, III, the matching conditions (3.21) are equivalent to simpler conditions:

SPECTRAL DISCRETIZATION

First spectral discretization

We define the discrete scalar product by

and we denote by L_{N} the Lagrange interpolation operator at the points (ξi,ξj), 0 ≤ i,j ≤ N in lP_{N} (Ω). We set

We assume that the data f belongs to C^{0}(Ω)^{2} and, for the sake of simplicity, that u and p
satisfy homogeneous boundary conditions, that is to say, we set φ = 0, U_{0} = 0 in (1.1)-(1.4).
From the variational formulation (3.2)-(3.3), we derive the following discrete problem:

Find u_{N} in X_{N} and p_{N }in lP_{N}(Ω)∩M, where M is defined by (3.1), such that

where the form b_{N} is defined by

We have a classical saddle point problem. We verify the inf-sup condition

where γ is a positive constant independent from N , by taking v_{N} = ∇q_{N }. Hence, we derive
the following theorem.

Theorem 4.1Let f be in C^{0}(Ω)^{2}. Then problem (4.3), (4.4) has a unique solution
(u_{N}, p_{N} ) satisfying

Next, we establish a theorem which implies the convergence of our discretization method.

Theorem 4.2Assume that the solution (u, p) of problem (4.3), (4.4) belongs to H^{s}(Ω)^{2} × H^{s+1}(Ω), s ≥ 0, and the data f belongs to H^{σ}(Ω)^{2}, σ > 1, where H^{s}(Ω), for non-integral values of s, is defined in (2.1). Then, the following estimate holds

Proof. From the abstract error estimate for the approximation of saddle-point problems (see [7, Chap. IV]), we derive the following estimate:

where V_{N} is defined by

Moreover, we recall (see [11, Chap. II, (1.16)]) that

Hence, we derive that, for all v_{N−1} and f_{N−1} in lP_{N−1}(Ω)^{2} and all q_{N−1} ∈ lP_{N−1}(Ω) ∩ M,

Then, we choose v_{N−1} = II _{N−1}u (resp. f_{N−1}= II_{N−1}f), that is to say the orthogonal
projection of u (resp. f) on lP_{N-1}(Ω)^{2} in L^{2}(Ω)^{2} and q_{N-1} = II^{1,Γ}_{2}_{N-1}p, where II^{1,Γ}_{2 N-1}p is the orthogonal projection of p on lP_{N−1}(Ω) ∩ M in H^{1}(Ω). It remains to prove the estimate,
for any m ≥ 1,

On the one hand, this result is obvious for m = 1. On the other hand, for m ≥ 2, we have (see [7, Chap. III]),

Then, an interpolation argument (see [11, TH 1.4, page 6]) gives (4.11). Finally, the result follows from (4.10), (4.11) and the classic estimate for the orthogonal projection on
lP_{N−1}(Ω) in L^{2}(Ω).

Remark 4.3 With the choice X_{N} = lP_{N}(Ω)^{2}, problem (4.3, (4.4) can be interpreted as a
collocation scheme. Indeed, by integrating by parts in the discrete bilinear form b_{N} with
respect to one of the two variables for each of the two terms of bN (this process being
allowed by the precision of the quadrature rule), and choosing as test functions the Lagrange
polynomials associated with the grid points of Ξ_{N}, it is easily seen that (4.3), (4.4) is
equivalent to these to set of equations for u_{N} in lP_{N}(Ω)^{2} and p_{N }in lP_{N}(Ω) ∩ M:

Second spectral discretization

In order to improve the approximation of the condition div u = 0, we can try to
decrease the dimension of the space X_{N} . So, we choose

We note that, in this case the forms b(., .) and b_{N} (., .) are equal on X_{N} × lP_{N} (Ω). It does
not appear spurious modes for the pressure, as we can see in the following lemma.

Lemma 4.4 Let Z_{N} be the space

Then Z_{N} = {0}.

Proof. Let q_{N} be in Z_{N }. Since ∂q_{N}/∂x is a polynomial of degree ≤ N − 1 with respect to x, which is orthogonal to lP_{N−1}(Ω), we can write:

with α_{N} and β_{N} in lP_{N }(Ω). In the same way with y in place of x, we have:

with γ_{N} and δ_{N} in lP_{N} (Ω). Hence, we derive

where λ and _{μ} are real numbers. But, since L_{N} (1) = 0, the condition:

implies λ = μ = 0.

We have to study the following discrete problem:

Find u_{N} in lP_{N−1}(Ω)^{2} and p_{N} in lP_{N}(Ω)∩M such that

For the inf-sup condition, the choice v_{N }= ∇q_{N} is no longer available, since v_{N }must be in lP_{N−1}(Ω)^{2}. In fact, in the next lemma, we find a inf-sup constant depending on N.

Lemma 4.5 There exists a constant c > 0 independent on N such that

On the one hand, we note that

where a_{N} is a polynomial of degree ≤ N=−1. Then, we have

which implies, using (2.7) and the inverse inequality (2.10),

On the other hand, in view of (2.8), we can write

where r_{N} (x, y) is a polynomial of lP_{N} (Ω) of degree < N − 1 with respect to y. The orthogonality properties imply

By combining this inequality with (4.18), we obtain

In the same way, we have the analogous inequality for ∂q^{*}_{N}/∂_{y} . Thus, we obtain

Next, the equality (4.16) yields

which implies, owing to (2.7) and (2.9),

It remains to estimate |α_{N} |. First, we have, owing to (4.16),

But, ∀y ∈ [−1, 1], q_{N} (1, y) = 0, which implies

If we set , we derive

and, therefore, thanks to a discrete Cauchy-Schwarz inequality

Then, in view of (4.19), we obtain

Finally, (4.20), (4.21) and (4.22) yield

which, in view of (4.17), ends the proof.

The bilinear form (.,.)_{N} and b(.,.) satisfy Brezzi’s conditions with respect to lP_{N−1}(Ω)^{2}
and lP_{N }(Ω) ∩ M (the bilinear form (., .)_{N} is continuous on lP_{N −1} (Ω)^{2 }and elliptic on
lP_{N−1}(Ω), the bilinear form b(.,.) is continuous on lP_{N−1}(Ω) × (lP_{N}(Ω) ∩ M) and verifies
the “ inf-sup condition”), see [7, Theorem 2.3, pages 116,117], whence the theorem.

Theorem 4.6Let f be in C ^{0} (Ω)^{2} . Then problem (4.13), (4.14) has a unique solution
(u_{N} , p_{N} ) satisfying

We establish the convergence of this second discretization in the following theorem.

Theorem 4.7Assume that the solution (u, p) of problem (4.13), (4.14) belongs to H^{s}(Ω)^{2}
× H^{s+1}(Ω), s ≥ 0, and the data f belongs to H^{σ(}Ω)^{2}, σ > 1. Then, the following estimate
holds

Proof. The abstract error estimate, analogous to (4.39) but with much simplification because of the exactness of the quadrature formulas, yields

where V_{N} is now the space

It remains to estimate the term $\underset{{W}_{N}\mathrm{\in}{V}_{N}}{\text{i}\text{n}\text{f}}$∥u − w_{N} ∥L^{2}(Ω)^{2} . Since u is such that div u = 0 and
u . n|Γ1 = 0, there exist a unique ψ in H1(Ω) (see [11, Chap. I and 4]) such that

u = curl ψ and ψ = 0 on Γ^{1}.

Moreover, if u belongs to H^{s}(Ω)^{2}, we have
∥ψ∥_{H}^{s+1}_{(Ω)} ≤ c∥u∥_{H}^{s}_{(Ω)}^{2}.
Let us define the operator π ̃_{N}^{1} (see [7, Chap. II]) on H^{1}(Λ) by

where the function φ ̃ stands for

Note that the definition of π ̃^{1}_{N} is available, because φ ̃ belongs to H^{1}_{0}(Λ), and that π^{1,0}_{N}φ ̃ and φ coincide in −1 and 1. In [8, Section 7], the following estimate is proven, for all r ≥ 1
and all function φ in H^{r}(Λ):

Assuming s ≥ 1, we set

Since we can verify that R_{N−1} (u) belongs to V_{N} and, in view of (4.26), that

Finally, we derive, for s ≥ 1,

Since we have an interpolation argument gives the result of approximation in V_{N} for any s ≥ 0 and the estimate of the theorem follows.

Third spectral discretization

The third discretization comes from the variational formulation (3.13), (3.14). We
define the space X_{N} by

Let M_{N} be a subspace of lP_{N}( Ω) that we shall set later. Then, we consider the following
discrete problem:

Find u_{N} in X_{N}_{}and p_{N} in M_{N }such that

where the form b∗_{N} is define by

In order to choose M_{N} , we begin to identify the spurious modes for the pressure. These
spurious modes for the pressure are derived by elimination from those of classic Stokes
problem (see [7, Chap. IV]). In particular, we can verify: ∀v_{N} ∈ X_{N} , b∗_{N} (v_{N} , q_{N} ) = 0, for
q_{N} (x, y) = L_{N} (x) or L_{N} (x)L_{N} (y). We obtain the following lemma.

Lemma 4.8 Let Z_{N}^{∗} be the space

Then Z_{N}^{∗} is spanned by (L_{N} (x), L_{N} (x)L_{N }(y)).

Finally, let M_{N} stand for the orthogonal complement of Z_{N}^{∗} for the scalar product in
L^{2}(Ω) or for the scalar product (.,.)N, owing to (2.11). The inf-sup condition is given in
the next lemma.

Lemma 4.9There exists a constant c > 0 independent from N such that

Proof. Any function q_{N} in M_{N}_{}has the expansion

With the convention L_{−1} = 0, we choose w_{N} = (w_{N} , z_{N} ) with

and

Then, in view of (2.8), we have

since w_{N} .n_{|Γ}^{1} = z_{N|Γ}^{1} =0. As in[6,Chap. IV], we prove

Hence, we derive, in the same way as in [8, Section 24]

Next, setting w_{N}∗ (x, y) = w_{N} (x, y) − q_{0,0}L_{1}(x) − q_{0,N} L_{1}(x)(L_{N} (y) − L_{N−2}(y)), we note that
w_{N}^{∗} (±1,y) = 0 for −1 ≤ y ≤ 1. Then, the Poincar ́e-Friedrichs inequality, applied with
respect to x or y, yields

Hence, owing to (4.35) and the estimate

which is derived from (4.34), we obtain

Finally (4.35) and (4.36) imply

and, in view of (4.33), the choice v_{N} = w_{N} in b∗_{N} (v_{N} , q_{N} ) is available and gives the inf-sup
condition (4.30).

From the previous lemma, we derive the following theorem.

Theorem 4.10Let f be in C^{0}(Ω)^{2}. Then problem (4.27), (4.28) has a unique solution
(u_{N} , p_{N} ) satisfying

Sketch of the proof. From(4.27), we derive (u_{N} , u_{N} )_{N} = (L_{N}f , u_{N} )_{N} . Then, Owing to
div u_{N} = 0 and (2.13), we derive ∥u_{N} ∥_{H (div ;Ω)} ≤ 3∥L_{N}f ∥_{L2 (Ω)} . Next, the inf-sup condition
(4.30) and (4.27) imply

Hence, in view of (2.13), we obtain

which ends the proof.

Next, in the same way as in Theorem 4.2, we prove an optimal error estimate.

Theorem 4.11 Assume that the solution (u,p) of problem (4.27), (4.28) belongs to H^{s}(Ω)^{2}
×H^{s}(Ω), s ≥ 0, and the data f belongs to H^{σ}(Ω)^{2}, σ > 1. Then, the following estimate holds.

Sketch of the proof. Again, from the abstract error estimate for the approximation of saddle-point problems (see [7, Chap. IV]), we derive the following estimate:

where V_{N} is defined by

Moreover, we still have (see [11, CH. II, (1.16)])

We end the proof in the same way as in Theorem 4.2.

Remark 4.12 As for the first discretization, problem (4.27, (4.28) can be interpreted as
a collocation scheme. In the same way, we prove that (4.27), (4.28) is equivalent to the
set of equations for u_{N} in X_{N }and p_{N} in M_{N}:

Therefore, the discrete solution u_{N} is exactly divergence-free, which is important for
some applications.

Numerical Results

The convergence of the methods corresponding to the first and third discretizations
were tested in a problem of the type (1.1)-(1.4), with homogeneous boundary conditions.
Precisely, we tested the convergence of these methods to the exact solution u(x,y) =
(πx^{2} cos(πy), −2x sin(πy) and p(x, y) = y sin(πx), which means that we studied the convergence of these methods for Problem (1.1)-(1.4), with

and homogeneous boundary conditions. In addition, we tested the convergence to 0 of the divergence for both methods.

We shall use the Lagrange polynomials. We denote l_{r} the Lagrange polynomial associated to the Gauss-Lobatto point ξ_{r}, 0 ≤ r ≤ N, the expression of which is

The derivative l_{r}′ verifies the following equalities

Uzawa's algorithm

Problem (4.3), (4.4), Problem (4.13), (4.14) and Problem (4.27), (4.28) are equivalent to a linear system of the type :

The unkowns are the vectors U and P which represent respectively the velocity and the
pressure values on a given grid points. The data f is representated by the vector F on
the same grid points. The diagonal matrix M is the weight matrix, while the matrix D is
associated to the form b_{N} for the first and second spectral discretizations and to the form
b^{∗}_{N} for the third spectral discretization and D^{T} is the transposed matrix of D.

Uzawa’s algorithm consists in rewriting the first equation of system (5.4) as:
U = F − M^{−1}DP and substituting in the second equation. We obtain a new equation for the pressure P:

Next, we solve this symetric system either directly if the matrix D^{T} M^{−1}D is invertible or
by diagonalizing the matrix D^{T}M^{−1}D if not, because the spurious modes correspond to
eigenvalues of the matrix equal to zero. Next, we compute the velocity via the formula

and its divergence by multiplying U on the left by the matrix D^{T} . Finally, we test the
convergence to 0 of its divergence.

Implementation of the rst discretization

We take (l_{j}(x)l_{k}(y),0), (0,l_{j}(x)l_{k}(y)), 0 ≤ i,j ≤ N as a basis of X_{N} and
l_{r}(x)l_{s}(y),1≤r≤N−1,0≤s≤N as a basis of P_{N}(Ω)∩M. The unknowns are
the velocity u_{N} = (u^{1}_{N},u^{2}_{N}) and the pressure p_{N}. For j = 1,2, for 0 ≤ r,s ≤ N, we
denote,where f=(f_{1},f_{2}) represents the data, and, for 1≤r≤N−1, for 0≤s≤N, we denote P^{N}_{r,s}_{}=p_{N} (ξ_{r},ξ_{s}). So, we have

Let us define the (N + 1)^{2} × (N + 1)^{2} diagonal matrix with

and the 2(N + 1)^{2} × 2(N + 1)^{2} diagonal matrix

By setting v_{N}(x,y) = (l_{j}_{(}x)l_{k}(y),0), for 0 ≤ j,k ≤ N and, next, v_{N}(x,y) = (0,l_{j}(x)l_{k}(y)
in (4.3), (4.4), we obtain the matrix system (5.4) where the column matrices U, P and F
are defined by

and, with (.,.)_{N} defined by (4.1), the 2(N +1)^{2} ×(N^{2} −1) matrix D by

where (j, k) represents the row index, (r, s) the column index, 0 ≤ i, j, s ≤ N ,
1≤r≤N−1, for the matrices D^{1,N} and D^{2,N}. Note that U and F are 2(N+1)^{2}×1
matrices and P is a (N^{2} − 1) × 1 matrix.

Proposition 5.1 The (N^{2} − 1) × (N^{2} − 1) square matrix D^{T} M^{−1}D is invertible.

Proof. We just have to prove that the rank of the matrix D is N^{2} − 1. Let us assume
that the rank of the matrix D is strictly smaller than N^{2} − 1. Then, there exist a sequence of real number (q_{r,s}),1≤r≤N−1,0≤s≤N, where all the real numbers q_{r,s} are not equal to zero, such that

where Dr,s are the column vectors of the matrix D. Setting the previous equality is equivalent to

which is in contradiction with the property that there is no spurious mode.

We have to compute the matrix D^{T}M^{−1}D = (b_{(t,u),(r,s)}), 1 ≤ r, t ≤ N −1, 0 ≤ s, u ≤ N.

Owing to the previous expression of the matrices D and M, we obtain

Next, we change the numbering for the matrix D^{T} M^{−1}D. Let us define the mapping φ by

Note that φ is a one to one mapping from {1,...,N −1}×{0,...,N} to {1,...,N^{2} −1}
and we note

Note that
ψ_{1}(i)-1 and ψ_{2}(i) are respectively the quotient and the remainder of the euclidian division of i − 1 by N + 1. Then, we can denote

where (t, u) = (ψ_{1}(i), ψ_{2}(i)) and (r, s) = (ψ_{1}(j), ψ_{2}(j)).

Computing the elements a_{i,j} of the matrix D^{T}M^{−1}D yields

Thus, most of elements of the matrix D^{T} M^{−1}D are equal to zero.

Next, we determine the column matrix D^{T} F = (c_{i}_{,1}) by

Hence, owing to (5.2) and (5.3), we compute the list [l_{k}′(ξ_{m}), 0 ≤ k,m ≤ N], which
allows us to determine the elements of the matrices D^{T} M^{−1}D and D^{T} F. Since the matrix
D^{T} M^{−1}D is invertible, the equation (D^{T} M^{−1}D)X = D^{T} F has a unique solution X. Then
we derive easily the column matrix P =(P^{N}_{r,s})such that P^{N}_{r,s} =Xφ(r,s)
,1≤r≤N−1,
0 ≤ s ≤ N and, next, the column matrix U, thanks to the relation U = F − M^{−1}DP.
Setting P^{N}_{0,t} = P^{N}_{N,t} = 0 for 0 ≤ t ≤ N in accordance with the boundary conditions, we obtain

Finally, we derive

We give the values of (div u_{N}, div u_{N})N^{1/2} , (p−p_{N},p−p_{N})N^{1/2} and (u−u_{N},u−u_{N})N^{1/2} for
N between 4 and 21.

We shall comment these results on comparing them with these ones of the third discretization.

Implementation of the first discretization

We only sketch the method because the first and third discretizations are rather similar, but we shall point out the differences between both discretizations.

First, the space X_{N} is not the same, because boundary conditions are taken in account.
Thus, we take (l_{j}(x)l_{k}(y),0), (0,l_{s}(x)l_{r}(y)), 0 ≤ i,j,s ≤ N, 1 ≤ r ≤ N − 1 as a basis of
X_{N}.

Second, we shall compute the pressure p_{N} as the orthogonal projection on M_{N} of an
element of lP_{N}(Ω) and we take l_{r(}x)l_{s}(y), 0 ≤ r,s ≤ N as a basis of lP_{N}(Ω). For 0 ≤ r,s ≤
N, we denote U^{1,N}_{r,s} =u^{1}_{N} (ξ_{r},ξ_{s}),P^{N}_{r,s} =p_{N}(ξ_{r},ξ_{s}) and for 0≤r≤N,1≤s≤N−1, we denote U^{2,}^{N}_{r,s} = u^{2}_{N} (ξ_{r} ,ξ_{s} ). Let us define the matrix M ̃_{1} = M ̃, where M ̃ is given by (5.7), the(N^{2} −1)×(N^{2} −1)matrix M ̃_{2} =(m_{(j,k),(r,s)}) with

and the 2N(N + 1) × 2N(N + 1) diagonal matrix

In the same way as the first discretization, we derive the following matrix system

where the column matrices U, P and F are defined by

and the 2N(N +1)×(N +1)^{2} matrix D∗ by

We have to compute the matrix D^{T}∗M^{−1}∗D∗ =(b∗
_{(t,u),(r,s)} ),0≤r,s,t,u≤N. In the same way as the first discretization, we obtain

Next, as in the first discretization, we change the numbering for the matrix D^{T}∗M^{−1}
∗D∗.
Let us define the mapping φ∗ by

Note that φ∗ is a one to one mapping from {0,...,N}^{2} to {1,...,(N +1)^{2}} and we note

Note that ψ_{1}∗(i) and ψ_{2}∗(i) are respectively the quotient and the remainder of the euclidian
division of i−1 by N +1. Then, we can denote

where (t, u) = (ψ_{1}∗(i), ψ_{2}∗(i)) and (r, s) = (ψ_{1}∗(j), ψ_{2}∗(j)).

Computing the elements a^{∗}_{i,j }of the matrix D^{T}∗ M^{−1}∗D∗ yields

Next, we determine the column matrix
D∗^{T} F = (c^{∗}_{i,1}) by

Now, we deal with the main difference between both discretizations. In the third discretization, the matrix D^{T}_{*}M^{-1}_{*}D_{*} is not invertible and the computation of the pressure is more complicated. First, let be a spurious mode, that is to say
an element of the space Z_{N}^{∗} defined in Lemma 4.8. In the same way as in the proof of
Proposition 5.1, considering that the matrices D^{T}_{*}M^{−1}_{*}D_{*} and D_{*} have the same rank, we
obtain the following equivalences

where Q is the column vector (q
_{t,u}) _{0≤t,u≤N} . We can consider the matrix D^{T}_{*}M^{−1}_{*}D_{*} as the matrix of a linear mapping f from the vector space lP_{N}(Ω) into itself equipped with the
basis B = (l_{i}(x)l_{j}(y))_{0≤i,j≤N} .Therefore, Z_{N}^{∗} is the eigenspace associated to the eigenvalue
equal to 0 of the linear mapping f or, equivalently, of its matrix D^{T}_{*}M^{−1}_{*}D_{*} . Note that this basis is orthonormal for the scalar product (.,.)_{N}. We can diagonalize the positive symmetric matrix D^{T}_{*}M^{−1}_{*}D_{*} and, thus, there exist a diagonal matrix Λ and an orthogonal matrix R such that D^{T}∗ M^{-1}∗ D∗ = RΛR^{-1} and such that the diagonal elements of Λ, that is to say (λ_{i,i})_{1≤i≤(N+1)}^{2} are in increasing order. Note that the matrix Λ is the matrix of f in
a basis B′ = (h_{i})_{1≤i≤(}N_{+1)}^{2} of lP_{N} (Ω), which is orthonormal for the scalar product (., .)_{N} , the elements of which are the eigenvectors of the matrix D^{T}_{*}M^{−1}_{*}D_{*} . Let P′ be the column matrix the elements of which are the components of the pressure p_{N} in the new basis B'
of lP_{N}(Ω). We have P′ =R^{T}P =R^{−1}P and the following equivalence

Since Z_{N}^{∗} is a two dimensions space, B_{0}′ = (h_{1} , h_{2} ) is the basis of Z_{N}^{∗} and (h_{i} )_{3≤i≤(N +1)2} is
a basis of M_{N} = (Z_{N}^{∗} )^{⊥}. Therefore, we determine the column vectors P′ and P by

and we have In the same way as for the first discretization, we derive the column matrix U by the relation U = F − M^{−1}_{*}D_{*}P, which gives, for
0 ≤ j, k, t ≤ N , 1 ≤ u ≤ N − 1,

Finally, considering the boundary conditions, we set U^{2,N}_{t,0} = U^{2,}^{N}_{t,N} = 0, 0≤t≤N , and we obtain the same formulas as in the first discretization for (div u_{N} , div u_{N})_{N} ,
(u − u_{N} , u − u_{N})_{N} and (p - p_{N}, p - p_{N})_{N} . We also give the values of ( div u_{N} , div u_{N} )_{N}^{1/2} ,
(p−p_{N} ,p−p_{N})_{N}^{1/2} and (u- u_{N}, u - u_{N})_{N}^{1/2} for N between 4 and 21.

Conclusion

Continuous problem and discretizations

We studied the Darcy problem by defining two equivalent variational formulations
corresponding to two different bilinear forms b and b^{∗}. The first one requests a more regular
pressure, which is bounded in H^{1}(Ω). The second one is less classic and was introduced
by A. Quarteroni and A. Valli (see [14]). This second variational formulation is important
because it allowed us to constuct a spectral method which gives a divergence-free discrete
solution.

Moreover, by studying a mixed problem of Dirichlet-Neumann for the Laplace operator,
we proved regularity results with the solution (u,p) in H^{1}(Ω)^{2} ×H^{2(}Ω) as long as the data
is regular enough.

The first variational formulation led us to two spectral discretizations. In the first one, the discrete solution is not divergence-free, which is a disadvantage for some applications. However, we obtain a fully optimal error estimate for the velocity and the pressure. In the second discretization, it does not appear spurious modes for the pressure. Moreover, because of the exactness of the quadrature formulas, the error estimate is easier to obtain. However, the error estimate for the pressure is not optimal, because the inf-sup constant depends on N.

The second variational formulation led to a third spectral discretization. There are spurious modes, which complicate the study, but the discrete velocity is divergence-free, which is important when the system is a stage of solving of a time-dependent problem, and the error estimate for the velocity and the pressure is fully optimal. In conclusion, this third discretization is the best discretization and we will only use it hereafter.

Comparison of both spectral discretizations

Tables 4.1 and 4.2 test the convergence of respectively the first and third discretizations to the exact solutions. For the first discretization, we see the convergence to zero of (div u_{N}, div u_{N})_{N}^{1/2} , about 10^{−1} for N = 8 and about 10^{−12} for N = 19. Concerning the third discretization, we see that the quantity (div u_{N} , div u_{N} )_{N}^{1/2} is very small for all values
of N, about 10^{−16 }for N = 4 and about 10^{−12} for N = 20. Thus, we verify that, in the
third discretization, the discrete solution u_{N} is exactly divergence-free, which is important,
as we saw previously. This property appears clearly in Figure 4.3, which represents the logarithm to the basis 10 of the quantity (div u_{N} , div u_{N} )_{N}^{1/2} as a function of N , for even N
from 4 to 20, for both discretizations.

Regarding the velocity, the discrete solution u_{N} converges fast to the exact solution u
for both discretizations. For the pressure, we also obtain fast convergence, except for odd
N in the third discretization.

References

Achdou, Y., Bernardi, C. & Coquel, F. “A priori a posteriori analysis of finite volume discretizations of Darcy’s equations”, Numer. Math. 96 (2003), 17-42.

Azaïez, M., Bernardi, C. & Grundmann, M. “Méthodes spectrales pour les équations du milieu poreux”, East West Journal of Numerical Analysis, No. 2, 1994.

Babuska, I. "The finite element method with lagrangian multipliers”, Numer. Math. 20, pp. 179-192 (1973).

Bègue, C., Conca, C., Murat, F. & Pironneau, O. “Les équations de Stokes et de Navier-Stokes avec des conditions aux limites sur la pression”, Nonlinear Partial Differen- tial Equations and their Applications, Coll'ege de France Seminar, vol. IX (1988).

Bernard, J.M. "Spectral discretizations of the Stokes equations with non standard boundary conditions”, Journal of Scientific Computing, 20(3), 355-377 (2004).

Bernardi, C., Canuto, C. & Maday, Y. "Spectral approximations of the Stokes equa- tions with boundary conditions on the pressure”, SIAM J. Numer. Anal. 28, No 2 (1991), p. 333-362.

Bernardi, C. & Maday, Y. Approximations spectrales de problèmes aux limites ellip- tiques, Springer-Verlag, Paris, 1992.

Bernardi, C. & Maday, Y. Spectral Methods, Handbook of Numerical Analysis, Vol. V, Paris (1997).

Brezzi, F. "On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers”, R.A.I.R.O., Anal. Mumer. R2, PP. 129-151 (1974).

Chorin, A.J. " Numerical solution of the Navier-Stokes equations”, Math. Comput. 22 (1968), 745-762.

Girault, V. & Raviart, P. A. Finite Element Methods for Navier-Stokes Equations,Springer-Verlag, Berlin (1986).

Grisvard, P. "Elliptic Problems in Nonsmooth Domains”, Pitman Monographs and Studies in Mathematics, 24, Boston (1985).

Necas, J. les Méthodes Directes en Théorie des Equations Elliptiques, Masson, Paris (1967).

Quarteroni, A. & Valli, A. Domain decomposition methods for partial differential equations, Oxford science publication (1999).

Temam, R. "Une méthode d’approximation de la solution des équations de Navier- Stokes”, Bull. Soc. Math. France 98 (1968), 115-152.