Mixed discontinuous Galerkin approximation of the elasticity eigenproblem
Felipe Lepe Salim Meddahi David Mora Rodolfo Rodríguez
math.NA 65N30, 65N12, 65N15, 74B10

MengeROS: a Crowd Simulation Tool for Autonomous Robot Navigation
    Anoop Aroor1, Susan L. Epstein1,2 and Raj Korpan1

Department of Computer Science

The Graduate Center1 and Hunter College2 of The City University of New York

New York, NY 10065

aaroor@gradcenter.cuny.edu, susan.epstein@hunter.cuny.edu, rkorpan@gradcenter.cuny.edu

    December 30, 2023
===============================================================================================================================================================================================================================================================================




























We introduce a discontinuous Galerkin method for the mixed formulation of the elasticity eigenproblem with reduced symmetry. The analysis of the resulting discrete eigenproblem does not fit in the standard spectral approximation framework since the underlying source operator is not compact and the scheme is nonconforming. We show that the proposed scheme provides a correct approximation of the spectrum and prove  asymptotic error estimates for the eigenvalues and the eigenfunctions. Finally, we provide several  numerical tests to illustrate the performance of the method and confirm the theoretical results.



Keywords:
Mixed elasticity equations, spectral problems,  finite elements,  discontinuous Galerkin methods, 
error estimates.


AMS Subject Classification: 65N30, 65N12, 65N15,  74B10




myheadings
plain
S. MEDDAHIA  discontinuous Galerkin method for an  
elasticity  eigenproblem




§ INTRODUCTION



We present a discontinuous Galerkin (DG) approximation of the linearized vibrations of an elastic structure. In many applications, the displacement field is not necessarily the variable of  primary interest. We consider here the dual-mixed formulation of the elasticity eigenproblem because it delivers a direct finite element approximation of the Cauchy stress tensor and it permits to deal safely with nearly incompressible materials.
 
A mixed finite element approximation of the eigenvalue elasticity problem with reduced symmetry has been analyzed in <cit.>. It consists in a formulation that only maintains  the stress tensor as primary unknown, besides the rotation whose role is the weak imposition of the symmetry restriction. It is shown that a discretization based on the lowest order Arnorld-Falk-Winther element provides  a correct spectral approximation and quasi optimal asymptotic error estimates for the eigenvalues and the eigenfunctions. 
 
 
The ability of DG methods handle  efficiently hp-adaptive strategies make them suitable for the numerical simulation of physical systems related to elastodynamics.  Our aim here is to introduce an interior penalty discontinuous Galerkin version for the H(div)-conforming finite element space employed in <cit.>. The k^th-order of this method amounts to approximate the Cauchy stress tensor and the rotation by discontinuous finite element spaces of degree k and k-1 respectively. We point out that an H(curl)-based interior penalty discontinuous Galerkin method has also been introduced in <cit.> for the Maxwell eigensystem. The DG approximation we are considering here may be regarded as its counterpart in the H(div)-setting. As in <cit.>, our analysis requires conforming meshes, but the DG method still permits one to employ different polynomial element  orders in the same triangulation. A further advantage of this DG scheme is that it allows to implement high-order elements in a mixed formulation by using standard shape functions. Let us remark that the DG method has also been analyzed in <cit.> for the Laplace operator.
 
 It is well known that the underlying source operator corresponding to mixed formulations is generally not compact. In our case, this operator admits a non physical zero eigenvalue  whose eigenspace is  infinite dimensional. It is then essential to use a scheme that is safe from the pollution that may appear in the form of spurious eigenvalues interspersed among the physically relevant ones.  It turns out (cf. <cit.>) that, for mixed eigenvalue problems, the conditions guarantying the convergence of the source problem does not necessarily a correct spectral approximation (as it happens for compact operators <cit.>).

It has been shown in <cit.> that DG methods can also benefit from the general theory developed in <cit.> to deal with the spectral numerical analysis of non-compact operators. We follow here the same strategy, combined with techniques from <cit.>, to prove that our numerical scheme is spurious free. We also establish asymptotic error estimates for the eigenvalues and eigenfunctions. We treat with special care the analysis of the limit problem obtained when the Lamé coefficient tends to infinity.

We end this section with some of the notations that we will  use below. Given
any Hilbert space V, let V^n and V^ denote, respectively,
the space of vectors and tensors of order n (n= 2, 3)  with
entries in V. In particular,  is the identity matrix of
^ and 0 denotes a generic null vector or tensor. 
Given :=(τ_ij) and :=(σ_ij)∈^, 
we define as usual the transpose tensor ^:=(τ_ji), 
the trace :=∑_i=1^nτ_ii, the deviatoric tensor 
^:=-1/n(), and the
tensor inner product ::=∑_i,j=1^nτ_ijσ_ij. 

Let Ø be a polyhedral Lipschitz bounded domain of ^n with
boundary . For s≥ 0, ·_s,Ø stands indistinctly
for the norm of the Hilbertian Sobolev spaces , ^n or
^, with the convention ^̋0(Ø):=. We also define for
s≥ 0 the Hilbert space 
:=∈^: ∈^n, whose norm
is given by ^2_
:=_s,Ø^2+^2_s,Ø and denote
:=^̋0(𝐝𝐢𝐯,Ø).

Henceforth, we denote by C generic constants independent of the discretization
parameter, which may take different values at different places.




§ THE MODEL PROBLEM



Let Ø⊂^n (n=2,3) be an open bounded Lipschitz 
polygon/polyhedron representing an elastic body. We denote by  the outward unit normal 
vector to   and assume that  =Γ_D∪Γ_N, with int(Γ_D)∩int(Γ_N) = ∅.  
The solid is supposed to be isotropic
and linearly elastic with  mass density ρ and Lamé constants μ
and λ. We assume that the structure is fixed at Γ_D≠∅ and free
of stress on Γ_N. 
We can combine the constitutive law 

    ^-1=() Ø,

and the equilibrium equation 

    ω^2   =  ρ^-1Ø,

to eliminate either the displacement field  or the Cauchy stress tensor  from the global spectral 
formulation of the elasticity problem. Here,  
():=1/2[∇+(∇)^] is  the linearized strain tensor,  and 
: ^→^ is the Hooke operator, which is given in terms of the Lamé coefficients 
λ and μ by

    :=λ() + 2μ    ∀ ∈^.

Opting for the elimination of the displacement  and  maintaining the stress tensor 
 as a main variable leads to the following dual mixed formulation 
of the elasticity eigenproblem:  
Find :Ø→^ symmetric,  :Ø→^ skew symmetric 
and ω∈ℝ such that,

    [      -∇(ρ^-1)             = ω^2 ( ^-1 + )            Ø,;                           =                        Γ_D,;                           =                        Γ_N. ]

We notice that the additional variable :=1/2[∇-(∇)^] is the rotation. 
It acts as a Lagrange multiplier for the symmetry restriction. We also point out that 
the displacement can be recovered and also post-processed at the discrete level by using identity (<ref>).


Taking into account that the Neumann boundary condition becomes essential in the mixed formulation,  we consider the closed subspace  of  given by

    :={∈:  = on Γ_N}.

 The rotation  will be sought in the space

    :=∈^: ^=-.



We introduce the symmetric bilinear forms 

    , := ∫_Ø^-1: + ∫_Ø: + ∫_Ø:

and

    ,  := ∫_Øρ^-1· + ,

and denote the Hilbertian product norm on ×Ł^2(Ø)^ by

    ^2 := ^2_ + ^2_0,Ø.


The variational formulation of the eigenvalue problem (<ref>) is given  as follows 
in terms of  κ:=1 + ω^2   (see <cit.> for more details):
Find κ∈ and ≠ (,)∈× such that

    (,),(,) = κ (,),(,)  ∀ (,)∈×.

 


We notice that the bilinear form 

    (, )_,  := ∫_Øρ^-1· +  ∫_Ø^-1:

also defines an inner product on . Moreover,  the following well-known result establishes that the norm 
induced by (·, ·)_, is equivalent to ·_ uniformly  in the Lamé coefficient λ.

There exist constants c_2≥ c_1>0 independent of λ such that 

    c_1 _≤_,≤ c_2 _  ∀∈,

where _,:= √((,))_,.


The bound from above follows immediately from the fact that 

    ∫_Ø^-1: = 1/2μ∫_Ø^:^ +1/n(nλ + 2μ)∫_Ø()()

is bounded by a constant independent of λ. The left inequality may be found, for example, in <cit.>.


As a consequence of Proposition <ref>,  
there exists a constant M>0 independent of λ such that 

    |, | ≤ M    ∀, ∈×.


 There exists a constant α>0, depending on ρ, μ and Ø (but not on λ), such that 
 
    sup_∈×, /≥α  ∀∈×.



It follows from Proposition <ref> that  

    (, 0), (, 0) = (,)_,≥ C_1^2 ^2_    ∀∈,

with C_1>0 independent of λ. 
On the other hand, there exists a constant β>0 depending only on Ø 
(see, for instance, <cit.>) such that 
 
    sup_∈∫_Ø:/_≥β_0,Ø    ∀∈.

Consequently, the Babuška-Brezzi theory  
shows  that, for any bounded linear form L∈ℒ(×), the 
problem: find ∈× such that 

    ,  = L(,)    ∀∈×

is well-posed, which proves  (<ref>).


We deduce from Proposition <ref> and from the symmetry of 
A(·, ·) that the operator : [Ł^2(Ø)^]^2 →× 
defined for any (, ) ∈ [Ł^2(Ø)^]^2, by 

    (, ),  = (, ),   ∀∈×

is well-defined and symmetric with respect to A(·, ·). Moreover, there exists a constant C>0 independent of 
λ such that 

    (, )≤ C (, )_0,Ø  ∀ (, ) ∈ [Ł^2(Ø)^]^2.

It is clear that  (κ,) 
is a solution of (<ref>) if and only if
( η=1/κ, ) is an eigenpair for . Let

    :=∈: =  Ø.

From the definition of , it is clear that 
|_×:×⟶×
reduces to the identity. Thus, η=1 is an eigenvalue of  with eigenspace ×. We introduce the orthogonal subspace to × in
× with respect to the bilinear form B,

    [×]^_B
    :={ (,)∈×:
      (,), (,)=0   ∀ (,)∈×}.




The subspace [×]^_B is invariant for
, i.e.,

    ([×]^_B)
    ⊂[×]^_B.

Moreover, we have the direct and stable decomposition

    × = [×]⊕[×]^_B.



See Lemma 3.3 and Lemma 3.4 of <cit.>.

We deduce from Lemma <ref> that there exists a unique projection : ×→×
with range [×]^_B and kernel ×
associated to the splitting (<ref>).

Let us consider the elasticity problem posed in 
Ø with a volume load in Ł^2(Ø)^n and with homogeneous Dirichlet and Neumann boundary conditions on Γ_D and  Γ_N, respectively. According to <cit.>, there exists ∈(0,1) that depends on Ø, λ and μ such that the displacement field that solves this problem belongs to ^̋1+s(Ø)^n for all s∈ (0,).
The following result shows that  and ∘ are regularizing operators. 



For all s∈ (0, ),
(×) ⊂^×^ and  
((×)) ⊂{^×^: ∈^̋1(Ø)^n}.
Moreover, there exists a constant C>0  such that 

    _^×^≤ C _0,Ø  ∀∈×

and

    ∘_×^≤ C _0,Ø  ∀∈×.



Estimate (<ref>) is proved in <cit.> and (<ref>) follows as a consequence of (<ref>), see <cit.>. 

We point out that, in principle, the exponent  and the constant C in (<ref>) depend on the Lamé coefficient λ. However, we know that (<ref>) also holds true when λ=+∞ (see the Appendix). Hence, it is natural to expect (<ref>) to be satisfied uniformly in λ. However, to the best of authors' knowledge, such a result is not available in the literature. For this reason, from now on we make the following assumption.

There exist ∈ (0,1) and C_0>0 independent of λ such that 

    _^̋s(Ø)^×^̋s(Ø)^≤C_0 _0,Ø  ∀∈×,  ∀ s∈ (0, ).


This would immediately imply the existence of C_1>0 independent of λ such that 

    ∘_×^≤C_1 _0,Ø  ∀∈×,  ∀ s∈ (0, ).

The next result gives the spectral characterization for the solution operator .

 The spectrum () of  decomposes as follows 
 
    () = 0, 1∪η_k_k∈ℕ

where η_k_k⊂ (0,1) is a real sequence of  
finite-multiplicity eigenvalues of   which converges to 0.  The ascent of each of
these eigenvalues is 1 and the corresponding eigenfunctions lie in (×).
Moreover, η=1 is an infinite-multiplicity eigenvalue of  with associated eigenspace 
× and  η=0 is not 
an eigenvalue. 


See <cit.>.


We end this section by providing a bound of the resolvent ( z - )^-1. 

	If z ∉(), 
	there exists a constant C>0 independent of λ and z such that 
	
    (z-) ( ,  )≥  C  (z,  () ) ( ,  )  ∀ ( ,  )∈×,
  
	where (z,  () ) represents the distance between z and  
	the spectrum  of  in the complex plane, which in principle depends on λ.
	

 
 See Proposition 2.4 in <cit.>.
 






§ A DISCONTINUOUS GALERKIN DISCRETIZATION



We consider  shape regular affine meshes 𝒯_h that subdivide the domain Ω̅ into  
triangles/tetrahedra K of diameter h_K. The parameter h:= max_K∈_h{h_K} 
represents the mesh size of _h. Hereafter, given an integer m≥ 0 and a domain 
D⊂ℝ^n, _m(D) denotes the space of polynomials of degree at most m on D.


We say that a closed subset F⊂Ω is an interior edge/face if F has a positive (n-1)-dimensional 
measure and if there are distinct elements K and K' such that F =K̅∩K̅'. A closed 
subset F⊂Ω is a boundary edge/face if
there exists K∈_h such that F is an edge/face of K and F =  K̅∩∂Ω. 
We consider the set _h^0 of interior edges/faces and the set _h^∂ of boundary edges/faces.
We assume that the boundary mesh _h^∂ is compatible with the partition = _D∪_N, i.e., 

    ⋃_F∈_h^D F = _D    and    ⋃_F∈_h^N F = _N,

where _h^D:= F∈_h^∂;    F⊂_D and 
_h^N:= F∈_h^∂;    F⊂_N.
We denote   

    _h := _h^0∪_h^∂    and    ^*_h:= _h^0∪_h^N,

and for any element K∈_h, we introduce the set 
 
    (K):= F∈_h;   F⊂∂ K

 of edges/faces composing the boundary of K.
 The space of piecewise polynomial functions of degree at most m relatively to _h is denoted by  

    _m(_h) := v∈ L^2(Ø);    v|_K ∈_m(K),  ∀ K∈_h .

For any k≥ 1, we consider the finite element spaces

    _h := _k(_h)^    _h^c := _h ∩    and    _h := _k-1(_h)^∩.


Let us now recall some well-known properties of the Brezzi-Douglas-Marini (BDM) 
mixed finite element <cit.>. For t>1/2, the tensorial version of the BDM-interpolation operator 
Π_h: ^̋t(Ø)^→_h^c,  
satisfies the following classical error estimate, see <cit.>,

    - Π_h _0,Ø≤ C h^min(t, k+1)_t,Ø    ∀∈^̋t(Ø)^,    t>1/2.

For less regular tensorial fields we also have the following error estimate

    - Π_h _0,Ø≤ C h^t (_t,Ø + _)    ∀∈^̋t(Ø)^∩,    t∈ (0, 1/2].

Moreover, thanks to the commutativity property, if ∈^̋t(Ø)^n, then

    ( - Π_h ) _0,Ø =  - ℛ_h _0,Ø≤ C h^min(t, k)_t,Ø,

where ℛ_h is the ^n-orthogonal projection onto _k-1(_h)^n. Finally, 
we denote by 𝒮_h: →_h the orthogonal
projector with respect to the ^-norm. 
It is well-known that, for any t>0, we have

    -𝒮_h_0,Ø≤ C h^min(t, k)_t,Ø    ∀∈^̋t(Ø)^∩.
 
For the analysis we need to decompose adequately the space ^c_h×_h.
We consider,

    _h = {∈^c_h;   = 0  }⊂.


There exists a projection _h: _h^c ×_h →_h^c ×_h with kernel 
_h ×_h such that for all s∈ (0,), there exists a constant C independent of h and λ such that

    ( - _h)≤ C  h^s _h_0,Ø  ∀∈_h^c ×_h.



See the proof of estimate (ii) of Lemma 4.2 from  <cit.> 








For any t≥ 0, we consider the broken Sobolev space 

    ^̋t(_h):=
     ∈Ł^2(Ø)^n;   |_K∈^̋t(K)^n  ∀ K∈_h.

For each :=_K∈^̋t(_h)^n and 
:= _K∈^̋t(_h)^
the components _K and _K  represent the restrictions |_K and |_K. 
When no confusion arises, the restrictions of these functions will be written 
without any subscript. We will also need the space given on the skeletons of the triangulations _h  by 

    Ł^2(_h):= ∏_F∈ℱ_hŁ^2(F).

Similarly, the components χ_F 
of χ := χ_F∈Ł^2(_h)  
coincide with the restrictions χ|_F and we denote

    ∫__hχ := ∑_F∈_h∫_F χ_F  and  χ^2_0,_h:= ∫__hχ^2, 
         ∀χ∈Ł^2(_h).

Similarly, χ^2_0,^*_h:= ∑_F∈^*_h∫_Fχ_F^2 for all χ∈Ł^2(^*_h):= ∏_F∈ℱ^*_hŁ^2(F). 

From now on, h_∈Ł^2(_h) is the piecewise constant function 
defined by h_|_F := h_F for all F ∈_h with h_F denoting the 
diameter of edge/face F.

Given a vector valued function ∈^̋t(_h)^n, with t>1/2, 
we define averages ∈Ł^2(_h)^n and jumps ∈Ł^2(_h) 
by

    _F := (_K + _K')/2   and  _F := _K ·_K + _K'·_K'  ∀ F ∈(K)∩(K'),

where _K is the outward unit normal vector to ∂ K. On the boundary of Ø we use the following 
conventions for averages and jumps:

    _F := _K    and  _F := _K ·  ∀ F ∈(K)∩.


Similarly, for matrix valued functions ∈^̋t(_h)^, we define ∈Ł^2(_h)^ and 
∈Ł^2(_h)^n by

    _F := (_K + _K')/2   and  _F := 
     _K _K + _K'_K'  ∀ F ∈(K)∩(K')

and on the boundary of Ω we set 

    _F := _K    and  _F := 
     _K   ∀ F ∈(K)∩.



Given ∈_h we define _h ∈Ł^2(Ø)^n by 
_h |_K =  (|_K) for all K∈_h and 
endow (h) :=  + _h with the seminorm

    ||^2_(h) := _h ^2_0,Ø + h_^-1/2^2_0,^*_h

and the norm

    ^2_(h) := ||^2_(h) + ^2_0,Ø.

For the sake of simplicity, we will also use the notation 

    (, )^2_DG : =  ^2_(h) + ^2_0,Ø.

The following result will be used in the sequel to ultimately derive a method free of spurious modes. Since according to Proposition <ref> the spectrum of  lies in the unit disk 𝔻:={ z∈ℂ:  |z|≤ 1}, we restrict our attention to this subset of the complex plane.



There exists a constant C>0 independent of h and λ such that for all  z ∈𝔻∖() with |z|≤ 1, there holds

    (z  - ) _DG≥ C(z,())|z| _DG  ∀∈(h)×.



We introduce

    (^*, ^*):= ∈×

and notice that 

    (z  - )(^*, ^*) =  (z  - ) .

By virtue of Proposition <ref> and the boundedness of :  [Ł^2(Ø)^]^2 →× we have that 

    C(z,()) (^*, ^*)≤(z  - )(^*, ^*)≤ (z  - ) 
    ≤(z  - ) _0 
    ≤(z  - ) _DG .

Finally, by the triangle inequality, 

    _DG≤ |z|^-1(^*, ^*) +    |z|^-1(z  - ) _DG
    ≤   |z|^-1( 1+C(z,())) (z  - ) _DG
    ≤    |z|^-1( C(z,())+C(z,()))(z  - ) _DG.

Hence, 

    C|z|((z,())/+(z,()))_DG≤(z  - ) _DG.

Since (z,())≤ |z|≤ 1  and ≤ C' (with C' independent of λ), we derive from the above estimate that 

    C|z|/1+C'(z,())_DG≤(z  - ) _DG,

and the result follows.























If E is a compact subset of 𝔻∖(), we deduce from Lemma <ref> that there exists a 
constant C>0 independent of h and λ such that, for all z∈ E,

    (z  - )^-1_ℒ((h)×, (h)×)≤C/(E,())|z|.



Let us now introduce the discrete counterpart of (<ref>). Given a parameter _S>0, we introduce  the symmetric bilinear form 

    ,:= 
     ∫_Øρ^-1_h ·_h  + , 
    
    + ∫_^*_h_S h_^-1 ·-∫_^*_h( ρ^-1_h· 
    + ρ^-1_h·)

and consider the DG method: Find κ_h∈ and 0≠∈_h×_h such that 

    , (_h,_h) = κ_h , (_h,_h)    ∀ (_h,_h)∈_h×_h.

We notice that, as it is usually the case for DG methods, the essential boundary condition is directly incorporated 
within the scheme. 

A straightforward application of the Cauchy-Schwarz  inequality shows that, for all , ∈^̋t(, _h)× (t>1/2), there exists a constant M^*>0 independent of h and λ such that 

    |,| ≤ M^* ^*_DG ^*_DG,

where 

    ^*_DG:= ( ^2_DG +  
    h_^1/2^2_0,^*_h)^1/2.

Moreover, we deduce from the discrete trace inequality (see <cit.>)

    h^1/2_v_0,_h≤ C v_0,Ø  ∀  v∈_k(_h),

that for all ∈^̋t(, _h)× (t>1/2), and (_h,_h) ∈_h×_h,

    |,(_h,_h)| ≤ M_DG^*_DG (_h,_h)_DG,

with M_DG>0 is independent of h and λ.






























































§ THE DG-DISCRETE SOURCE OPERATOR


















The following discrete projection operator from the DG-space _h onto the (̋,Ø)-conforming mixed finite element space ^c is essential in the forthcoming analysis.


There exists a projection ℐ_h: _h →_h^c such that 
the norm equivalence 

    C _(h)≤( ℐ_h ^2_ + 
     h_^-1/2^2_0,^*_h)^1/2≤C̅_(h)

holds true on _h with constants C>0 and C̅>0 independent of h. Moreover, we have that 

    _h (- ℐ_h )^2_0,Ø + ∑_K∈_h h_K^-2- ℐ_h ^2_0,K≤ C_0  
     h_^-1/2^2_0,^*_h,

 with C_0>0 independent of h.


See <cit.>.

We can prove, with the aid of this result, that the bilinear form A_h satisfies the following inf-sup condition that ensures the stability of our DG method.

There exists a positive parameter _S^* such that, for all _S≥_S^*,

    sup_∈_h×_h, /_DG≥α_DG_DG  ∀∈_h×_h
 
with α_DG>0 independent of h and λ.


It is shown in <cit.> that 
 there exists a constant α_A^c>0 independent of h and λ  such that 
 
    sup_∈_h^c×_h, /≥α_A^c  ∀∈_h^c×_h.

 It follows that there exists an 
  operator Θ_h: ^c_h×_h →^c_h×_h satisfying
  
    , Θ_h  = α_A^c^2   and  Θ_h ≤

  for all ∈^c_h×_h.
  
  Given ∈_h×_h, the decomposition _h = _h^c + _h, with   
  _h^c := ℐ_h _h and _h := _h - ℐ_h _h, and (<ref>) yield   
  
    , Θ_h +  
      = α_A^c   ^2 +
    ,  + 
      , Θ_h  + 
      , .
 
  By the Cauchy-Schwarz inequality, 
    , = ρ^-1_h _h_0,Ø^2 + 
       _S h_^-1/2_h_0,^*_h^2 + 
       ∫_Ø^-1_h:_h 
    
       - 2 ∫_^*_hρ^-1_h _h·_h  ≥  _S h_^-1/2_h_0,^*_h^2 
    
         -2 ρ^-1h_^1/2_h _h_0,^*_hh_^-1/2_h_0,^*_h

and we deduce from (<ref>) and (<ref>) that

    , ≥    (_S -C_1 ) h_^-1/2_h_0,^*_h^2 ,

with a constant C_1 independent of h and λ. 

We proceed similarly for the terms in the right-hand side of (<ref>).
Indeed, it is straightforward  that 

    , ≥ -ρ^-1_h^c_0,Ω_h _h_0,Ω - C_2  
     _h_0,Ω(_h^c_0,Ω+ _h_0,Ω)- 
    ρ^-1h_^1/2^c_h_0,^*_hh_^-1/2_h_0,^*_h,

and using again (<ref>) and  (<ref>)   we obtain 

    , ≥ 
     -C_3 h_^-1/2_h_0,^*_h≥
    
     - α_A^c/4^2 -  C_4h_^-1/2_h_0,^*_h^2

with C_4>0 independent of h and λ. Similar estimates  lead to  

    , Θ_h ≥ 
      -C_5  h_^-1/2_h_0,^*_hΘ_h ≥
    
      -C_5 h_^-1/2_h_0,^*_h,

where the last inequality follows from  (<ref>). We conclude that there exists 
C_6>0 independent of h and λ such that 

    , Θ_h ≥ 
      - α_D^c/4^2 -  C_6 h_^-1/2_h_0,^*_h^2.

We then have shown that,

    , Θ_h + ≥α_A^c  /2^2 + 
    (_S  -  
     C_7 ) 
     h_^-1/2_h_0,_h^2,

with C_7 := C_1+ C_4+C_6. Consequently,  
if _S > _S^*:=C_7 + α_A^c/2,  

    , Θ_h + ≥α_A^c  /2( ^2 + 
     h_^-1/2^2_0,^*_h),

and thanks to (<ref>),   we conclude that there exists α_DG>0 such that, 

    , Θ_h + ≥α_DG_DG(Θ_h  + _DG),
     
which gives (<ref>).

In the sequel, we assume that the stabilization parameter is big enough _S > _S^* so that the inf-sup condition (<ref>) is guaranteed. The first consequence of this inf-sup condition is that the operator _h: Ł^2(Ø)^×Ł^2(Ø)^→_h×_h 
characterized, for any (, ) ∈ [Ł^2(Ø)^]^2, by 

    _h(, ), (_h,_h) = (, ), (_h,_h)  ∀ (_h,_h)∈_h×_h

is well-defined, symmetric with respect to A_h(·, ·) and there exists a constant C>0 independent of 
λ and h such that 

    _h (, )_DG≤ C (, )_0,Ø  ∀ (, ) ∈ [Ł^2(Ø)^]^2.



We observe that if (κ_h,)∈ℝ×_h× is a solution of problem (<ref>) if and only if (μ_h,), with μ_h=1/(1+κ_h) is an eigenpair  of _h, i.e.

    _h=1/1+κ_h.





Analogously to the continuous case, we prove that the discrete resolvent associated to the discrete operator _h is bounded.

Assume that (, ):=(, )∈^̋t(, Ø)×^̋t(Ø)^ for some t>1/2.
Then, 

    ( - _h)(, )_DG≤(1 + M_DGα_DG) 
    inf_∈_h×_h(, ) - ^*_DG.





Moreover, the error estimate  

    ( - _h)(, )_DG≤  C   h^min(t, k) ( _^̋t(, Ø) +  _^̋t(Ø)^),

holds true with a constant  C>0 independent of h and λ.


We first notice that the DG approximation (<ref>) is consistent with regards to its continuous counterpart (<ref>) in the sense that 

    (-_h)(, ),  = 0   ∀∈_h×_h.

Indeed, by definition, 
 
    (, ),  = ∫_Øρ^-1·_h _h + (, ), 
     -
    ∫_^*_hρ^-1·_h.

 It is straightforward to deduce from (<ref>) 
 
 
    ∇(ρ^-1) = ^-1(- ) + -   and   
     (- ^)/2 = ( - ^)/2.

Moreover, an  integration by parts  yields

    ∫_Øρ^-1·_h _h  = 
     -∑_K∈_h∫_K ∇(ρ^-1): _h + ∑_K∈_h∫_∂ Kρ^-1·_h_K
    
     = -∑_K∈_h∫_K ∇(ρ^-1): _h + ∫_^*_hρ^-1·_h.

Substituting back the last identity and (<ref>) into  (<ref>) 
we obtain 

    (, ),  = (, ),   ∀∈_h×_h

and (<ref>) follows.


The Céa estimate (<ref>) follows now in the usual way by taking advantage of (<ref>), the inf-sup condition (<ref>), estimate  (<ref>), and the triangle inequality.

It follows from (<ref>) that  

    ( - _h)(, )_DG≤(1 + M_DGα_DG) 
     (, ) - (Π_h, 𝒮_h)^*_DG.

Using the interpolation error estimates (<ref>), (<ref>) and (<ref>) we immediately obtain 

    (, ) - (Π_h, 𝒮_h)_DG = (, ) - (Π_h, 𝒮_h)≤ 
    C_0   h^min(t, k) ( _^̋t(, Ø) +  _^̋t(Ø)^).

Moreover, we notice that 

    h_^1/2 (-Π_h)_0,^*_h≤∑_K∈_h∑_F∈(K) h_F (-Π_h) ^2_0,F.

Under the regularity hypotheses on  , the commuting diagram property satisfied by Π_h, the trace theorem and standard scaling arguments give 

    h_F^1/2 (-Π_h) _0,F = h_F^1/2- ℛ_K _0,F≤ C_2 h_K^min(t,k)_t,K

for all F∈(K), where the Ł^2(K)-orthogonal projection ℛ_K:= ℛ_h|_K 
onto _k-1(K) is applied componentwise. It follows that 

    h_^1/2 (-Π_h)_0,^*_h≤ C_3 h_K^min(t,k)( ∑_K∈_h_t,K^2 )^1/2≤ 
    	C_3 h_K^min(t,k)_t,Ø.

Combining  (<ref>) and (<ref>) with (<ref>) proves the asymptotic error estimate (<ref>).



For all s∈ (0,), there exists a constant  C>0 independent of h and λ, such that for all ∈×

    ( - _h)_DG≤  C   h^s _0,Ω.



The result is a consequence of Theorem <ref> by noticing that, by virtue of Lemma <ref> and Assumption <ref>, ∘⊂{(,)∈[^̋s(Ø)^n× n]^2: ∈^̋1(Ø)^n} for all s∈(0,).



For all s∈ (0,), there exists a constant C>0 independent of h and λ such that 

    ( - _h)_DG≤ C   h^s  _DG  ∀∈_h×_h.



For any _h∈_h we consider the splitting  _h = _h^c + _h with 
_h^c:=ℐ_h _h∈_h^c. We have that 

    ( - _h) = ( - _h)(_h,) + ( - _h)(_h^c, _h)
     = ( - _h)(_h,) +
    ( - _h)_h(_h^c, _h),

where the last identity is due to the fact that ( - _h)(_h^c, _h)∈_h×_h 
and  - _h vanishes identically on this subspace. It follows that 

    ( - _h) = ( - _h)(_h,) +
    ( - _h)(_h - )(_h^c, _h) + ( - _h)(_h^c, _h),

and the triangle inequality together with (<ref>) and  (<ref>) yield

    ( - _h)_DG≤( - _h)(_h,)_DG +
    ( - _h)(_h - )(_h^c, _h)_DG
     + ( - _h)(_h^c, _h)_DG≤(_ℒ([Ł^2(Ø)^]^2, ×) + 
    _h_ℒ([Ł^2(Ø)^]^2, _h×_h))
    ( _h_0,Ø  + (_h - )(_h^c, _h)) + 
     ( - _h)(_h^c, _h)_DG.

Using (<ref>), Lemma <ref>, Assumption <ref> and Lemma <ref> we have that

    _h_0,Ø≤ C h _h_(h),


    (_h - )(_h^c, _h)≤ C h^s _h^c_0,Ø≤ C h^s _h_(h)

and 

    ( - _h)(_h^c, _h)_DG≤ C h^s _h^c_0,Ø≤ C h^s _h_(h)

respectively, which gives the result.





§ SPECTRAL CORRECTNESS OF THE DG METHOD



The convergence analysis follows the same steps introduced in <cit.>, we only need to adapt 
it to the DG context, cf. also <cit.>.  

For the sake of brevity, we
will denote in  this section :=×, _h:=_h×_h
and (h) := (h)×. Moreover, when no confusion can arise, we
will use indistinctly , , etc. to denote elements in 
and, analogously, _h, _h, etc. for those in _h.
Finally, we will use ·_ℒ(_h,  (h))
to denote the norm of an operator restricted to the discrete
subspace _h; namely, if :(h)→(h), then

    _ℒ(_h, (h)):=sup_≠_h∈_h_h_DG/_h_DG.



If z ∈𝔻∖(), there exists h_0>0 such that if h≤ h_0,

    (z  - _h) _h _DG≥ C(z,())|z| _h_DG  ∀_h ∈_h.

with C>0 independent of h and λ.


It follows from 

    (z  - _h) _h = (z  - ) _h + ( - _h ) _h

and Lemma <ref> that 

    (z  - _h) _h_DG≥(C(z,())|z| -  - _h_ℒ(_h, (h))) _h_DG

and the result follows from Lemma <ref>.














































 
If z ∈𝔻∖(), there exists h_0>0 such that if h≤ h_0,

    (z  - _h) _DG≥ C (z,())|z|^2 _DG  ∀∈(h),

with C>0 independent of h and λ.



Given ∈(h) we let 

    _h^* = _h ∈_h.

We deduce from the identity 

    (z  - _h) _h^* = _h (z - _h)

and from Lemma <ref> that

    C(z,())|z| _h^* _DG≤(z  - _h) _h^*_DG≤_h_ℒ((h), _h)(z  - _h)_DG.

This and the triangle inequality leads to

    _DG≤ |z|^-1^*_h_DG + |z|^-1 (z  - _h)_DG
    ≤|z|^-1( 1 + _h_ℒ((h), _h)/C(z,())|z|) (z  - _h)_DG.
    
    ≤|z|^-1( C(z,())|z|+_h_ℒ((h), _h)/C(z,())|z|) (z  - _h)_DG.

Hence,

    C|z|(C(z,())|z|/_h_ℒ((h), _h)+C(z,())|z|)_DG≤(z  - _h) _DG.

Now, using that (z,())≤ |z|≤ 1  and _h_ℒ((h), _h)≤ C' (with C' independent of λ), from the estimate above we derive

    C|z|^2 (z,())_DG≤(z  - ) _DG,

and the result follows. 





If E is a compact subset of 𝔻∖() and h is small enough, we deduce from Lemma <ref> that (z-_h):(h)→(h) is invertible for all z∈ E. Hence, E⊂𝔻\(_h). Consequently, for h small enough, the numerical method does not introduce spurious eigenvalues. Moreover, we have that there exists a constant C>0 independent of h and λ such that, for all z∈ E,

    (z  - _h )^-1_ℒ((h),(h))≤C/(E,())|z|^2.








For
∈(h) and 𝔼 and 𝔽 closed subspaces of (h), we set
δ(,𝔼):=inf_∈𝔼-_DG,
δ(𝔼,𝔽):=sup_∈𝔼: =1δ(,𝔽),
and (𝔼,𝔽):=maxδ(𝔼,𝔽),δ(𝔽,𝔼),
the latter being the so called gap between subspaces 𝔼 and
𝔽.

Given an isolated eigenvalue κ≠ 1  of , we define 

    _κ:=1/2(κ,()∖{κ}).

 It follows that  the closed disk D_κ:={z∈ℂ:   |z-κ|≤_κ} of the complex plane, with center κ and boundary γ is such that   
D_κ∩() = κ.  We deduce from Remark <ref> that the operator 
:=1/2π i∫_γ(z-)^-1  dz:(h)⟶(h) is well-defined and bounded uniformly in 
h. Moreover, |_ is a spectral projection in  onto the (finite dimensional)  eigenspace () corresponding to the eigenvalue κ of . In fact, 

    ((h)) = ().

To prove this, let κ^*∈ D_κ be an eigenvalue of : (h)→(h) and 
^*∈(h) be the corresponding eigenfunction. Since κ^*≠ 0 and ((h))⊂, we actually have that ^*∈. Then, necessarily,  κ^*=κ and taking into account that () is the eigenspace associated with κ we deduce (<ref>). 

Similarly, we deduce from Remark <ref> that, for h small enough, the operator 
_h:=1/2π i∫_γ(z-_h)^-1  dz:(h)⟶(h) is also well-defined and bounded uniformly in h. Moreover, _h|__h is a projector in _h onto the eigenspace _h(_h)  corresponding to the eigenvalues of _h: _h →_h contained in γ. The same arguments as above show that we also have, 

    _h((h)) = _h(_h).

Our aim now is to compare _h(_h) to () in terms of the gap . In order to do that,  we assume the following regularity assumption ()⊂^̋t(,Ø)×^̋t(Ø)^n× n with t>s. 


There exists C>0, independent of h and λ, such that

    - _h_ℒ(_h, (h))≤C/_κ - _h_ℒ(_h, (h)).



We deduce from the identity 

    (z-)^-1-(z-_h)^-1 = (z-)^-1( - _h)(z-_h)^-1

that, for any _h∈_h,


    ( - _h)_h_DG≤1/2π∫_γ[(z-)^-1-(z-_h)^-1]_h_DG|dz|
    
    =1/2π∫_γ[(z-)^-1( - _h)(z-_h)^-1]_h_DG|dz|
    ≤1/2π∫_γ(z-)^-1_ℒ((h), (h)) - _h_ℒ(_h, (h))(z-_h)^-1_ℒ(_h, _h)_h_DG|dz|

and the result follows from Lemmas <ref> and <ref>, the definition (<ref>) and the fact that for all z∈γ, |z|≥κ-_κ≥1/2κ.




There exists a constant C>0 independent of h and λ such that

    ((), _h(_h)) ≤ C (  - _h_ℒ(_h, (h))/_κ +
    δ((), _h)).



As _h is a projector,  for h sufficiently small,  we have that _h_h=_h for all _h∈_h(_h).  It follows from (<ref>) that _h ∈(), which leads to

    δ(_h,())
    ≤_h_h-_h_DG≤_h-_ℒ(_h, (h))_h_DG

for all _h∈_h(_h). We deduce from (<ref>)  that 

    δ( _h(_h), ())≤C/_κ - _h_ℒ(_h, (h)).


On the other hand,  as =
for all ∈(), for h small enough and _h∈_h,

    -_h_h_DG≤(-_h)_DG
    +(-_h)_h_DG≤
    _ℒ( (h), (h))(-_h)_DG
    +(-_h)_ℒ(_h, (h))_h_DG
    ≤(_h_ℒ((h), (h))+2_ℒ((h), (h))) - _h_DG
    +-_h_ℒ(_h, (h))_DG.

Consequently,

    δ(, _h(_h)) ≤ C (δ(, _h) + -_h_ℒ(_h, (h)) )

for all ∈() such that _DG = 1 and using that the eigenspace () is finite dimensional we deduce that

    δ((), _h(_h)) ≤ C (δ((), _h) + -_h_ℒ(_h, (h)) )

and the result follows from the last estimate and (<ref>).


Let κ≠ 1 be an eigenvalue of  of algebraic multiplicity m  and 
let D_κ be a closed disk in the complex plane centered at κ with  boundary γ such that 
D_κ∩() = κ. Let κ_1, h, …, κ_m(h), h be the eigenvalues of 
_h: _h →_h  lying in  D_κ and repeated according to their algebraic multiplicity. 
Then, we have that m(h) = m for h sufficiently small and 

    lim_h→ 0max_1≤ i ≤ m |κ -  κ_i, h| =0.

Moreover, if () is the eigenspace corresponding to κ and _h(_h) is the
_h-invariant subspace of _h spanned by the eigenspaces  corresponding to
κ_i, h, i = 1,…, m then

    lim_h → 0((), _h(_h)) = 0.



We deduce from Lemma <ref> that 

    lim_h→ 0 - _h_ℒ(_h, (h)) = 0.
 
Moreover, as () ⊂^̋t(, Ø)×^̋t(Ø)^, it follows from (<ref>) that




    lim_h→ 0δ((), _h) = 0.

Hence, by virtue of Theorem <ref>, we have that 

    lim_h→ 0((), _h(_h)) = 0,

and, as a consequence, () and _h(_h) have the same dimension provided h is sufficiently small. Finally, being κ an isolated eigenvalue and the radius of the circle γ arbitrary,  we deduce that 

    lim_h→ 0max_1≤ i ≤ m |κ -  κ_i, h| =0.




§ ASYMPTOTIC ERROR ESTIMATES



Along this section we fix a particular eigenvalue κ≠ 1 of . We wish to obtain error 
estimates for the eigenfunctions and the eigenvalues in terms of the quantity

    δ^* (() , _h):= sup_∈(), =1inf__h∈_h - _h^*_DG.
 









For h small enough, there exists a constant C independent of h such that 

    ((), _h(_h) )≤C/_κδ^* (() , _h).

 

As ((h))=() and _h((h))=_h(_h), it is equivalent to show that 

    (((h)), _h((h)))≤C/_κδ^* (() , _h).

We consider here again the disk D_κ centered at κ with radius _κ and boundary γ. We first notice that for all z∈γ 

    (z-)^-1-(z-_h)^-1 = (z-_h)^-1( - _h)(z-)^-1,

which implies

    ( - _h)|_()≤1/2π∫_γ(z-)^-1-(z-_h)^-1|_()|dz|
    
    =1/2π∫_γ(z-_h)^-1( - _h)(z-)^-1|_()|dz|
    ≤1/2π∫_γ(z-_h)^-1_ℒ((h), (h))( - _h)|_()_ℒ(, (h))(z-)^-1_ℒ(, (h))|dz|
    ≤C/_κ( - _h)|_()_ℒ(, (h))






Now, on the one hand, it is clear that

    δ(((h)), _h((h))) ≤ ( - _h)|_()_ℒ(, (h)).

On the other hand,  (<ref>), the Céa estimate given by Theorem <ref> and the fact that  
() is finite dimensional yield 

    ( - _h)|_()_ℒ(, (h))≤C/_κδ^* (() , _h),

which proves that 

    δ(((h)), _h((h)))   ≤C/_κδ^* (() , _h).

Consequently, as ()⊂^̋t(, Ø)×^̋t(Ø)^, we have that 

    lim_h→ 0δ(((h)) , _h((h))) = 0.

It is shown in <cit.> that (<ref>) 
implies that, for h small enough, Λ_h:=_h|_(): () →_h((h)) 
is bijective and Λ_h^-1 exists and is uniformly bounded with respect to 
h. Furthermore, it holds that,

    sup__h ∈_h((h)), _h_DG = 1Λ_h^-1 - _DG≤ 
    2 sup_∈((h)), _DG = 1Λ_h  - _DG.

Hence,

    δ(_h((h)), ((h))) ≤sup__h ∈_h((h)), _h_DG = 1_h - Λ_h^-1_DG≤  2 sup_∈(), _DG = 1 - _h _DG,

and (<ref>) shows that we also have δ( _h((h)), ((h)))≤C/_κδ^* (() , _h), 
and the result follows from this last estimate and (<ref>).






































Assume that  () ⊂^̋t(,Ø)×^̋t(Ø)^, then there exists C>0 independent of h and λ such that 

    (_h(_h), ())≤C/_κ h^min{t, k}.

Moreover, there exists C'>0 independent of h such that

    max_1≤ i≤ m |κ - κ_i, h|≤C'/_κ  h^2min{t, k}



Using the estimate (<ref>) from the last theorem and proceeding as in the proof of (<ref>) we immediately obtain  (<ref>).

Let κ_1, h, ⋯, κ_m, h be the eigenvalues of 
_h: _h →_h  lying in  D_κ and repeated according to their algebraic multiplicity.
 We denote by 
_i,h the eigenfunction  corresponding to κ_i,h and satisfying 
_i,h_DG=1. We know from Theorem <ref> that, if h is sufficiently small,

    δ(_i,h,())≤C/_κδ^*((), _h).

Then, there exists an eigenfunction :=(, ) ∈() satisfying


    _i,h-_DG = δ(_i,h,()) ≤(_h(_h), ()) ≤C/_κδ^*((), _h)→ 0   as h→ 0,

which proves that _DG is bounded from below and above by constant independent of h.
Proceeding as in the proof of the consistency property in Theorem <ref> we readily  obtain that

    A_h(, _h) = κ B(, _h)

for all _h:=(_h, _h)∈_h.
With the aid of (<ref>), it is easy to show that the identity

    A_h(-_i,h,-_i,h)
    - κ B(-_i,h,-_i,h)
    =(κ_i,h-κ) B(_i,h,_i,h)

holds true. Now, according to Lemma 3.6 of <cit.>, for any ∈(), ≠ 0, it holds that B(,)>0.Thus, since  () is finite-dimensional,
there exists c>0, independent of h, such that B(,)≥ c_DG.
This proves that B(_ih,_ih)≥c/2  for h sufficiently
small. We obtain from (<ref>) that  

    c/2κ_i,h-κ≤A_h(-_i,h,-_i,h) + |κ| B(-_i,h,-_i,h)≤ C (-_i,h_DG^*)^2.

Since :=(,) and _i,h:=(_h,_h), and  by definition of  ·_DG^* we have

    -_i,h_DG^*:=(,)-(_h,_h)_DG^*=(,)-(_h,_h)_DG+h_^1/2(-_h)__h^*.

It follows from Theorem <ref>, Lemma <ref> and the interpolation error estimates (<ref>)-(<ref>) that 

    (,)-(_h,_h)_DG≤ C_0 ((), _h(_h)) ≤ C_1 h^min{t, k}( 1 + _^̋t(, Ø) +  _^̋t(Ø)^).

On the other hand, 

    h_^1/2(-_h)__h^*≤h_^1/2(-Π_h)__h^*+h_^1/2(Π_h-_h)__h^*

and it follows from (<ref>) that 

    h_^1/2(-Π_h)__h^*≤ C_2 h_K^min(t,k)_t,Ø.

Finally,  using (<ref>), (<ref>) and (<ref>) yield 

    h_^1/2(Π_h-_h)__h^*   ≤ C_3(Π_h-_h)_0,Ø
       ≤ C_3
    ((Π_h-)_0,Ø+_h(-_h)_0,Ø)
       ≤ C_3 ((Π_h-)_0,Ø+(,)-(_h,_h)_DG) 
       ≤ C_4h^min{t,k}( 1 + _^̋t(, Ø) +  _^̋t(Ø)^).

Combining (<ref>), (<ref>)-(<ref>) and (<ref>), we obtain (<ref>).

 

In the proof provided above for the error estimate (<ref>) the constant C' is not independent of λ. Indeed, according to the proof of Lemma 3.6 from <cit.>, we have that 

    B((,),(,))=∫_Ø𝒞^-1:≥min{n/nλ+2μ,1/2μ}^2_0,Ø≥ 0.

Therefore,  the constant c in the proof above tends to zero when λ goes to infinity. However, the numerical experiments presented below suggest that (<ref>) holds true uniformly in λ. 


We notice that there is in  (<ref>) and (<ref>) a hidden reliance on λ through the constant _κ:=1/2(κ,()∖{κ}) because () depends on λ. The constant _κ measures the deterioration of the error estimates given in Theorem <ref> when the eigenvalue κ is too close to the accumulation point 0. 



We point out that, thanks to Lemma <ref>, we always have that () ⊂(,) ∈ [^]^2: ∈^̋1(Ø)^n for all s∈ (0,). Consequently, the error estimates given in Theorem <ref> will always hold true for any t∈(0,) even if ≤ 1/2. However, it may happen that some eigenspaces satisfy the regularity assumption of the theorem with t ≥. 












































































































































§ NUMERICAL RESULTS


We present a series of numerical experiments to solve the elasticity eigenproblem in mixed form with the discontinuous Galerkin scheme (<ref>).  All the numerical results have been obtained by using the FEniCS Problem Solving Environment <cit.>. For simplicity we consider a two-dimensional model problem. We choose Ω = (0,1)×(0,1),  ρ = 1, and a Young modulus E=1. We will let the Poisson ratio ν take different values in (0,1/2]. We recall that the Lamé coefficients are related to E and  ν by

    λ:=Eν/(1+ν)(1-2ν)    and    μ:=E/2(1+ν).

The limit problem corresponding to λ=∞ is obtained by taking ν = 1/2. In all our experiments we used uniform meshes with the symmetry pattern shown in Figure <ref>. The refinement parameter N represents the number of elements on each edge.




In the first tests we are concerned with the determination of a reliable  stabilization parameter _S. We know that the spectral correctness of the method can only be guaranteed if _S is sufficiently large (Proposition <ref>) and if the meshsize h is sufficiently small (cf. Remark <ref>). In a first stage, we fix the refinement level to N=8 and report in Tables <ref>, <ref> and <ref> the 10 smallest vibration frequencies computed for different values of _S. The polynomial degrees are given by k=3, 4, 5, respectively. The boxed numbers are spurious eigenvalues. We observe that they emerge at random positions when we vary _S and k and they disappear completely when _S is sufficiently large. 

 






Next, we present in Table <ref> different approximations of the first 10 vibration frequencies corresponding to N = 8, 16, 32, 64, obtained with _S =20 and a polynomial degree k = 3.   We notice that as the level of refinement increases the lower frequencies are progressively cleaned from spurious modes. We conclude that our method provides a correct approximation of the spectrum as long as N and _S are large enough. In the forthcoming tests we will take _S=1000. We point out that the previous tests have been carried out with a Poisson ratio ν=0.35, but similar results were obtained for values ranging from 0.35 to 0.5. 



The subsequent numerical tests are aimed to determine the convergence rate of the scheme. With the boundary conditions considered in our model problem, it turns out that (cf. <cit.> and the references therein) the regularity exponents  defined in Lemma <ref> are given by Table <ref> for different values of the Poisson ratio ν. 



We present in Tables <ref>, <ref> and <ref> (corresponding to the polynomial degrees k=2,3,4, respectively) the first two vibration frequencies computed on a series of nested  meshes for a range of Poisson ratios given by ν=0.35, 0.49, 0.5. We also report in these tables  an estimate  of the order of convergence α, as well as more accurate approximations of the vibration frequencies obtained by means of the least-squares fitting technique explained in <cit.>.  Comparing with the exponents given in Table <ref>, we observe that our method provides a double order of convergence for the vibration frequencies. Namely, in all cases  we have α≃ 2, which corresponds to the the worst possible order of convergence. The eigenfunctions corresponding to higher natural frequencies are oscillating  but they can be more regular (see Remark <ref>), which justifies the use of high polynomial orders of approximation.  Finally, we point out that the method is clearly locking-free.









§ APPENDIX. THE LIMIT PROBLEM


As was shown in the previous section, the proposed method works fine also for the limit problem (λ=+∞), namely, for perfectly incompressible elasticity. In this appendix, we will establish a spectral characterization in this case. Also, we will prove that the eigenvalues of the nearly incompressible elasticity problem converge to those of the incompressible elasticity problem as λ→∞.

In the limit case λ=+∞, the bilinear forms A and B change in their definitions, since the term where λ appears in (<ref>) vanishes. Therefore,  the limit eigenvalue problem reads as follows: Find κ∈ℝ and (,)∈× such that

    A_∞((,),(,))=κ B_∞((,),(,))    ∀(,)∈×

with

    B_∞((,),(,)):=1/2μ∫_Ø^:^+∫_Ø:+∫_Ø:

and

    A_∞((,),(,)):=∫_Øρ^-1· +B_∞((,),(,))

for all (,),(,)∈×.

It is easy to check that A_∞ is a bounded bilinear form. Moreover, the arguments used in the proofs of Propositions <ref> and <ref> hold true for λ=+∞, so that  A_∞ satisfies the following inf-sup condition:

    sup_∈×A_∞((,),(,))/≥α(,)    ∀(_,) ∈×.

In consequence, we are in a position to introduce a solution operator for the limit eigenvalue problem. Let _∞: [Ł^2(Ø)^]^2 →× be
defined for any (, ) ∈ [Ł^2(Ø)^]^2  by

    A_∞(_∞(,),(,))= B_∞((,),(,))    ∀(,)∈×.


It is easy to check that μ is a non-zero eigenvalue of _∞ with eigenfunction (_∞,_∞) if and only of κ = 1/μ is a non-vanishing eigenvalue of problem (<ref>) with the same eigenfunction.

Our first goal is to prove that the operators  defined by (<ref>) converges to _∞ as λ goes to infinity. To recall that  actually depends on λ,   in what follows we will denote it by _λ.

Before proving the convergence of _λ to _∞, we will characterize the spectrum of _∞. 
Let  be defined as in (<ref>) and 

    [×]^_B_∞
    :={ (,)∈×:
     B_∞( (,), (,))=0   ∀ (,)∈×}.

We observe that _∞|_×:×→× reduces to the identity, so that μ=1 is  an eigenvalue of _∞. Moreover, its associated eigenspace is precisely ×.

Let us introduce the following operator which will play a role similar to that of  in the limit problem:

    _∞:×   →×,
    
          (,)   ↦_∞:=(, ).

where (,(,))∈×[Ł^2(Ø)^n×] is the solution of the following problem:

    1/2μ   ∫_Ø^:^+∫_Ø·+∫_Ø:=0    ∀∈,
       ∫_Ø·+∫_Ø:=∫_Ø·    ∀(,)∈Ł^2(Ø)^n×.


The previous problem is well posed, since the ellipticity of ∫_Ø^:^ in the corresponding kernel is established in Lemma 2.3 of <cit.> and the following inf-sup condition holds true (see <cit.>):


    sup_∈∫_Ø·+∫_Ø:/_(̋,Ø)≥β(_0,Ø+_0,Ø)    ∀(,)∈Ł^2(Ø)^n×.


We observe that problem (<ref>)–(<ref>) is a dual mixed formulation with weakly imposed symmetry
of the following incompressible elasticity problem with volumetric force density -

    -   =-in Ø,
    1/2μ^   =ε()in Ø,
       =on Γ_N,
       =on Γ_D.


It is easy to check that (,)∈(̋,Ø)×^̋1(Ø)^n satisfies (<ref>)–(<ref>) if and only if (,(,))∈×[Ł^2(Ø)^n×] is the solution of (<ref>)–(<ref>) with =1/2[∇-(∇)^].

Now, by resorting to the relation between the incompressible elasticity and the Stokes problems, we conclude that there exists _∞∈ (0,1) depending only on Ø and μ (see for instance <cit.>) such that, for all s∈ (0,_∞)  the solution  of (<ref>)–(<ref>) belongs to ^̋1+s(Ø)^n and the following estimate hold true 

    _1+s,Ø≤ C_0,Ø,

with a constant C independent of .

The following lemma is a consequence of this regularity result.

For all s∈ (0,) and (,)∈×, if (,(,)) is the solution of (<ref>)–(<ref>), then  ∈^̋s(Ø)^n× n, ∈^̋1+s(Ø)^n× n, ∈^̋s(Ø)^n× n and 

    _s,Ø+_1+s,Ø+_s,Ø≤ C_0,Ø,

with a constant C independent of . Consequently, _∞(×)⊂^̋s(Ø)^n× n×^̋s(Ø)^.


We observe that _∞ is idempotent  and that (_∞)=×. Moreover, being _∞  a projector, the orthogonal decomposition ×=(×)⊕_∞(×) holds true. On the other hand, _∞(×) is an invariant space of _∞ (see Proposition A.1 in <cit.>).  


For all s∈ (0,) 

    _∞(_∞(×))⊂(^*,^*)∈^̋s(Ø)^×^̋s(Ø)^: ^*∈^̋1(Ø)^n,

and there exists C>0 such that for all (,)∈_∞(×), if (^*,^*)=_∞(,), then

    ^*_s,Ø+^*_1,Ø+^*_s,Ø≤ C(,).

Moreover, _∞|__∞(×):_∞(×)→_∞(×) is a compact operator.


Let (,)∈_∞(×) and (^*,^*)=_∞(,). Hence, we have 

    ∫_Øρ^-1^*·+1/2μ∫_Ø^*:^+∫_Ø^*:=1/2μ∫_Ø^:^+∫_Ø:  ∀∈,
       ∫_Ø^*:=∫_Ø:    ∀∈.

Then, testing the first equation of the system above with ∈𝒟(Ø)^⊂, we have that

    -ρ^-1∇(^*)+1/2μ^*+^*=1/2μ^+.

Hence, since ρ and μ are constants, we conclude that ^*∈^̋1(Ø)^n.

Since _∞(×) is invariant with  respect to _∞, applying  Lemma <ref> we obtain directly (<ref>). On the other hand, (<ref>) is a consequence of Lemma <ref>. Finally,  the compactness of _∞|__∞(×) is a consequence of the following compact embedding

    (^*,^*)∈^̋s(Ø)^×^̋s(Ø)^:^*∈^̋1(Ø)^n↪×,

which allow us to conclude the proof.
 
Now we are in position to establish a spectral characterization for             _∞.

The spectrum of _∞ decomposes as follows: (_∞)=0,1∪μ_k_k∈ℕ, where:


  (i) μ=1 is an infinite-multiplicity eigenvalue of _∞ and its associated eigenspace is ×.

  (ii) μ=0 is an eigenvalue of _∞ and its associated eigenspace is 𝒵×, where

    𝒵:=∈:^=0=q:q∈^̋1(Ø)andq=0  onΓ_N.


  (iii) μ_k_k∈ℕ⊂(0,1) is a sequence of nondefective finite-multiplicity eigenvalues of _∞ which converge to zero and the corresponding eigenspaces lie in _∞(×).



It is enough to follow the steps of Theorem 3.5 from <cit.>.



Now we are in position to establish the following convergence result.

There exists a constant C>0 such that 

    (_λ-_∞)((,))≤C/λ(,)_0,Ø    ∀(,)∈[Ł^2(Ø)^]^2.



Let (,)∈[Ł^2(Ø)^]^2 and let (_λ,_λ):=_λ(,) and (_∞,_∞):=_∞(,). Then, from (<ref>) and the definition of  we have

    ∫_Øρ^-1_λ·+1/2μ∫_Ø_λ^:^+   1/n(nλ+2μ)∫_Ø(_λ)()+∫_Ø_λ:
       =1/2μ∫_Ø^:^+1/n(nλ+2μ)∫_Ø()()+∫_Ø:,
    ∫_Ø_λ:=∫_Ø:.

Whereas

    ∫_Øρ^-1_∞·+1/2μ∫_Ø_∞^:^+∫_Ø_∞:=1/2μ∫_Ø^:^+∫_Ø:  ∀∈,
       ∫_Ø_∞:=∫_Ø:    ∀∈.

Subtracting the  above equations we have 

    ∫_Øρ^-1   (_λ-_∞)·+1/2μ∫_Ø(_λ^-_∞^):^
       +∫_Ø(_λ-_∞):
    =1/n(nλ+2μ)∫_Ø(-_λ)()    ∀∈,
    ∫_Ø(_λ-   _∞):=0    ∀∈.

Testing this equation with :=_λ-_∞ and :=_λ-_∞ we have 






    ρ^-1(_λ-_∞)_0,Ø^2+1/2μ_λ^-_∞^_0,Ø^2 =   1/n(nλ+2μ)∫_Ø(()-(_λ))(_λ-_∞)
    ≤   1/n(nλ+2μ)∫_Ø()-(_λ)_0,Ø (_λ-_∞)_0,Ø
    ≤   1/nλ+2μ(_0,Ø+_λ_0,Ø)_λ-_∞_0,Ø
    ≤   C/nλ(,)_0,Ø_λ-_∞_0,Ø,

where we have used (<ref>) to bound _λ_0,Ø. Moreover

    min{ρ^-1,1/2μ}_C_ρ,μ((_λ-_∞)_0,Ø^2+_λ^-_∞^_0,Ø^2) ≤C/nλ(,)_0,Ø_λ-_∞_0,Ø.

 We observe that  (_λ-_∞)∈ is symmetric due to equation (<ref>). Then, we resort to the following estimate (see <cit.> for instance)

    C_λ-_∞_0,Ø^2≤_λ^-_∞^_0,Ø^2+(_λ-_∞)_0,Ø^2

with C>0 to deduce that

    C_λ-_∞_(̋,Ø)≤( _λ^-_∞^_0,Ø^2+(_λ-_∞)_0,Ø^2)^1/2.

Hence 

    (_λ-_∞)_0,Ø^2 +   _λ^-_∞^_0,Ø^2 
       ≤C_ρ,μ/nλ(,)_0,Ø( ^-_∞^_0,Ø^2+(_λ-_∞)_0,Ø^2)^1/2

and, finally,

    _λ-_∞_(̋,Ø)≤C/λ(,)_0,Ø,

with C  a positive constant depending on ρ, μ and n.

On the other hand, taking into account the inf-sup condition (<ref>),  (<ref>),  Cauchy-Schwarz inequality, (<ref>) and (<ref>), we have 
 
    β   _λ-_∞_0,Ø
       ≤sup_∈1/n(nλ+2μ)∫_Ø(_λ-_∞)()-∫_Øρ^-1(_λ-_∞)·-1/2μ∫_Ø(_λ^-^_∞):^/_
       ≤sup_∈C/nλ+2μ(,)_0,Ø_0,Ø+ρ^-1(_λ-_∞)_0,Ø_0,Ø+1/2μ_λ^-^_∞_0,Ø^_0,Ø/_
       ≤C/λ(,)_0,Ø.


Hence, the proof follows by combining (<ref>) and (<ref>). 




Now we are in a position to establish the following result.

Let μ_∞>0 be an eigenvalue of _∞ of multiplicity m. Let D be any disc of the complex plane centered at  μ_∞ and containing no other element of the spectrum of _∞. Then, for λ large enough, D contains exactly m eigenvalues of _λ (repeated according to their respective multiplicities). Consequently, each eigenvalue μ_∞>0 of _∞ is a limit of eigenvalues μ of _λ, as λ goes to infinity.
















9


Antonietti
P. F. Antonietti, A. Buffa, and I. Perugia,
Discontinuous Galerkin approximation of the Laplace
              eigenproblem, 
Comput. Methods Appl. Mech. Engrg., 195 (2006), pp. 3483–3503.


















BO
I. Babuška and J. Osborn,
Eigenvalue Problems, 
in Handbook of Numerical Analysis, Vol. II, 
P. G. Ciarlet and J. L. Lions, eds., 
North-Holland, Amsterdam, 1991, pp. 641–787.

ActaBoffi
D. Boffi, 
Finite element approximation of eigenvalue problems,
Acta Numerica, 19, (2010), pp. 1–120.

BoffiBrezziFortinBook
D. Boffi, F. Brezzi, and M. Fortin,
Mixed Finite Element Methods and Applications.
Springer Series in Computational Mathematics, 44. Springer, Heidelberg, 2013.

BoffiBrezziFortin
D. Boffi, F. Brezzi, and M. Fortin,
Reduced symmetry elements in linear elasticity,
Comm. Pure Appl. Anal., 8 (2009), pp. 1–28.

BoffiBrezziGastaldi(a)
D. Boffi, F. Brezzi, and L. Gastaldi, 
On the problem of spurious eigenvalues in the approximation of
linear elliptic problems in mixed form, 
Math. Comp., 69 (2000), pp. 121–140.


BDM
F. Brezzi, J. Douglas, Jr., and L.D. Marini,
Two families of mixed finite elements for second order elliptic problems,
Numer. Math., 47 (1985), pp. 217–235.











 
BuffaPerugia
A. Buffa and I. Perugia,
Discontinuous Galerkin approximation of the Maxwell eigenproblem,
SIAM J. Numer. Anal., 44 (2006), pp. 2198–2226.






D 
M. Dauge,
Elliptic Boundary Value Problems on Corner Domains,
Lecture Notes in Mathematics, 1341, Springer, Berlin, 1988.

DNR1
J. Descloux, N. Nassif, and J. Rappaz,
On spectral approximation. Part 1: The problem of convergence,
RAIRO Anal. Numér., 12 (1978), pp. 97–112.

DNR2
J. Descloux, N. Nassif, and J. Rappaz,
On spectral approximation. Part 2: Error estimates for the Galerkin
method,
RAIRO Anal. Numér., 12 (1978), pp. 113–119.






DiPietroErn
D. A. Di Pietro and A. Ern,
Mathematical Aspects of Discontinuous Galerkin Methods,
Springer-Verlag, Berlin Heidelberg, 2012.






fenics
 A. Logg,   K. A. Mardal,  G. N. Wells  et al.
 Automated Solution of Differential Equations by the Finite Element Method, 
 Springer, 2012.




 






girault-raviart
V. Girault and P.A. Raviart,
Finite Element Methods for Navier-Stokes Equations. Theory and Algorithms,
Berlin, Springer, (1986).







 
 grisvard
 P. Grisvard, 
 Probléms aux limites dans les polygones. Mode démploi,
 EDF, Bull. Dir. Etudes Rech. Ser. C, 1 (1986), pp. 21–59.






 

 
 




 
 MMT
 A. Márquez, S. Meddahi, and T. Tran,
 Analyses of mixed continuous and discontinuous Galerkin methods for the time harmonic 
 elasticity problem with reduced symmetry,
 SIAM J. Sci. Comput., 37 (2015) pp. 1909–1933.
 
MMR
S. Meddahi, D. Mora, and R. Rodríguez,
Finite element spectral analysis for the mixed formulation of
the elasticity equations,
SIAM J. Numer. Anal., 51 (2013) pp. 1041–1063. 
 
MMR-stokes
 S. Meddahi, D. Mora, and R. Rodríguez,
A finite element analysis of a pseudostress formulation for
the Stokes eigenvalue problem,
IMA J. Numer. Anal., 35 (2015) pp. 749–766.
 
 




































