Numer. Math. 74: 69–84 (1996)

Numerische Mathematik

c Springer-Verlag 1996

Condition number and diagonal preconditioning: comparison of the p-version and the spectral element methods Jean-Franc¸ois Maitre1 , Olivier Pourquier2 1 Equipe d’Analyse Num´ erique Lyon-Saint Etienne, C.N.R.S. U.R.A. 0740, Ecole Centrale de Lyon, Dept. M.I.S., B.P. 163, F-69131 Ecully Cedex, France; e-mail: [email protected] 2 I.V.P. de Lorient, 10 rue Jean Zay, F-56325 Lorient Cedex, France; e-mail: [email protected]

Received November 24, 1994 / Revised version received March 20, 1995

Summary. In the framework of adaptive methods, bases of hierarchical type are used in the p-version of the finite element method. We study the matrices corresponding to the commonly used basis, introduced by Babuˇska and Szabo, in the case of d -dimensional rectangular elements for 2nd order elliptic boundary value problems. For the internal nodes, we show that the condition number is equivalent to p 4(d −1) and to p 4d for the stiffness and mass matrix, respectively. Moreover, we show that the usual diagonal preconditioning divides in the previous orders the exponents of p by two. Finally, we compare these results with those obtained for spectral elements (nodal basis). Mathematics Subject Classification (1991): 65F35, 65M60, 65N22

1. Introduction In the p-version of the finite element method, it is the degree p of the polynomials which is increased (with a constant mesh) to improve accuracy, in contrast to the standard h-version. For the p-version, we have the following approximation result (for the proof, we refer to Babuˇska, Suri (1990) [3] or Zhuang (1988) [11] and Bernardi, Maday (1992) for spectral elements [5]): ∀u ∈ H k (Ω) ,

∃up ∈ Vp ,

ku − up k1,Ω ≤ cp 1−k kukk ,Ω ,

where Ω is the domain of the boundary value problem and Vp is the approximation space (continuous functions which are polynomials of degree p on each element).

70

J.-F. Maitre, O. Pourquier

We consider a second order elliptic problem of homogeneous Dirichlet type in a domain Ω of Rd discretized with finite elements which are affine images of the reference element [−1, 1]d . Denoting by {Ni }1≤i ≤n the chosen basis of Vp , we have ∀u ∈ Vp ,

u=

n X

αi (u)Ni .

i =1

For d = 2, we write u = uI + uS + uV where uI , respectively uS , uV corresponds to the functions associated to the interiors, respectively sides and vertices. The part uI corresponds to the so called “bubble modes” which vanish on the boundary. Let A and M be the stiffness and mass matrix, respectively: Z Z ∇(Ni )∇(Nj ) and Mi, j = Ni Nj . Ai, j = Ω

Ω

AI (respectively MI ) is the diagonal block of A (respectively M ) corresponding to the bubble modes. Actually, for the solution of the linear system with the matrix A, essentially two types of methods have been studied and used: – internal Gaussian elminiation where the system is replaced by an upper block triangular system with AI as a diagonal block (see Babuˇska, Elman and Markley (1992) [2]), – block diagonal preconditioning with AI as a diagonal block (see Babuˇska, Craig, Mandel and Pitk¨aranta (1991) [1]). Both methods require the solution of a system with the block diagonal matrix AI . Therefore it is of interest to study the condition number of AI , more precisely its behavior as a function of the degree p. This study reduces to that of the matrix corresponding to the bubble modes on the reference element, since AI is block diagonal where each block corresponds to the bubble modes of one element. The condition number κ(. ) which we study is the usual one, that is the maximal eigenvalue divided by the minimal eigenvalue: κ(. ) = λmax (. )/λmin (. ). 2. Condition numbers of the mass and stiffness matrices obtained with p-version finite elements In all the following, we are interested in the basis introduced by Babuˇska and Szabo, which has been used by different authors (see Babuˇska, Elman, Markley (1992) [2], Babuˇska, Szabo (1991) [4], Oden (1991) [7], Pavarino (1994) [9], Zhuang (1988) [11]). The d -dimensional one is obtained from the 1-dimensional one by tensor products. The basis functions on Λ = [−1, 1] are given by the following expressions ({Lj }0≤j ≤p−1 denote the Legendre polynomials on Λ): 1 1 (1 − x ); ϕ1 (x ) = (1 + x ); 2 Z x 2 1 ϕj (x ) = Lj −1 (t)dt , kLj −1 kO,[−1,1] −1 ϕ0 (x ) =

∀2 ≤ j ≤ p .

Condition number in p-version of the finite element methods

71

For the bubble moddes (i ≥ 2), the stiffness matrix A1 is the identity. From {ϕi }2≤i ≤p , the basis of Qp0 (Λ), one obtains the basis {ϕi1 ⊗ . . . ⊗ ϕid }2≤i1 ,...,id ≤p of the space Qp0 (Λd ) of polynomials of degree p in each variable and vanishing on the boundary of Λd , where (ϕi1 ⊗ . . . ⊗ ϕid )(xi1 , . . . , xid ) = Qd k =1 ϕik (xik ). Md (respectively Ad ) is the matrix associatedR to this basis for the scalar R product L2 (Λd ) (respectively H01 (Λd )) defined by Λd (. )(. ) (respectively Λd ∇(. ). ∇(. )). d With a tensorial numbering, we define the canonical scalar product of R(p−1) Pp by (u, v)d = i1 ,...,id =2 ui1 ,...,id vi1 ,...,id . Denoting by λmax (. )(λmin (. )) the maximum (minimum) eigenvalue, we can state the following result: Proposition 1. The extreme eigenvalues of the d -dimensional mass and stiffness matrix can be expressed as functions of those of the 1-dimensional mass matrix by (1)

λmax (Md ) = λmax (M1 )d ,

λmin (Md ) = λmin (M1 )d ,

(2)

λmax (Ad ) = d λmax (M1 )d −1 ,

∀d ∈ N∗

λmin (Ad ) = d λmin (M1 )d −1 ,

∀d ∈ N∗

Proof. (1) follows immediately from the equalities: Z (Md )i1 ,...,id , j1 ,..., jd =

Λd

(ϕi1 ⊗ . . . ⊗ ϕid )(ϕj1 ⊗ . . . ⊗ ϕjd ) =

d Y

(M1 )ik , jk

k =1

Qd which means (Md θ, θ)d = θi1 ,...,id k =1 (M1 )ik jk θj1 ,..., jd (with the convention of summation on the repeated indices, which will be used in the following), and of the following lemma, which can be considered as generalizing one given in [8]. Lemma 1. Let (B k )1≤k ≤d be d symmetric definite positive p −1×p −1 matrices, then d Y

λmax (B k )(θ, θ)d

≥

θi1 ,...,id

k =1

d Y

(B k )ik , jk θj1 ,..., jd

k =1

≥

d Y

λmin (B k )(θ, θ)d ,

∀θ ∈ R(p−1)

d

k =1

Proof. We denote by I the identity matrix on Rp−1 . We use the d successive linear variable transformations by (B i )1/2 , i = 1 to d , in order to simplify the central term of Lemma 1, and note θ(k ) the tensor obtained after k such transformations. With θ(1).j2 ,..., jd = (B 1 )1/2 θ.j2 ,..., jd , ∀2 ≤ j2 , . . . , jd ≤ p , we obtain:

72

J.-F. Maitre, O. Pourquier

θi1 ,...,id

d Y

(B k )ik , jk θj1 ,..., jd

k =1 −1/2

= (B 1 )i1 ,r θ(1)r,i2 ,...,id

d Y

−1/2

(B k )ik , jk (B 1 )j1 ,s θ(1)s, j2 ,..., jd .

k =1

From −1/2

−1/2

(B 1 )i1 ,r (B 1 )i1 , j1 (B 1 )j1 ,s

−1/2

−1/2

1/2

1/2

= (B 1 )i1 ,r (B 1 )i1 ,s = (B 1 )i1 ,r (B 1 )s,i1 = Is,r ,

we get θi1 ,...,id

d Y

(B k )ik ,jK θj1 ,..., jd = θ(1)i1 ,i2 ,...,id

k =1

d Y

(B k )ik ,jk θ(1)i1 , j2 ,..., jd0

k =2

since θ(1)r,i2 ,...,id Is,r θ(1)s, j2 ,...,jd = θ(1)i1 ,i2 ,...,id θ(1)i1 , j2 ,..., jd . For θ(2), we have θ(2)i1 ,., j3 ,..., jd = (B 2 )1/2 θ(1)i1 ,.,j3 ,...jd ,

for all 2 ≤ j2 , . . . , jd ≤ p ,

which gives: θi1 ,...,id

d Y

(B k )ik ,jk θj1 ,..., jd = θ(2)i1 ,i2 ,i3 ,...,id

k =1

d Y

(B k )ik , jk θ(2)i1 ,i2 , j3 ,..., jd .

k =3

We obtain finally θi1 ,...id

d Y

(B k )ik ,jk θj1 ,..., jd = (θ(d ), θ(d ))d .

k =1

Now, successively we return to (θ, θ)d : (θ(d ), θ(d ))d

(θ(d ), θ(d ))d

=

θ(d )i1 ,...,id θ(d )i1 ,...,id

=

(B d )id ,r θ(d − 1)i1 ,...,id −1 ,r (B d )id ,s θ(d − 1)i1 ,...,id −1 ,s

=

(B d )r,s θ(d − 1)1i ,...,id −1 ,r θ(d − 1)i1 ,...,id −1 ,s   (B d )θ(d − 1)i1 ,...,id −10. , θ(d − 1)i1 ,...id −10.

=

1/2

1/2

1

so that λmin (B d )(θ(d − 1), θ(d − 1))d

≤

(θ(d ), θ(d ))d

≤

λmax (B d )(θ(d − 1), θ(d − 1))d .

The final result follows by returning successively from d − 1 to 0, that is to tensor θ. u t

Condition number in p-version of the finite element methods

73

In fact, this lemma implies: λmax (Md ) ≤ (λmax (M1 ))d

and λmin (Md ) ≥ (λmin (M1 ))d .

To prove that the equalities are attained, one considers the special tensor θ˜ defined by d Y ˜θi1 ,...,id = vik , 2 ≤ i1 , . . . , id ≤ p , k =1

where v is an eigenvector of M1 associated to λmax (M1 ) or λmin (M1 ). For (2), we have the equalities: Z

∇ ϕi1 ⊗ . . . ⊗ ϕid . ∇ ϕj1 ⊗ . . . ⊗ ϕjd (Ad )i1 ,...,id , j1 ,...jd = Λd

=

d X

δik jk

k =1

d Y

(M1 )il jl

l =1,l/ =k

which means (Ad θ, θ)d

=

θi1 ,i2 ,...,id

d Y

(M1 )il jl θi1 ,j2 ,..., jd

l =2

+θi1 ,...,ik ,...,id

d Y

(M1 )il jl θj1 ,...,ik ,..., jd + . . .

l =1,l/ =k

+θi1 ,...,id −1 ,id

dY −1

(M1 )il jl θj1 ,..., jd −1 ,id .

l =1

Applying Lemma 1 to each of the d terms, we obtain the inequalities: λmax (Ad ) ≤ d (λmax (M1 ))d −1

and λmin (Ad ) ≥ d (λmin (M1 ))d −1 .

To prove that the equalities are attained, one considers the same particular θ˜ as above. u t Proposition 1 shows that to obtain estimates of the extreme eigenvalues of the mass and stiffness matrices for dimension d , one only has to get estimates for the extreme eigenvalues of the mass matrix in dimension 1. In the sequel we will use the notation f ≈ g, and say that f is equivalent to g, if there exist two positive constants C and C 0 independent of g, such that C g ≤ f ≤ C 0 g. Proposition 2. Let M1 be the mass matrix corresponding to the bubble modes on the reference element, then (3)

λmin (M1 ) ≈ p −4

(4)

λmax (M1 ) ≈ 1 .

74

J.-F. Maitre, O. Pourquier

Proof. 1st step. (lower bound of λmin , upper bound of λmax ). We recall the inverse inequality (see Bernardi, Maday (1992) [5]): ∀v ∈ Qp0 ([−1, 1]) ,

kvk20,[−1,1] ≥ Cp −4 |v|21,[−1,1] ,

which, with our notation, can be written as follows: (M1 α, α)1 ≥ Cp −4 (A1 α, α)1 ∀α ∈ Rp−1 ,

with v =

p X

α i ϕi .

i =2

Since A1 is the identity, i.e. (A1 α, α)1 = (α, α)1 , we obtain λmin (M1 ) ≥ Cp −4 . For the upper bound of λmax (M1 ), we can use the explicit knowledge of M1 , the only nonzero elements of which are:  2  , 2≤i ≤p (M1 )i ,i =    (2i + 1)(2i − 3)  −1 √ , 2≤i ≤p−2, (M1 )i ,i +2 =   (2i + 1) (2i − 1)(2i + 3)    (M1 )i ,i −2 = (M1 )i −2,i , 4 ≤ i ≤ p We obtain λmax (M1 ) ≤

2 5

+

√1 5 21

by Gerschgorin’s theorem.

Remark. We could have ignored the detailed knowledge of M1 and used Poincar´e’s inequality. 2nd step. (lower bound of λmax , upper bound of λmin ). For the upper bound of λmin (M1 ), we use the value of the Rayleigh quotient for a well chosen function (very oscillating). More precisely, we have the following result: Lemma 2. The function ψp = L00p+2 −

(p+3)(p+4) 0 Lp+1 4

2

of Qp0 (Λ) verifies: 2

40 {|ψp |1,[−1,1] } ∼ p 4 {kψp k0,[−1,1] } where the notation f (p) ∼ g(p) means

f (p) g(p)

→ 1 when p → +∞.

Proof. We note f (n) for the n th derivative of f . We need the following properties of the Legendre polynomials: Q l ,|l |
) 2 (c) kL(k m k0,Λ =

k m 4k −2 X (−1)l +1 + O(m 4k −3 ) for m ≥ 2k − 1 . 22k −2 (k − l )!(k + l − 1)! l =1

Condition number in p-version of the finite element methods

Z Λ

) (k −1) L(k m Lm−1

75

k −1 m 4k −4 X (−1)l +1 + O(m 4k −5 ) 22k −3 (k − l − 1)!(k + l − 1)!

=

l =1

for m ≥ max(2k − 2, k + 2).

(d)

We obtain ψp (1) = 0 and ψp (−1) = 0 with (a) and (b) respectively, so that ψp ∈ Qp0 (Λ). We have 2  (p + 3)(p + 4) 2 00 2 kL0p+1 k20,Λ kψp k0,Λ = kLp+2 k0,Λ + 4 Z (p + 3)(p + 4) − L00p+2 L0p+1 2 Λ and with (c) and (d):   1 1 p6 p6 1 p6 2 − + − + O(p 5 ) , kψp k0,Λ = 2 6 4 16 2 4 We have |ψp |21,Λ

that is kψp k20,Λ =

p6 + O(p 5 ) . 48

2 (p + 3)(p + 4) kL00p+1 k20,Λ 4 Z (p + 3)(p + 4) 00 − L000 p+2 Lp+1 2 Λ 2 kL000 p+2 k0,Λ +

=



and with (c) and (d):  10      1 1 p 1 p 10 1 1 p 10 1 1 − + + − − − + O(p 9 ) , |ψp |21,Λ = 12 24 120 16 2 6 64 6 24 16 p 10 9 1920 + O(p ). 2 40|ψp | → have p 4 kψp k1,Λ 2 0,Λ

so that |ψp |21,Λ = Thus we

1 when p → +∞.

t u

Remark. A similar function has already been used by Bernardi, Maday (1992) [5]. Lemma 2 gives a satisfactory bound for λmin (M1 ): with α˜ associated to ψp , we have |ψp |21,[−1,1] = α˜ T A1 α˜ = (α, ˜ α) ˜ 1 , and the existence of α˜ ∈ Rp−1 such that 1 −4 ˜ α) ˜ 1 ∼ 40 p (α, ˜ α) ˜ 1 , which implies the asymptotic upper bound p −4 /40. (M1 , α, t u From Propositions 1, 2, we deduce the following final result: Theorem 1. For Ad (respectively Md ) the stiffness (respectively mass) matrix corresponding to the bubble modes of the Babuˇska-Szabo basis on the reference element in d dimension, we have asymptotically in p, the degree of the polynomials ∀d ∈ N∗

(5)

λmax (Md ) ≈ 1 and λmin (Md ) ≈ p −4d ,

(6)

λmax (Ad ) ≈ 1 and λmin (Ad ) ≈ p −4(d −1) , κ(Md ) ≈ p

4d

and κ(Ad ) ≈ p

4(d −1)

,

∀d ∈ N∗

∀d ∈ N∗ .

76

J.-F. Maitre, O. Pourquier

3. Diagonal preconditioning of the mass matrix Denoting by ∆d the diagonal of the mass matrix Md , we have the following result: Proposition 3. The extreme eigenvalues of the d -dimensional mass matrix preconditioned by its diagonal can be expressed as functions of those of the same matrix in one dimension by λmax (∆−1 d Md )

(7)

λmin (∆−1 d (Md )

(8)

= =

d λmax (∆−1 1 M1 ) , d λmin (∆−1 1 M1 ) ,

∀d ∈ N∗

∀d ∈ N∗

−1/2 −1/2 Proof. The scaled matrix ∆d Md ∆d corresponds to the tensor M˜ d built from the basis {ϕ˜ i }2≤i ≤p with ϕ˜ i = ϕi /kϕi kL2 (Λ) . The proof of Proposition 1, which is independent of the basis, allows us to ˜ d ) = λmax (M ˜ 1 )d and λmin (M ˜ d ) = λmin (M˜ 1 )d , ∀d ∈ N∗ . show that: λmax (M Formulas (7) and (8) follow directly since M˜ d and ∆−1 d Md have the same eigenvalues. u t

Estimates of the extreme eigenvalues of ∆−1 1 M1 are given by the following result: Proposition 4. The extreme eigenvalues of the one-dimensional scaled mass matrix ∆−1 1 M1 verify −2 λmin (∆−1 1 M1 ) ≈ p

(9)

λmax (∆−1 1 M1 ) ≈ 1

(10)

Proof. The knowledge of M1 and ∆1 permits us to write down the elements of −1/2 −1/2 M˜ = ∆1 M1 ∆1 , explicitly:  (M˜ )i ,i = 1 fors2 ≤ i ≤ p ,        −1 2i + 5 2i − 3 ˜ for 2 ≤ i ≤ p − 2 , (M )i ,i +2 =  2 2i − 1 2i + 3    ˜ )i −2,i for 4 ≤ i ≤ p − 2 ˜ )i ,i −2 = (M (M By Gerschgorin’s theorem one has   p X  ˜ )i, j |, 4 ≤ i ≤ p − 2 ≤ 1 − [6/(2p − 1)(2p − 9)] max |(M   j =2,j / =i

and then

λmax (∆−1 1 M1 ) ≤ 2 .

The scaled matrix M˜ is diagonal dominant, while M1 is not, since 3 1 ˜ )i ,i +2 | ≤ 1 − < . |(M 2 (2i − 1)(2i + 3) 2

Condition number in p-version of the finite element methods

77

Thus we can use again Gerschgorin’s theorem to obtain: 2 λmin (∆−1 1 M1 ) ≥ 6/[(2p − 1)(2p − 9)] ≥ 3/[2p ] .

˜ ) ≥ max{(M ˜ )i ,i , 2 ≤ i ≤ p} and (M˜ )i ,i = 1. (10) results from λmax (M −1 (9) results from the following upper bound of λmin (∆−1 1 M1 ): λmin (∆1 M1 ) ≤ 2 λmin (M1 )/λmin (∆1 ), with λmin (∆1 ) = (M1 )p,p = (2p+1)(2p−3) and from equality (3). t u Combining results (7), (8) and (9), (10), we can state the final result: Theorem 2. For the scaled mass matrix ∆−1 d Md corresponding to the bubble modes of the Babuˇska-Szabo basis on the reference element in d dimension, we have asymptotically in p, the degree of the polynomials λmax (∆−1 d Md ) ≈ 1

(11)

−2d λmin (∆−1 d Md ) ≈ p

(12)

2d , κ(∆−1 d Md ) ≈ p

∀d ∈ N∗ .

4. Diagonal preconditioning of the stiffness matrix Denoting by Dd the diagonal of the stiffness matrix Ad , we have the following result: Proposition 5. The extreme eigenvalues of the stiffness matrix preconditioned by its diagonal can be expressed as functions of those of the mass matrix in one dimension: d −1 , λmax (Dd−1 Ad ) = λmax (∆−1 1 M1 )

(13)

λmin (Dd−1 Ad )

(14)

=

d −1 λmin (∆−1 1 M1 )

,

∀d ∈ N ∗ ∀d ∈ N ∗

Proof. From Z (Ad )i1 ,...,id , j1 ,..., jd

= =

Λd d X k =1



∇ ϕi1 ⊗ . . . ⊗ ϕid . ∇ ϕj1 ⊗ . . . ⊗ ϕjd δik jk

d Y

(M1 )il jl ,

l =1,l/ =k

we see that (Ad θ, θ)d is made up of d terms. Each term is of the type (Md −1 β, β)d −1 where β is for example βi2 ,...,id = θi1 ,i2 ,...,id , and can be bounded as follows: (Md −1 β, β)d −1 ≤ λmax (∆−1 d −1 Md −1 )(∆d −1 β, β)d −1 . Summing the d terms leads to the inequality:

78

J.-F. Maitre, O. Pourquier

(Ad θ, θ)d ≤ λmax (∆−1 d −1 Md −1 )

d X

d Y

θi1 ,...,id

k =1

(M1 )il il θi1 ,...,id .

l =1,l/ =k

Identifying the tensor corresponding to Dd in the above expression, we obtain: (Ad θ, θ)d ≤ λmax (∆−1 d −1 Md −1 )(Dd θ, θ)d . In the same way, we can prove (Ad θ, θ)d ≥ λmin (∆−1 d −1 Md −1 )(Dd θ, θ)d . The preceding inequalities provide an upper and a lower bound for λmax −1 (∆−1 1 M1 ) and λmin (∆1 M1 ), respectively. To prove that these bounds are attained, i.e. the equalities (13), (14), we exhibit a special θ˜ defined by θ˜i1 ,...,id =

d Y

vik /{(m)ik }1/2

k =1 −1/2

−1/2

where v is an eigenvector of the matrix ∆1 M1 ∆1 associated with λmax −1 (∆−1 M ) or λ (∆ M ), and (m) = (M ) . 1 min 1 i 1 i ,i k k k 1 1 We obtain now by using Proposition 3:   d h d d d h Y i X i Y Y ˜ θ) ˜d = vik /{(m)ik }1/2 vjq /{(m)jq }1/2 δ in jn (M1 )il jl (Ad θ,   n=1

k =1

 ˜ θ) ˜ d = λmax (∆−1 M1 ) d −1 (Ad θ, 1

d X



q=1

l =1,l/ =n

vi2n /min



n=1

d Y

v ik v ik

k =1,k/ =n

˜ θ) ˜ d = {λmax (∆−1 M1 )}d −1 (Dd θ, ˜ θ) ˜ d. Thus we have: (Ad θ, 1 t The same proof applied to λmin proves (13), (14). u Combining them with the estimates (9), (10) for d = 1, we obtain the final result: Theorem 3. For the scaled stiffness matrix Dd−1 Ad corresponding to the bubble modes of the Babuˇska-Szabo basis on the reference element in dimension d , we have (15) (16)

λmax (Dd−1 Ad ) ≈ 1

λmin (Dd−1 Ad ) ≈ p −2(d −1)

κ(Dd−1 Ad ) ≈ p 2(d −1) ,

∀d ∈ N∗ .

5. Remark on the condition number of the global stiffness matrix in 2-D The following result shows that, asymptotically in p, the condition number of C , the block-diagonal preconditioner, is of the same order of magnitude as that of AI :

Condition number in p-version of the finite element methods

79

Lemma 3. For the stiffness matrix A on Ω, and the block-diagonal preconditioning C containing AI as a block, we have κ(AI ) ≤ κ(C ) ≤ cκ(AI ) where c is independent of the degree p. Proof. For α = (αV , αS,1 , . . . , αS,Nc , αI,1 , . . . , αI,Ni ), where the mesh on Ω contains Nc element sides and Ni element interiors, we have: αT C α = (αV )T AV αV +

Nc X

(αS,k )T AS,k αS,k +

k =1

Ni X

(αI,k )T AI,k αI,k

k =1

where AV , AS,1 , . . . AS,Nc , AI,1 , . . . , AI,Ni are the diagonal blocks of C (AI contained the blocks AI,1 , . . . , AI,Ni ). First of all AV , corresponding to the vertices of the mesh, is independent of the degree p. With the affine mapping, each block AI,k has the same behavior in p as A2 , the 2-dimensional stiffness matrix built with the shape functions ϕi ⊗ ϕj for all 2 ≤ i , j ≤ p associated to the internal nodes of the reference element. We have: λmax (A2 ) ≈ 1 and λmin (A2 ) ≈ p −4 (see Sect. 2). With the affine mapping, each block ASk has the same behavior in p as Aˆ S , the matrix built with the shape functions ϕj ⊗ ϕ0 for all 2 ≤ j ≤ p associated with a side of the reference element. For Aˆ S , we prove below that the behavior of its extremal eigenvalues is independent of the degree p. For 2 ≤ i , j ≤ p, we have Z Z Z Z Z

0 0 ˆ ∇ N(1)i ∇ N(1)j = ϕi ϕj ϕ0 ϕ0 + ϕi ϕj ϕ00 ϕ00 , (AS )ij = Λ2

Λ

Λ

Λ

Λ

Z

so that (Aˆ S )ij = 2/3δij + (1/2)

Λ

ϕi ϕj .

R / {j − For Aˆ S , only three diagonals are non zero since Λ ϕi ϕj = 0 if i ∈ 2, j , j + 2}.R From | Λ ϕi ϕj | ≤ 25 (see Sect. 2), we obtain by using Gerschgorin’ s theorem   p X  |(Aˆ S )ij | , 2 ≤ i ≤ p λmax (Aˆ S ) ≤ max   j =2

and λmax (Aˆ S ) ≤ max

 p X 2 

j =2

 

  p Z X  1 δij , 2 ≤ i ≤ p + max ϕ ϕ , 2 ≤ i ≤ p i j 2    3 Λ j =2

We conclude that λmax (Aˆ S ) ≤ 23 + 12 3 25 , so that λmax (Aˆ S ) ≤ Naturally λmax (Aˆ S ) ≥ max((Aˆ S )kk ) = (Aˆ S )22 = 13 15 .

19 15 .

80

J.-F. Maitre, O. Pourquier

By Gerschgorin’s theorem, we have now,   p   X |(Aˆ S )ij |, 2 ≤ i ≤ p λmin (Aˆ S ) ≥ min |(Aˆ S )ii | −   j =2, j / =i   p  X   ≥ min |(Aˆ S )ii | − max |(Aˆ S )ij |, 2 ≤ i ≤ p   j =2, j / =i   Z   Z p   X 1 2 1 , 2 ≤ i ≤ p ≥ min + ϕ2i − max ϕ ϕ i j 2   3 2 Λ Λ j =2, j / =i

4 so that λmin (Aˆ S ) ≥ 23 − 12 2 25 and λmin (Aˆ S ) ≥ 15 . 13 ˆ ˆ Naturally λmin (AS ) ≤ (AS )22 = 15 . Because the greatest (respectively smallest) eigenvalue of C equals the maximum (respectively minimum) of the greatest (respectively smallest) eigenvalues of the blocks of C , Lemma 3 is proved since λmax (AV ), λmax (AS,k ), λmax (AI,k ), λmin (AV ), λmin (AS,k ) are equivalent to 1, i.e. independent of p, and λmin (AI,k ) is t equivalent to p −4 . u

Corollary 1. For the 2-D case (only), the condition number κ(A) of the stiffness matrix for the discrete homogeneous Dirchlet problem corresponding to a uniform square-mesh of Ω satisfies αp 4 ≤ κ(A) ≤ βp 4 (1 + [ln(p)]2 ) where α, β are constants independent of p. Proof. In the proof of Lemma 3, we have shown that the asymptotic behaviour in p of the extremal eigenvalues of C are: λmin (C ) ≈ p −4 and λmax (C ) ≈ 1. Using the result κ(C −1 A) ≤ γ(1+[ln(p)]2 ) proved by Babuˇska, Craig, Mandel and Pitk¨aranta (1991) [1], we obtain the upper bound since κ(A) ≤ κ(C )κ(C −1 A). The lower bound results from the fact that C is built from diagonal blocks of A. t u Corollary 2. For the 2-D case (only), the condition number κ(D −1 A) of the stiffness matrix with diagonal preconditioning for the discrete homogeneous Dirichlet problem corresponding to a uniform square-mesh of Ω verifies γp 2 ≤ κ(D −1 A) ≤ δp 2 (1 + [ln(p)]2 ] ; where γ, δ are constants independent of p. Proof. With a similar proof as that of Lemma 3, we obtain κ(D −1 AI ) ≤ κ(D −1 C ) ≤ cκ(D −1 AI ). Moreover, since κ((D −1 C )−1 D −1 A) = κ(C −1 A), we can use the results of Theorem 3. u t

Condition number in p-version of the finite element methods

81

6. Comparison with the spectral element matrices For the spectral element method, we study only the condition number for the matrix built with bubble functions on the reference element Λd = [−1, 1]d , and consider the nodal basis associated with the Gauss-Lobatto quadrature points. In 1-D, this gives p + 1 points, and polynomials of degree 2p − 1 are integrated exactly. Therefore the stiffness matrix As,1 is computed exactly and the mass matrix only approximately. For d > 1, the stiffness matrix As,d , corresponding R to the scalar product Λd ∇(. )∇(. ), is not computed exactly by the Gauss-Lobatto quadrature formula. The computed mass matrix Ds,d is diagonal since the shape functions are the Lagrange polynomials associated to the Gauss-Lobatto points (see Bernardi, Maday (1992) [5]). Denoting by (αij ), 2 ≤ i , j ≤ p, the coefficients of As,1 and by ρk , 2 ≤ k ≤ p, the coefficients of Ds,1 , we have: (As,d )i1 ,...,id , j1 ,..., jd =

d X

d Y

αik jk

k =1

(Ds,1 )il jl =

d X k =1

l =1,l/ =k

d Y

αik jk

ρi l δ i l j l .

l =1,l/ =k

It can be shown that the ρk satisfy (see Bernardi, Maday (1992) [5]): cp −2 ≤ ρk ≤ c 0 p −1 ,

∀2 ≤ k ≤ p .

Lemma 4. For the stiffness matrix of spectral element method, we have c1 p 1−d λmax (As,1 ) ≤ λmax (As,d ) and λmin (As,d ) ≤ c2 p 1−d λmin (As,1 ) . Proof. Since (As,d θ, θ)d = θi1 ,...,id

d X

αik jk

k =1

 d  Y 

ρil δil jl

l =1,l/ =k

  

θj1 ,..., jd ,

each of the d terms is identical to the first term up to a permutation. Consider one term: θi1 ,i2 ,...,id αi1 j1 θj1 ,i2 ,...,id =

p X

As,1 θ.,i2 ,...,id , θ.,i2 ,...,id

i2 ,...,id =2


1

.

Since the weights ρi are positive, we obtain: (As,d θ, θ)d ≥ θi1 ,i2 ,...,id αi1 j1 θj1 ,i2 ,...,id

d Y

ρi k .

k =2

For λmax (As,d ), we consider θ˜i1 ,...,id = vi1 , ∀2 ≤ i1 , . . . , id ≤ p, where v is an eigenvector associated to λmax (As,1 ), and obtain !d −1 p X ˜ θ) ˜ d ≥ (As,1 v, v)1 ρi . (As,d θ, i =2

82

J.-F. Maitre, O. Pourquier

Pp 4 ˜ θ) ˜d ≥ Since i =2 ρi = 2− p(p+1) (see Bernardi, Maday (1992) [5]), we have (As,d θ, d −1 ˜ θ) ˜ d = (p −1) (v, v)1 , we find that λmax (As,d ) ≥ k λmax (As,1 )(v, v)1 , and since (θ, C1 p 1−d λmax (As,1 ), which proves the first result. For λmin (As,d ), with βk = θk ,i2 ,...,id , we have (As,1 β, β)1

d Y

ρii ≤ k2 p 1−d (As,1 β, β)1

since

ρk ≤ c 0 p −1 , ∀2 ≤ k ≤ p .

l =2

The same inequality is satisfied by each term and we obtain for all θ: (As,d θ, θ)d ≤ k2 p 1−d θi1 ,...,id

d X k =1

αik jk

d Y

δil jl θj1 ,..., jd .

l =1,l/ =k

Qd For θ˜i1 ,...,id = k =1 wik where w is an eigenvector of As,1 associated with λmax (As,1 ), we have ˜ θ) ˜ ≤ k dp 1−d λ (A ) {(w, w) }d (As,dθ, d 2 min s,1 1 ˜ θ) ˜ . which implies the second result since {(w, w)1 }d = (θ, d

t u

Knowing that κ(As,1 ) ≈ p 3 (see Bernardi, Maday (1992) [5]), we can give the following result: Theorem 4. The stiffness matrix As,d for the spectral element method in Rd satisfies cp 3 ≤ κ(As,d ) , ∀d ∈ N . Remark 1. For d = 1, we know (see above) that κ(As,1 ) ≤ c 0 p 3 . For d = 2, 3, we have also κ(As,d ) ≤ c 0 p 3 (see above reference and Azaiev, Dauge and Maday (1994) [0]). Remark 2. Concerning the diagonal preconditioning, we emphasize that there is actually no theoretical result (for the stiffness matrix since the mass matrix is diagonal).

7. Numerical results We now present numerical results for the condition number with and without diagonal preconditioning, for the stiffness matrices of the p-version and the spectral elements (Fig. 1 for d = 2 and Fig. 2 for d = 3). They permit us to verify the results for which there is a proof and to provide conjectures for others. In the figures, continuous lines are for the p-version and dotted lines for the spectral elements. In each case, the lowest curve corresponds to the diagonal preconditioning. The figures represent κ(. ) as a function of p in logarithmic scales.

Condition number in p-version of the finite element methods

83

Fig. 1.

Fig. 2.

Comments on the figures. For computing reasons, our results are limited to p = 30 in 2D and to p = 9 in 3D. For the p-version, the behavior of κ(Ad ) (respectively κ(Dd−1 Ad )) as a function of p are approximately in accord with our theoretical results p 4(d −1) (respectively p 2(d −1) ). While diagonal preconditioning is very efficient for the p-version, this is not the case for the spectral element: in 2D, the behavior in p 3 seems to become in p 2 , only for very large p, and in 3D there is almost no effect of the diagonal preconditioning.

References 0. 1.

Azaiev, M., Dauge, M., Maday, Y. (1994): M´ethodes spectrales et des e´ l´ements spectraux. Publication 94-17 de l’I.R.M.A.R. et cours I.N.R.I.A., Rennes, France Babuˇska, I., Craig, A., Mandel, J., Pitk¨aranta, J. (1991): Efficient preconditioning for the pversion F.E.M. in two dimension. S.I.A.M. J. Numer. Anal. 28(3), 624–661

84

J.-F. Maitre, O. Pourquier

2.

Babuˇska, I., Elman, H.C., Markley, K. (1992): Parallel implementation of the h-p version of the finite element method on a shared-memory architecture. S.I.A.M. J. Sci. Stat. Comput. 13(6), 1433–1459 3. Babuˇska, I., Suri, M. (1990): The h-p version of the F.E.M., an overview. Comput. Methods Appl. Mech. Eng. 80, 5–26 4. Babuˇska, I., Szabo, B.A. (1991): Finite element analysis. Wiley Interscience, New York 5. Bernardi, C., Maday, Y. (1992): Approximations spectrales de probl`emes aux limites elliptiques. Math. & Appl., Springer Verlag 6. Maˆıtre, J.F., Pourquier, O. (1994): Conditionnements et pr´econditionnements diagonaux des syst`emes pour la p-version des m´ethodes d’´el´ements finis pour des probl`emes elliptiques d’ordre 2. C.R. Acad. Sci. Paris 318 S´erie I (6), 583–586 7. Oden, J.T. (1991): Theory and implementation of high-order adaptive h-p methods for the analysis of incompressible viscous flows. Computational non linear mechanics in aerospace eng. (AIAA progress in aeronautics and astronautics, S.N. Atluri) 8. Olsen, E.T., Douglas, J.: Condition numbers of matrices arising in the p-version of the finite element method. Numer. Math. (to appear) 9. Pavarino, L.F. (1994): Adaptive Schwarz methods for the p-version finite element method. Numer. Math. 66(4), 493–516 10. Pourquier, O. (1994): M´ethode des e´ l´ements finis de haut degr´e (p-version): estimation du conditionnement des matrices et construction de pr´econditionneurs. Thesis 94-35 (Ecole Centrale de Lyon) 11. Zhuang, G.W. (1988): Hierarchical elements, local mappings and the h-p version of the F.E.M. (I) and (II). J. Comput. Math. 6(1,2), 54–67, 142–155

This article was processed by the author using the LaTEX style file pljour1m from Springer-Verlag.