Introduction

We consider the following Darcy problem:

where Ω is the plane square ] — 1, 1[² and n = (n 1 ,n2) is the exterior unit normal to the boundary Γ = ∂Ω. The boundary Γ is divided into two parts: the horizontal portion Γ¹ = {(x,y)| — 1 < x < 1 ,y = ± 1} and the vertical portion Γ² = {(x,y)|x = ± 1, — 1 < y < 1}. As we can see, the boundary condition on Γ² are nonstandard, since we prescribe the value of the pressure on Γ². On the contrary, we have a classical condition on the portion Γ¹. 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 Γ¹ corresponds to rigid walls, whereas the entry or exit of the fluids takes place through Γ². The spectral dis­cretization 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. Se­cond, 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 discretiza­tion, 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₁ ,k₂) with k_i > 0, set |k| = k₁ + k₂ 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₀₀ (Γ') = {v_|Γ', v ∈ H1(Ω), ∀x ∈ ∂Ω \ Γ ', v_|∂Ω(x)=0}. We shall also be interested in the dual space of H^1/2₀₀ (Γ'),

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

For vanishing boundary values, we define:

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

We note x = (x,y) the generic point of the square Ω and we call Γ_I, Γ_II , Γ_III and Γ_IVthe 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 Γ_Jis 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² 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⁰_n(Δ) for the subspace P_n(Δ) ∩ H¹₀(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²( $\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 ¹(Ω; Γ') = {v ∈ H ¹(Ω)|v = 0 on Γ'} of H ¹(Ω) 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²(Ω)²and p in H¹(Ω) such that p − ϕ belongs to M and that

where the bilinear forms a and b are defined by

Theorem 3.1Let fbe in L²(Ω)², U₀ in H^−1/2(Γ¹) and ϕ in H^1/2(Γ²), where H^−1/2(Γ¹) and H^1/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 φ ̃|_Γ² = φ and || φ ̃||_1/2,Γ ≤ c|| φ ||_1/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 φ ̃|Γ₁₁ by

and φ ̃|_ΓIV by

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

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

and we set u₀ = ∇ψ. Applying Proposition 1.2 of [11, page 14], we derive that ψ belongs to H¹ (Ω), which implies that u₀ belongs to H (div ; Ω) with

since div u_o = ∆ψ = 0. Finally, u₀ verifies

Now, let us split p as: p = Φ + p ̃ with p ̃ in M and u as: u = u₀ + 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²(Ω)² 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²_(Ω) ≤ P|q|H¹_(Ω), 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²(Ω) such that u − u₀ belongs to X, where u₀ 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₀ + 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²₀(Ω) and, owing to a classic result (see [11, Chap. I]), there exists v₀ in H¹₀(Ω)², 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, U₀ in H^−1/2(Γ¹) and φ in H^1/2(Γ²), where H^−1/2(Γ¹) and H^1/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¹(Ω)², φ is in H^3/2(Γ²) and U₀ is in H^1/2(Γ¹). In addition, we assume matching conditions at the vertices of Γ (see [7, Chap I]):

Theorem 3.3For any data f in H¹(Ω)², φ in H^3/2(Γ²) and U₀ in H^1/2(Γ¹), where H^3/2(Γ²) and H^1/2(Γ¹) 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¹(Ω)² × H²(Ω).

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

Since the boundary between Γ¹ and Γ² 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²(Ω) (see [12]). Hence, we derive the regularity of p and u. ♦

Remark 3.4If (f . n−U₀) is Lipschitz-continuous on Γ_J or belongs to H¹(Γ_J ), J = II, IV and if φ belongs to C^1,1 (Γ_J) or to H² (Γ_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⁰(Ω)² and, for the sake of simplicity, that u and p satisfy homogeneous boundary conditions, that is to say, we set φ = 0, U₀ = 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⁰(Ω)². 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(Ω)² × H^s+1(Ω), s ≥ 0, and the data f belongs to H^σ(Ω)², σ > 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(Ω)² and all q_N−1 ∈ lP_N−1(Ω) ∩ M,

Then, we choose v_N−1 = II _N−1u (resp. f_N−1= II_N−1f), that is to say the orthogonal projection of u (resp. f) on lP_N-1(Ω)² in L²(Ω)² and q_N-1 = II^1,Γ_2N-1p, where II^1,Γ_{2 N-1}p is the orthogonal projection of p on lP_N−1(Ω) ∩ M in H¹(Ω). 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²(Ω).

Remark 4.3 With the choice X_N = lP_N(Ω)², 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(Ω)² 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(Ω)² 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(Ω)². 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(Ω)² and lP_N (Ω) ∩ M (the bilinear form (., .)_N is continuous on lP_{N −1} (Ω)² 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 ⁰ (Ω)² . 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(Ω)² × H^s+1(Ω), s ≥ 0, and the data f belongs to H^σ(Ω)², σ > 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²(Ω)² . 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 Γ¹.

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

where the function φ ̃ stands for

Note that the definition of π ̃¹_N is available, because φ ̃ belongs to H¹₀(Λ), 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_Nand 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²(Ω) 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_Nhas the expansion

With the convention L₋₁ = 0, we choose w_N = (w_N , z_N ) with

and

Then, in view of (2.8), we have

since w_N .n_|Γ¹ = z_N|Γ¹ =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,0L₁(x) − q_0,N L₁(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⁰(Ω)². 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_Nf , u_N )_N . Then, Owing to div u_N = 0 and (2.13), we derive ∥u_N ∥_{H (div ;Ω)} ≤ 3∥L_Nf ∥_{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(Ω)² ×H^s(Ω), s ≥ 0, and the data f belongs to H^σ(Ω)², σ > 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² 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⁻¹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⁻¹D is invertible or by diagonalizing the matrix D^TM⁻¹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¹_N,u²_N) and the pressure p_N. For j = 1,2, for 0 ≤ r,s ≤ N, we denote,where f=(f₁,f₂) 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)² × (N + 1)² diagonal matrix with

and the 2(N + 1)² × 2(N + 1)² 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)² ×(N² −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)²×1 matrices and P is a (N² − 1) × 1 matrix.

Proposition 5.1 The (N² − 1) × (N² − 1) square matrix D^T M⁻¹D is invertible.

Proof. We just have to prove that the rank of the matrix D is N² − 1. Let us assume that the rank of the matrix D is strictly smaller than N² − 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^TM⁻¹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⁻¹D. Let us define the mapping φ by

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

Note that ψ₁(i)-1 and ψ₂(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) = (ψ₁(i), ψ₂(i)) and (r, s) = (ψ₁(j), ψ₂(j)).

Computing the elements a_i,j of the matrix D^TM⁻¹D yields

Thus, most of elements of the matrix D^T M⁻¹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⁻¹D and D^T F. Since the matrix D^T M⁻¹D is invertible, the equation (D^T M⁻¹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⁻¹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¹_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²_N (ξ_r ,ξ_s ). Let us define the matrix M ̃₁ = M ̃, where M ̃ is given by (5.7), the(N² −1)×(N² −1)matrix M ̃₂ =(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)² matrix D∗ by

We have to compute the matrix D^T∗M⁻¹∗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⁻¹ ∗D∗. Let us define the mapping φ∗ by

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

Note that ψ₁∗(i) and ψ₂∗(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) = (ψ₁∗(i), ψ₂∗(i)) and (r, s) = (ψ₁∗(j), ψ₂∗(j)).

Computing the elements a^∗_i,j of the matrix D^T∗ M⁻¹∗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⁻¹_*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⁻¹_*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⁻¹_*D_* . Note that this basis is orthonormal for the scalar product (.,.)_N. We can diagonalize the positive symmetric matrix D^T_*M⁻¹_*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)}² are in increasing order. Note that the matrix Λ is the matrix of f in a basis B′ = (h_i)_1≤i≤(N₊₁₎² of lP_N (Ω), which is orthonormal for the scalar product (., .)_N , the elements of which are the eigenvectors of the matrix D^T_*M⁻¹_*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^TP =R⁻¹P and the following equivalence

Since Z_N^∗ is a two dimensions space, B₀′ = (h₁ , h₂ ) 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⁻¹_*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¹(Ω). 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¹(Ω)² ×H²⁽Ω) 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⁻¹ for N = 8 and about 10⁻¹² 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⁻¹⁶ for N = 4 and about 10⁻¹² 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.