LEAST-SQUARES FORMULATIONS FOR EIGENVALUE PROBLEMS ASSOCIATED WITH LINEAR ELASTICITY

We study the approximation of the spectrum of least-squares operators arising from linear elasticity. We consider a two-field (stress/displacement) and a three-field (stress/displacement/vorticity) formulation; other formulations might be analyzed with similar techniques. We prove a priori estimates and we confirm the theoretical results with simple two-dimensional numerical experiments.


Introduction
In this paper we continue the investigations started in [4] about the spectral properties of operators associated with finite element least-squares formulations. While [4] deals with the Poisson eigenvalue problem, here we consider linear elasticity in the stress/displacement formulation. We discuss a two-field least-squares formulation introduced in [9] and a three-field least-squares formulation studied in [6].
We show that under natural approximation properties of the involved finite element spaces the discretization of the eigenvalue problems converges optimally by using the standard Babuška-Osborn theory [3,7].
An important difference with respect to the case of the Poisson equation, is that for elasticity problem the resulting formulation is not symmetric. We will see that in most cases the eigenvalues are nevertheless real, but there are situations in which complex solutions are found.
As for the analysis presented in [4], it should be clear that the aim of this study is not to present a competitive method for the computation of elasticity eigenmodes. However, knowing the behavior of the spectrum of least-squares operators can be useful for several reasons; for instance, people interested in implementing leastsquares approximations for transient problems, could benefit from such information.
Some numerical tests confirm our theoretical results.

Problem setting
Given a polytopal domain Ω ∈ R d (d = 2, 3) with boundary ∂Ω divided into two parts Γ D and Γ N , we consider the stress/displacement formulation of linear elasticity as a system of first order as follows: find a symmetric d-by-d stress tensor σ and a displacement vectorfield u such that where the compliance tensor is defined in terms of the Lamé constants µ and λ Aτ = 1 2µ τ − λ 2µ + dλ tr(τ )I , the symmetric gradient is defined as and n denotes the outward unit normal vector to Γ N . The formulation we are discussing has the important property to be robust in the incompressible limit. Other formulations could be studied as well. We refer the interested reader, for instance, to the introduction of [6] for an historical perspective.
A least-squares formulation of (1) was considered in [9] by minimizing the functional and X N the subset of X corresponding to the boundary condition τ n = 0 on Γ N .
Remark 1. This formulation does not require a symmetry constraint in the definition of X N since the constitutive equation implies automatically that the solution satisfies σ = σ . We shall refer to this formulation as the two-field formulation.
Other two formulations have been introduced in [6] which make use of three and four fields, respectively, by the introduction of the vorticity and the pressure as new unknowns. We describe the three-field formulation, similar considerations could be derived for the four-field formulation.
The three-field formulation seeks a minimizer of the functional and where as τ = τ − τ /2 is the skew-symmetric part of τ .
We are interested in the eigenvalue problem corresponding to (1), that is, we are looking for ω ∈ R such that for non-vanishing u and for some σ it holds Thanks to the regularity properties of the solution of (1), the eigenvalue problem (2) is compact so that its eigenvalues form an increasing sequence 0 < ω 1 ≤ ω 2 ≤ · · · ≤ ω i ≤ · · · and the eigenspaces are finite dimensional. As usual, we repeat the eigenvalues according to their multiplicity, so that each eigenvalue ω i corresponds to a onedimensional eigenspace E i .
In order to approximate the eigenmodes, we are going to generalize the ideas of [4] to the two-and three-field formulations presented above. In particular, we are not writing directly a least-squares formulation of the eigenvalue problem (which would lead to a non-linear problem), but we study the spectrum of the operators associated with the least-squares source formulations.
2.1. Eigenvalue problem associated with the two-field formulation. The minimization of the functional F(τ , v, f ) gives rise to the following variational formulation: find σ ∈ X N and u ∈ H 1 0,D (Ω) d such that The corresponding eigenvalue problems is obtained after replacing f with ωu, that is: find ω ∈ C such that for a non-vanishing u ∈ H 1 0,D (Ω) d and for some σ ∈ X N we have The structure of the eigenvalue problem (4) in terms of operators is similar to the one described in [4] in the case of the FOSLS formulation for the Poisson equation. Namely, by natural definition of the operators, we are led to the following form: One important difference between this formulation and the one presented in [4] for the Poisson equation is that in our case the operators B and −D are not the same. It follows that it is not possible to show in general that (5) corresponds to a symmetric problem. This fact has important consequences for the numerical approximation. We will see in Section 5 that most of the computed eigenvalues are real, but there are exceptions.

2.2.
Eigenvalue problem associated with the three-field formulation. The variational formulation associated with the minimization of the functional G is obtained by seeking σ ∈ X N , u ∈ H 1 0,D (Ω) d , and ψ ∈L 2 (Ω) such that Also in this case we consider the eigenvalue problem by replacing f with ωu in the right hand side. The problem reads: find ω ∈ C such that for a non-vanishing u ∈ H 1 0,D (Ω) d and for some σ ∈ X N and ψ ∈L 2 (Ω) it holds The operator form of the eigenvalue problem (7) involves 3-by-3 block operators as follows: Also in this case the system is not symmetric and the numerical approximation presented in Section 5 will show that some computed eigenvalues may be complex.

Numerical approximation
As it is apparent from formulations (4) and (7), the numerical approximation will lead to generalized eigenvalue problems of the form Ax = ωBx where the matrix B is singular. From the algebraic point of view, as we observed in [4], the solution of this problem satisfies the following properties.
(1) If the matrix B is invertible then our system is equivalent to the standard eigenvalue problem B −1 Ax = ωx. (2) If K = ker A ∩ ker B is not trivial then the eigenvalue problem is degenerate and vectors in K do not correspond to any eigenvalue of our original system. (3) If the matrix B has a non-trivial kernel which does not contain any nonzero vector of ker(A) then it is conventionally assumed that our system has an eigenvalue ω = ∞ with eigenspace equal to ker(B). (4) If B is singular and A is not (which is the most common situation in our framework) then it may be convenient to switch the roles of the two matrices and to consider the problem Then (γ, x) with γ = 0 corresponds to the eigenmode (∞, x) of the original system; the remaining eigenmodes are (ω, x) with ω = 1/γ.
In the examples we will discuss, we are not going to deal with Case (3) although that situation would open intriguing and interesting new scenarios, similar to what was presented, for instance, in [10]. More aspects of the numerical implementation will be mentioned in Section 5.
3.1. Approximation of the two-field formulation. Given finite dimensional subspaces Σ h ⊂ X N and U h ⊂ H 1 0,D (Ω) d , the Galerkin approximation of (4) reads: find ω h ∈ C such that for a non-vanishing u h ∈ U h and for some σ h ∈ Σ h we have The structure of this problem is analogous of (5) with the natural definition of the involved matrices. The following proposition is the analogue of [4, Prop. 6] and characterizes the number of significant eigenvalues of (9). Proposition 1. The solution of the generalized eigenvalue problem associated with the formulation (9) includes ω = ∞ with multiplicity equal to dim(Σ h )+dim(ker(D)) and a number of complex eigenvalues (counted with their multiplicity) equal to rank(D).
Remark 2. Since the eigenvalue problem stated in (9) considers eigenfunctions with u h = 0, the total number of eigenvalues of the problem we are interested in, is equal to dim U h ; the corresponding values are the rank(D) complex solutions and possibly ω = ∞ if D is not full rank.

3.2.
Approximation of the three-field formulation. Given finite dimensional subspaces Σ h ⊂ X N , U h ⊂ H 1 0,D (Ω) d , and Φ h ⊂L 2 (Ω), the Galerkin approximation of (7) is: find ω h ∈ C such that for a non-vanishing u h ∈ U h and for some σ h ∈ Σ h and ψ h ∈ Φ h we have The matrix structure of this problem corresponds to the one of (7) and the following proposition, in analogy of Proposition 1, gives the characterization of the discrete eigenvalues.

A priori error estimates
It is quite standard to obtain a priori error estimates for eigenvalue problems by studying the convergence of the discrete solution operator towards the continuous one [7]. In our framework, this boils down to showing L 2 (Ω)-estimates for the discretization of (3) and (6), respectively, when the right hand side f is in L 2 (Ω) d . For brevity, we omit the details and we refer the interested reader to [4]. 4.1. L 2 (Ω)-estimate for the two-field formulation. The L 2 (Ω)-estimate for the two-field formulation (3) is the natural generalization of what has been presented in [4] in the case of the Poisson problem. The original idea comes from [2, Sec. 7] (see also [8]).
Let (σ, u) ∈ X N × H 1 0,D (Ω) 2 solve (3) and (σ h , u h ) ∈ Σ h × U h solve the corresponding discrete problem; we consider the dual problem: find s ∈ X N and p ∈ H 1 0,D (Ω) d such that for all τ h ∈ Σ h and v h ∈ U h , where we used the error equations corresponding to (3). It follows We observe explicitly that we cannot obtain any rate of convergence out of the term σ − σ h div since div σ is only in L 2 (Ω) if f is not more regular. On the other hand, the required uniform convergence comes from the (minimal) regularity of the dual problem (11) so that we get u − u h 2 0 ≤ Ch s u − u h 0 f 0 for some positive s as long as the finite element spaces Σ h and U h satisfy the following approximation properties Remark 4. The power s appearing in the estimate of u − u h 0 , which is limited by f ∈ L 2 (Ω) d , is not related to the rate of convergence of the numerical scheme, but guarantees the uniform convergence that implies the correct approximation of the spectrum. The optimal convergence of the scheme is shown in the next theorem by using the full regularity of the eigenspaces.

Theorem 3.
Let ω i = ω i+1 = · · · = ω i+m−1 be an eigenvalue of (4) and denote by E = i+m−1 j=i E j the corresponding eigenspace. If the discrete spaces Σ h and U h satisfy the approximation properties (12) then for h small enough (so that dim U h ≥ i+m−1) the m discrete eigenvalues ω i,h = ω i+1,h = · · · = ω i+m−1,h of (9) converge to ω i ; let us denote by E h the sum of the discrete eigenspaces. Then the following error estimates hold true Proof. The proof is standard (see [4,7]); the only delicate part consists in showing the double order of convergence for the eigenvalues, since the discrete problem is not symmetric. The result can be obtained by considering the adjoint problem, performing the corresponding analysis for its approximation (with ρ(h) replaced by ρ * (h)) and by using the standard Babuška-Osborn theory [3] from which we can conclude that the error of the eigenvalues is bounded by ρ(h)ρ * (h) (note that the continuous eigenvalues have ascent multiplicity equal to one since the continuous problem is symmetric).
4.2. L 2 (Ω)-estimate for the three-field formulation. The L 2 (Ω) error estimate for the discretization of the three-field formulation (7) has been presented in [6,Theorem 3] in the case of a convex domain so that the H 2 regularity of a generalized Stokes problem could be used (see [8]). Since we are not interested in optimal estimates, but only in the uniform convergence in the spirit or Remark 4, the arguments of [8] and [6] could be adapted in order to get the estimate where u solves (6) and u h is the corresponding discrete solution, provided the following approximation properties are satisfied This allows to state the following convergence theorem which is the analogous of Theorem 3 in this situation. In analogy to (12) we make the requirements on the approximation properties of our spaces explicit.

Theorem 4.
Let ω i = ω i+1 = · · · = ω i+m−1 be an eigenvalue of (7) and denote by E = i+m−1 j=i E j the corresponding eigenspace. If the discrete spaces Σ h , U h , and Φ h satisfy the approximation properties (13) then for h small enough (so that dim U h ≥ i + m − 1) the m discrete eigenvalues ω i,h = ω i+1,h = · · · = ω i+m−1,h of (10) converge to ω i ; let us denote by E h the sum of the discrete eigenspaces. Then the following error estimates hold true

Numerical results
Our numerical results confirm the theoretical investigations of the previous sections.
The numerical solution of a generalized eigenvalue problem such as those arising from the discretization of (4) and (7) can be challenging and in this paper we do not focus on the efficiency of the solver but rather on the accuracy of the obtained results.
More specifically, we have to solve a nonsymmetric generalized eigenvalue problem with several infinite eigenvalues. In our computations, we followed two main strategies and we compared the two in order to make sure that no artifact was introduced by the numerical method. We assembled the matrices in FEniCS [1]. Then our first approach is to solve the eigenvalue problem with SLEPc [12] by reversing the matrix on the left-hand side and the one on the right-hand side (see Case (4) at the beginning of Section 3); we used as options "shift-and-invert" with target "magnitude". The second approach consists in exporting the matrices to Matlab and solve with the built-in command "eigs".
More advanced numerical experiments are in progress, which will assess other solvers and compare their performances.

5.1.
Numerical results on the square. We start by taking the domain equal to the unit square Ω =]0, 1[ 2 . In this case a pretty accurate estimate of the first eigenvalue is given by ω = 52.344691168 if we consider which corresponds to a Stokes problem, and homogeneous Dirichlet boundary conditions everywhere, that is Γ N = ∅ (see, for instance, [11]). Since we also want to investigate the symmetry of the formulation, we consider three different mesh sequences: a symmetric and uniform mesh (CROSSED), a non-symmetric and uniform mesh (RIGHT), and a non structured non symmetric mesh (NONUNIF). The meshes are indexed with respect to the number N of subdivisions of the square's sides. The three meshes for N = 4 are plotted in Figure 1.  Table 3. Eigenvalues of the three-field formulation on the unit square: one more refinement for the non uniform mesh In Table 1 we report the numerical results corresponding to computations performed with the two-field scheme (9) when using the three mesh sequences with the level of refinement varying from N = 4 to N = 12. We considered a second order scheme based on Raviart-Thomas elements for the approximation of Σ h (satisfying assumption (12)). The estimate of Theorem 3 predicts fourth order of convergence for the eigenvalues which is confirmed by our tests. It is interesting to observe that the fist computed eigenvalue is always real in our tests.
We performed the same computations for the three-field formulation presented in (10). Also in this case we use a second order scheme based on Raviart-Thomas elements (satisfying assumption (13)). The corresponding results are reported in Table 2. It turns out that also in this case the first computed eigenvalue is always real and that the convergence properties match the theoretical results with the exception of the convergence rate for N = 12 on the non-uniform mesh. In order to check better this phenomenon, we computed one more refinement which is reported in Table 3. It is clear that the overall convergence matches the expected fourth order rate and that the non uniform behavior is due to the fact that the mesh sequence is not structured.
Since the eigenvalue problems corresponding to the considered formulations are not symmetric, it is interesting to investigate whether the computed eigenvalues are complex or real. From the results that we are going show, it is clear that in some cases the computed eigenvalues have a non-vanishing imaginary part; this implies that in general the discrete formulations cannot be reduced to symmetric problems. On the other hand, in most cases the computed eigenvalues are real; this occurs, in particular on the symmetric meshes.
In Table 4 we report the first five eigenvalues computed with the three-field scheme on the non-uniform mesh for N = 10 (similar results could be presented   for the two-field scheme as well). It is apparent that the second and the third eigenvalue are complex conjugate. For the sake of comparison, Table 5 shows the same results on the mesh for N = 11 and we can see that in this case all the first five eigenvalues are real. It is also interesting to observe that the second and the third eigenvalues on the mesh for N = 10 have the same real part, while the corresponding eigenvalues on the mesh for N = 11 are real and different from each other. Clearly they are an approximation of the double eigenvalue ω 2 = ω 3 . In any case, this is perfectly compatible with the statements of Theorems 3 and 4 where the difference |ω j −ω j,h | has to be understood in the sense of the distance in the complex plane.

5.2.
Numerical results on the L-shaped domain. We conclude our numerical tests with computations on the L-shaped domain where the re-entrant corner causes singularities in the solution. We consider a uniform mesh (see Figure 2 for the case N = 4).
An estimate of the first eigenvalues in this case has been provided in [5] and corresponds to ω = 32.13269464746. Tables 6 and 7 report on the computed approximation with the two-and three-field approximations (9) and (10), respectively. As expected, the singularity of the eigenfunction corresponding to the first eigenvalue causes a reduction of the order of convergence.  Also in this case the first computed eigenvalue is real, however, also in presence of a symmetric mesh, it turns out that there might be complex eigenvalues. For instance, Table 8 reports some higher eigenvalues computed for N = 8. Also in this case the complex eigenvalues converge towards real number according to the statement of Theorems 3 and 4; moreover, the presence of complex eigenvalues depends on the mesh and on the level of refinements. The effect of the mesh on complex eigenvalues will be the object of further investigation.