    # Least-squares for linear elasticity eigenvalue problem

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.

## Authors

02/19/2020

### First order least-squares formulations for eigenvalue problems

In this paper we discuss spectral properties of operators associated wit...
03/23/2022

### Spectral analysis of a mixed method for linear elasticity

The purpose of this paper is to analyze a mixed method for linear elasti...
07/13/2020

### Basis functions for residual stresses

We consider arbitrary preexisting residual stress states in arbitrarily ...
08/09/2021

### On the spectrum of an operator associated with least-squares finite elements for linear elasticity

In this paper we provide some more details on the numerical analysis and...
08/18/2021

### On the spectrum of the finite element approximation of a three field formulation for linear elasticity

We continue the investigation on the spectrum of operators arising from ...
12/11/2020

### DPG approximation of eigenvalue problems

In this paper, the discontinuous Petrov–Galerkin approximation of the La...
04/22/2020

### Approximation of the Axisymmetric Elasticity Equations with Weak Symmetry

In this article we consider the linear elasticity problem in an axisymme...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

In this paper we continue the investigations started in  about the spectral properties of operators associated with finite element least-squares formulations. While  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  and a three-field least-squares formulation studied in .

We show that under natural approximation properties of the involved finite element spaces the discretization of the eigenvalues problems converges optimally by using the standard Babuška–Osborn theory [2, 6].

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 , 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 least-squares approximation for transient problems, could benefit from such information.

Some numerical tests confirm our theoretical results.

## 2. Problem setting

Given a polytopal domain () with boundary divided into two parts and , we consider the stress/displacement formulation of linear elasticity as a system of first order as follows: find a symmetric -by-

stress tensor

and a displacement vectorfield such that

 (1) Aσ––−ε––(u)=0 in Ω divσ––=−f in Ω u=0 on ΓD σ––n=0 on ΓN,

where the compliance tensor is defined in terms of the Lamé constants and

 Aτ––=12μ(τ––−λ2μ+dλtr(τ––)I–),

the symmetric gradient is defined as

 ε––(v)=12(∇–––v+(∇–––v)⊤),

and

denotes the outward unit normal vector to

.

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  for an historical perspective.

A least-squares formulation of (1) was considered in  by minimizing the functional

 F(τ––,v;f)=∥Aσ––−ε––(u)∥20+∥divσ––+f∥20

in , with

 X––={H(div;Ω)dif ΓN≠∅{τ––∈H(div;Ω)d:∫Ωtr(τ––)dx=0}if ΓN=∅

and the subset of corresponding to the boundary condition on .

###### Remark 1.

This formulation does not require a symmetry constraint in the definition of 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  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

 G(τ––,v,φ;f)=∥Aτ––−∇–––v+(−1)dχ––φ∥20+∥divσ––+f∥20+∥asτ––∥20

in with

 ¯L2(Ω)={L2(Ω)if ΓN≠∅{φ∈L2(Ω):∫Ωφdx=0}if ΓN=∅,

where

 χ––=⎧⎪⎨⎪⎩(0−110)if d=2(χ––1,χ––2,χ––3)if d=3

with

 χ––1=⎛⎜⎝00000−1010⎞⎟⎠χ––2=⎛⎜⎝001000−100⎞⎟⎠χ––3=⎛⎜⎝0−10100000⎞⎟⎠

and where

is the skew-symmetric part of

.

We are interested in the eigenvalue problem corresponding to (1), that is, we are looking for such that for non vanishing and for some it holds

 (2) Aσ––−ε––(u)=0 in Ω divσ––=−ωu in Ω u=0 on ΓD σ––n=0 on ΓN.

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

corresponds to a one-dimensional eigenspace .

In order to approximate the eigenmodes, we are going to generalize the ideas of  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 gives rise to the following variational formulation: find and such that

 (3) {(Aσ––,Aτ––)+(divσ––,divτ––)−(Aτ––,ε––(u))=−(f,divτ––)∀τ––∈X––N−(Aσ––,ε––(v)+(ε––(u),ε––(v))=0∀v∈H10,D(Ω)d.

The corresponding eigenvalue problems is obtained after replacing with , that is: find such that for a non vanishing and for some we have

 (4) {(Aσ––,Aτ––)+(divσ––,divτ––)−(Aτ––,ε––(u))=−ω(u,divτ––)∀τ––∈X––N−(Aσ––,ε––(v)+(ε––(u),ε––(v))=0∀v∈H10,D(Ω)d.

The structure of the eigenvalue problem (4) in terms of operators is similar to the one described in  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:

 (5) (AB⊤BC)(xy)=ω(0D00)(xy).

One important difference between this formulation and the one presented in  for the Poisson equation is that in our case the operators and 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 is obtained by seeking , , and such that

 (6) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(Aσ––,Aτ––)+(divσ––,divτ––)+(as(σ––),as(τ––))−(Aτ––,ε––(u))+(−1)d(Aτ––,χ––ψ)=−(f,divτ––)∀τ––∈X––N−(Aσ––,ε––(v)+(ε––(u),ε––(v))−(−1)d(χ––ψ,∇–––v)=0∀v∈H10,D(Ω)d(−1)d(Aσ––,χ––φ)−(−1)d(χ––φ,∇–––u)+(χ––ψ,χ––φ)=0∀φ∈¯L2(Ω).

Also in this case we consider the eigenvalue problem by replacing with in the right hand side. The problems reads: find such that for a non vanishing and for some and it holds

 (7) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(Aσ––,Aτ––)+(divσ––,divτ––)+(as(σ––),as(τ––))−(Aτ––,ε––(u))+(−1)d(Aτ––,χ––ψ)=−ω(u,divτ––)∀τ––∈X––N−(Aσ––,ε––(v)+(ε––(u),ε––(v))−(−1)d(χ––ψ,∇–––v)=0∀v∈H10,D(Ω)d(−1)d(Aσ––,χ––φ)−(−1)d(χ––φ,∇–––u)+(χ––ψ,χ––φ)=0∀φ∈¯L2(Ω).

The operator form of the eigenvalue problem (7) involves -by- block operators as follows:

 (8)

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.

## 3. Numerical approximation

As it is apparent from formulations (4) and (7), the numerical approximation will lead to generalized eigenvalue problems of the form where the matrix is singular. From the algebraic point of view, as we observed in , the solution of this problem satisfies the following properties.

1. If the matrix is invertible then our system is equivalent to the standard eigenvalue problem .

2. If is not trivial then the eigenvalue problem is degenerate and vectors in do not correspond to any eigenvalue of our original system.

3. If the matrix has a non-trivial kernel which does not contain any nonzero vector of then it is conventionally assumed that our system has an eigenvalue with eigenspace equal to .

4. If is singular and 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

 Bx=γAx

Then with corresponds to the eigenmode of the original system; the remaining eigenmodes are with .

### 3.1. Approximation of the two-field formulation

Given finite dimensional subspaces and , the Galerkin approximation of (4) reads: find such that for a non vanishing and for some we have

 (9) {(Aσ––h,Aτ––)+(divσ––h,divτ––)−(Aτ––,ε––(uh))=−ωh(uh,divτ––)∀τ––∈Σh−(Aσ––h,ε––(v)+(ε––(uh),ε––(v))=0∀v∈Uh.

The structure of this problem is analogous of (5) with the natural definition of the involved matrices. The following proposition is the analogue of [3, 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) include with multiplicity equal to and a number of complex eigenvalues (counted with their multiplicity) equal to .

###### Remark 2.

Since the eigenvalue problem stated in (9

) considers eigenfunctions with

, the total number of eigenvalues of the problem we are interested in, is equal to ; the corresponding values are the complex solutions and possibly if is not full rank.

### 3.2. Approximation of the three-field formulation

Given finite dimensional subspaces , , and , the Galerkin approximation of (7) is: find such that for a non vanishing and for some and we have

 (10) ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩(Aσ––h,Aτ––)+(divσ––h,divτ––)+(as(σ––h),as(τ––))−(Aτ––,ε––(uh))+(−1)d(Aτ––,χ––ψh)=−ωh(uh,divτ––)∀τ––∈Σh−(Aσ––h,ε––(v)+(ε––(uh),ε––(v))−(−1)d(χ––ψh,∇–––v)=0∀v∈Uh(−1)d(Aσ––h,χ––φ)−(−1)d(χ––φ,∇–––uh)+(χ––ψh,χ––φ)=0∀φ∈Φh.

The matrix structure of this problem corresponds to (7) and the following proposition, in analogy of Proposition 1, gives the characterization of the discrete eigenvalues.

###### Proposition 2.

The solution of the generalized eigenvalue problem associated with the formulation (10) consists of and a number of complex eigenvalues (counted with their multiplicity) equal to .

###### Remark 3.

As in Remark 2, the eigenvalues of formulation (10) are complex values and possibly if is not full rank.

## 4. 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 . In our framework, this boils down to showing estimates for the discretization of (3) and (6), respectively, when the right hand side is in . For brevity, we omit the details and we refer the interested reader to .

### 4.1. L2(Ω) estimate for the two-field formulation

The estimate for the two-field formulation (3) is the natural generalization of what has been presented in  in the case of the Poisson problem. The original idea comes from [1, Sec. 7] (see also ).

Let solve (3) and solve the corresponding discrete problem; we consider the dual problem: find and such that

 (11) {(As–,At–)+(divs–,divt–)−(At–,ε––(p))=0∀t–∈X––N−(As–,ε––(q)+(ε––(p),ε––(q))=(u−uh,q)∀q∈H10,D(Ω)d.

Taking as test functions and we have

 ∥u−uh∥20 =(As–,A(σ––−σ––h))+(divs–,div(σ––−σ––h)) −(ε––(p),A(σ––−σ––h))−(As–,ε––(u−uh)) +(ε––(p),ε––(u−uh)) =(A(s–−τ––h),A(σ––−σ––h))+(div(s–−τ––h),div(σ––−σ––h)) −(ε––(p−vh),A(σ––−σ––h))−(A(s–−τ––h),ε––(u−uh)) +(ε––(p−vh),ε––(u−uh))

for all and , where we used the error equations corresponding to (3). It follows

 ∥u−uh∥20 ≤C(∥s–−τ––h∥div+∥p−vh∥1)(∥σ––−σ––h∥div+∥u−uh∥1) ≤C(∥s–−τ––h∥div+∥p−vh∥1)∥f∥0.

We observe explicitly that we cannot obtain any rate of convergence out of the term since is only in if 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−uh∥20≤Chs∥u−uh∥0∥f∥0

for some positive as long as the finite element spaces and satisfies the following approximation properties

 (12) infτ––∈Σh∥s–−τ––∥div≤Chs∥s–∥s+∥divs–∥s infv∈Uh∥p−vh∥1≤Chs∥p∥1+s.
###### Remark 4.

The power appearing in the estimate of , which is limited by , 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 be an eigenvalue of (4) and denote by the corresponding eigenspace. If the discrete spaces and satisfy the approximation properties (12) then for small enough (so that ) the discrete eigenvalues of (9) converge to ; let us denote by the sum of the discrete eigenspaces. Then the following error estimates hold true

 δ(E,Eh)≤Cρ(h) |ωj−ωj,h|≤ρ(h)2 j=i,…,i+m−1

with

 ρ(h)=supu∈E∥u∥=1infτ––∈Σhv∈Uh(∥ε–(u)−τ––∥div+∥u−v∥1).
###### Proof.

The proof is standard (see [3, 6]); the only delicate part consists in showing the double order or 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 replaced by ) and by using the standard Babuška–Osborn theory  from which we can conclude that the error of the eigenvalues is bounded by (note that the continuous eigenvalues have ascent multiplicity equal to one since the continuous problem is symmetric). ∎

### 4.2. L2(Ω)-estimate for the three-field formulation

The error estimate for the discretization of the three-field formulation (7) has been presented in [5, Theorem 3] in the case of a convex domain so that the regularity of a generalized Stokes problem could be used (see ). Since we are not interested in optimal estimates, but only in the uniform convergence in the spirit or Remark 4, the arguments of  and  could be adapted in order to get the estimate

 ∥u−uh∥0≤Chs∥f∥0

where solves (6) and is the corresponding discrete solution, provided the following approximation properties are satisfied

 (13) infτ––∈Σh∥s–−τ––∥div≤Chs∥s–∥s+∥divs–∥s infv∈Uh∥p−vh∥1≤Chs∥p∥1+s infφ∈Φh∥ψ−φ∥0≤Chs∥ψ∥s.

This allows to state the following convergence theorem which is the analogous of Theorem 3 in this situation. In analogy to (12) we make explicit the requirements on the approximation properties of our spaces.

###### Theorem 4.

Let be an eigenvalue of (7) and denote by the corresponding eigenspace. If the discrete spaces , , and satisfy the approximation properties (13) then for small enough (so that ) the discrete eigenvalues of (10) converge to ; let us denote by the sum of the discrete eigenspaces. Then the following error estimates hold true

 δ(E,Eh)≤Cρ(h) |ωj−ωj,h|≤ρ(h)2 j=i,…,i+m−1

with

 ρ(h)=supu∈E∥u∥=1infτ––∈Σhv∈Uhφ∈Φh(∥–∇(u)−τ––∥div+∥u−v∥1+∥curl(u)−φ∥0).

## 5. 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.

### 5.1. Numerical results on the square

We start by taking the domain equal to the unit square . In this case a pretty accurate estimate of the first eigenvalue is given by if we consider

 Aτ––=12(τ––−12tr(τ––)I–),

which corresponds to a Stokes problem, and homogeneous Dirichlet boundary conditions everywhere, that is (see, for instance, ). 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 of subdivisions of the square’s sides. The three meshes for are plotted in Figure 1.

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 to . We considered a second order scheme based on Raviart–Thomas elements for the approximation of (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 with the theoretical results with the exception of the convergence rate for 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 with 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 (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 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 have the same real part, while the corresponding eigenvalues on the mesh for are real and different from each other. Clearly they are approximation of the double eigenvalue . In any case, this is perfectly compatible with the statements of Theorems 3 and 4 where the difference 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 ).

An estimate of the first eigenvalues in this case has been provided in  and corresponds to . 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 . 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.

## References

•  Douglas N. Arnold, Daniele Boffi, and Richard S. Falk, Quadrilateral finite elements, SIAM J. Numer. Anal. 42 (2005), no. 6, 2429–2451. MR 2139400
•  Ivo Babuška and John Osborn, Eigenvalue problems, Handbook of numerical analysis, Vol. II, Handb. Numer. Anal., II, North-Holland, Amsterdam, 1991, pp. 641–787.
•  Fleurianne Bertrand and Daniele Boffi, First order least-squares formulations for eigenvalue problems, arXiv:2002.08145 [math.NA], 2020.
•  Fleurianne Bertrand, Daniele Boffi, and Rui Ma, An adaptive finite element scheme for the Hellinger–Reissner elasticity mixed eigenvalue problem, submitted, 2020.
•  Fleurianne Bertrand, Zhiqiang Cai, and Eun Young Park, Least-squares methods for elasticity and Stokes equations with weakly imposed symmetry, Comput. Methods Appl. Math. 19 (2019), no. 3, 415–430. MR 3977481
•  Daniele Boffi, Finite element approximation of eigenvalue problems, Acta Numer. 19 (2010), 1–120. MR 2652780
•  Zhiqiang Cai and Jaeun Ku, The norm error estimates for the div least-squares method, SIAM J. Numer. Anal. 44 (2006), no. 4, 1721–1734. MR 2257124
•  Zhiqiang Cai and Gerhard Starke, Least-squares methods for linear elasticity, SIAM J. Numer. Anal. 42 (2004), no. 2, 826–842. MR 2084237
•  Joscha Gedicke and Arbaz Khan, Arnold-Winther mixed finite elements for Stokes eigenvalue problems, SIAM J. Sci. Comput. 40 (2018), no. 5, A3449–A3469. MR 3864690