Adaptive Morley element algorithms for the biharmonic eigenvalue problem

This paper is devoted to the adaptive Morley element algorithms for a biharmonic eigenvalue problem in Rn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb{R}^{n}$\end{document} (n≥2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$n\geq2$\end{document}). We combine the Morley element method with the shifted-inverse iteration including Rayleigh quotient iteration and the inverse iteration with fixed shift to propose multigrid discretization schemes in an adaptive fashion. We establish an inequality on Rayleigh quotient and use it to prove the efficiency of the adaptive algorithms. Numerical experiments show that these algorithms are efficient and can get the optimal convergence rate.


Introduction
Biharmonic equation/eigenvalue problem plays an important role in elastic mechanics. In 1968, Morley designed a famous non-conforming element called the Morley element [1] to solve biharmonic equation (plate bending problem). The Morley element was extended to arbitrarily dimensions by Wang and Xu [2] in 2006. For biharmonic equation, the a priori/a posteriori error estimate was studied in [3][4][5][6] and the convergence and optimality of the adaptive Morley element method was proved in [7,8]. The Morley element has been employed to solve the biharmonic eigenvalue problem, including the vibration of a plate; and [9] studied its a priori error estimate. [10,11] studied a posteriori error estimate and the adaptive method, [12] adopted a new method dispensing with any additional regularity assumption to study the error estimates and adaptive algorithms. This paper further studies the adaptive Morley element method and has the following features: 1. The adaptive finite element methods, which were first proposed by Babuska and Rheinboldt [13], have gained an extensive attention in academia. More and more researchers entered this field and obtained many good results, most of which have been systemically summarized in [5,[14][15][16]. And [10,12] have employed the adaptive Morley element algorithms for the biharmonic eigenvalue problem based on solving directly the original eigenvalue problem a(u, v) = λb (u, v) in each iteration. In this paper, we establish the adaptive Morley element algorithms based on the shifted-inverse iteration including Rayleigh quotient iteration and the inverse iteration with fixed shift to solve the biharmonic eigenvalue problem. The shifted-inverse iteration method based on the multigrid discretizations has been studied in-depth (see [17] and the references therein), but they did not involve the Morley element. With our method, the solution of an original eigenvalue problem is reduced to the solution of an eigenvalue problem on a much coarser grid and the solution of a series of linear algebraic equations on finer and finer grids. Therefore, our method is more efficient than the method in [10,12]. 2. For fourth order equations in R 3 , it is difficult to employ a conforming element. For instance, Zenicek constructed a conforming tetrahedral finite element with 9 degree of polynomials and 220 nodal parameters [5], while the Morley tetrahedral element [2] has only 10 nodal parameters. Based on [4], we comply with the adaptive Morley element computation for the biharmonic eigenvalue problem in R 3 . Numerical results indicate that the adaptive algorithms are very efficient. 3. A family of good adaptive meshes should satisfy h = O(h α min ), where h is the mesh size, h min is the diameter of the smallest element, and α is the regularity index of the biharmonic equation over the domain with reentrant corner (see [18]). However, we find through the numerical computation that h h α min will become bigger and bigger when the iteration increases for the standard adaptive algorithm. Thus, referring to [19], we combine the standard local refined adaptive algorithm with uniformly refined algorithm to give new algorithms.
It is easy to verify that a(u, v) is a symmetric, continuous, and H 2 0 ( )-elliptic bilinear form. Let u a = √ a(u, u), then the norms u a , u 2 , and |u| 2 are equivalent. We assume that π h = {κ} is a regular simplex partition of and satisfies = κ (see [20]). Let h κ be the diameter of κ, and h = max{h κ : κ ∈ π h } be the mesh size of π h (h < 1), h min = min{h κ : κ ∈ π h }. Let ε h = {F} denote the set of faces ((n -1)-simplexes) of π h , and let ε h = {l} denote the set of faces (n -2)-simplexes of π h . When n = 2, l = z is a vertex of κ, and 1 meas(l) l v = v(z). Let π h (κ) denote the set of all elements sharing a common face with the element κ. Let κ + and κbe any two n-simplexes with a face F in common such that the unit outward normal to κat F corresponds to γ F . We denote the jump of v across the face F by And the jump on boundary faces is simply given by the trace of the function on each face.
In the papers [2,5], the Morley element space is defined by where P 2 (κ) denotes the space of polynomials of degree less than or equal to 2 on κ.
Define the interpolation operator I h : From Lemma 8 in [2], we know that | · | 2,h is equivalent to · 2,h , · 2,h is a norm in S h , and a h (·, ·) is a uniformly S h -elliptic bilinear form, and · h = a h (·, ·) 1 2 is a norm in S h . And the following equality holds for any w ∈ H 2 0 ( ): The discrete form of (2.2) reads: The corresponding boundary value problem of (2.1) is From [18], we know that where α ∈ ( 1 2 , 1) for the domain with reentrant corner, and α = 1 for the convex domain in R 2 .
The weak form of (2.4) and its discrete form are to find w ∈ H 2 0 ( ) such that and to find w h ∈ S h such that Define the solution operators T : L 2 ( ) → H 2 0 ( ) ⊂ L 2 ( ) and T h : L 2 ( ) → S h as follows: Then T, T h : L 2 ( ) → L 2 ( ) are self-adjoint and compact. It is well known that the eigenvalue problem (2.1) has countably many eigenvalues, which are real and positive diverging to +∞. Suppose that λ and λ h are the kth eigenvalue of (2.2) and (2.3), respectively, the algebraic multiplicity of λ is equal to q, λ = λ k = λ k+1 = · · · = λ k+q-1 . Let M(λ) be the space spanned by all eigenfunctions corresponding to λ and M h (λ) be the direct sum of eigenspaces corresponding to all eigenvalues of (2.3) that converge to λ. LetM(λ) = {u : u ∈ M(λ), u h = 1}. Now we introduce the following quantity: The saturation condition was analyzed in [21][22][23], especially, it was analyzed in [22] for very general cases. According to this condition, we can make the following assumption: where C 1 and C 2 are independent of mesh parameters.
Due to the generalized Poincare-Friedrichs inequality, Theorem 3 in [24] and a h (u -I h u, v) = 0, ∀v ∈ S h (see [5]), we deduce for any w ∈ S h , u ∈ H 2 0 ( ) where C 3 is a positive constant independent of mesh parameters. From (2.5) we have the following estimate using the Cauchy-Schwarz inequality: For any g ∈ L 2 ( ), T h g ∈ S h satisfying Suppose w ∈ H 2+r ( ), r ∈ ( 1 2 , 1], then we have the following estimate: Using the trace inequality [5] proves the above estimate under the case r = 1. Using the arguments in [5], we can obtain the above estimate under the case r = ( 1 2 , 1] (also see [11]). We can derive the following Lemma 2.1 from Lemma 2.3 in [25].
where constants C 5 and C 6 are positive and only depend on λ.
The following inequality on Rayleigh quotient plays an important role.
which together with u By Lemma 2.5 in [26], we get Hence, from inequalities (2.7)-(2.9) we deduce We get the results that we need.

The shifted-inverse iteration based on multigrid discretization
Let {S h i } ∞ 0 be a family Morley element spaces, h 0 = H. Refer to the references [17], we present the following calculation schemes. Scheme 1 (Rayleigh quotient iteration based on multigrid discretizations) Given the iteration times l. Step Step 2.
Step 3. Solve a linear system on S h i : Find u ∈ S h i such that Step 4. Compute the Rayleigh quotient: Step 5. If i = l, then output (λ h l , u h l ), stop; else, i ⇐ i + 1, and return to Step 3.
Scheme 2 (The inverse iteration with fixed shift based on multigrid discretizations) Given the iteration times l and i 0 .
Step 6. Solve a linear system on S h i : Find u ∈ S h i such that Step 7. Compute the Rayleigh quotient Step 8. If i = l, then output (λ h l , u h l ), stop; else, i ⇐ i + 1, and return to Step 6.
Strictly speaking, the above a h (·, ·) and · h should be written as a h i (·, ·) and · h i . For the sake of simplicity, we write a h i (·, ·) and · h i as a h (·, ·) and · h , in this paper.

The theoretical analysis
In this section, we will prove the convergence of (λ h l , u h l ) derived from Scheme 1/ Scheme 2, and that the constants appearing in the error estimates are not only independent of mesh parameter but also iterative times l.
In the following discussion, let (λ k , u k ) and (λ k,h , u k,h ) denote the kth eigenpair of (2.2) and (2.3), respectively, and μ k = 1 Our analysis is based on the following Lemma 4.1 (see Lemma 4.1 in [17]).
then the following inequality holds: Next, we will use the proof method in [17] to analyze the error of Schemes 1-2. Let δ 0 be a positive constant satisfying the following inequalities: where λ 0 is an approximate eigenvalue of λ k , u h l is an approximate eigenfunction obtained by Scheme 1 or Scheme 2, and ρ is the separation constant of the eigenvalue μ k = 1 λ k . Condition 4.1 plays a key role in proving Theorem 4.1, by which we can prove Theorems 4.2-4.3. In the proof of Theorems 4.2-4.3, we can deduce that Condition 4.1 holds when the mesh size H is appropriately small. However, it is difficult to verify the condition whether the mesh size H is appropriately small or not. And it seems to be a necessary condition in many papers on the convergence and error estimates of the finite element method for eigenvalue problem. But numerical experiments in Sect. 6 present a satisfying practical performance for our algorithms, which shows that it is unnecessary for the mesh size H to be appropriately small, even though the theory is not complete.
The following Theorems 4.
Proof We use Lemma 4.1 to complete the proof. Select μ 0 = 1 λ 0 and u 0 = Then, by (2.6) and (2.8), we have Using the triangle inequality, (4.7), (2.11), Condition 4.1 and (4.3), we get Hence, the conditions in Lemma 4.1 are verified. By (2.5) we see that Step 3 in Scheme 1 or Step 6 in Scheme 2 is equivalent to the following: Then Step 3 in Scheme 1 or Step 6 in Scheme 2 is equivalent to From (4.4), (2.13) and (4.5), we derive that Let the eigenvectors {u j,h l } k+q-1 k be an orthogonal basis of M h l (λ k ) with respect to a h (·, ·). Denote Hence, substituting (4.8) and (4.9) into (4.1), we obtain Using (2.10), we deduce that Noting that the constants C 3 , C 5 , C 6 , C 7 and ρ are independent of mesh parameters and iterative times l, and u h l-1 k -ū h ≤ δ 0 , |λ 0λ k | ≤ δ 0 and δ h l (λ k ) ≤ δ 0 , by (4.10) and (4.2), we know that there exists a positive constant C 0 that is independent of mesh parameters and l such that (4.6) holds. And we can have C 0 ≥ C 5 .
We need the following two conditions (see 18. For a nonsmooth eigenfunction, the condition could be met when the local refinement is done near the singular point.
Proof The proof is completed by using induction and Theorem 4.1 with λ 0 = λ h i 0 k . Note that δ H (λ k ) → 0 (H → 0), then there is a proper small H 0 > 0 such that if H ≤ H 0 , Lemma 2.1 and the following inequalities hold: When l = i 0 + 1, by Theorem 4.2 we know that there exists u k ∈ M(λ k ) such that Suppose that Theorem 4.3 holds for l -1, i.e. there existsū ∈ M(λ k ) such that Then we infer from (4.17) that the conditions of Theorem 4.1 hold; therefore, for l, we can get which together with (4.18), we get (4.15). Substituting (4.15) into the inequality (2.12), we get (4.16).
Remark For some adaptive local refined grids used usually, (2.9) can be expressed as therefore r in the theorems of this paper can take 1.

Adaptive algorithms
In this section, referring to [10,17,27], we present six algorithms. We denote Algorithm 1 in [10] as Algorithm 1 in this paper, and Algorithms 2-3 are established based on Schemes 1-2, respectively. Then we combine Algorithms 1-3 with a uniformly refined algorithm to get Algorithms 1M-3M, respectively. And the a posterior error estimator in the following algorithms comes from [4], that is where w h is the finite element approximate solution of (2.4), τ F is the tangential vector and γ F the unit outward normal on F ∈ ε h . In the following algorithms, we have to provide an initial shape regular triangulation π h 0 and a parameter θ ∈ (0, 1). Also, from [10,11] we know that replacing w h with u h and replacing f with λ h u h in (5.1), we can obtain the error estimator of Algorithms 1 and 1M. By Lemma 4.1 we can deduce that replacing w h with u h and replacing f with λ h u h in (5.1), we can obtain the error estimator of Algorithms 2-3 and Algorithms 2M-3M.
Step 4. Compute the local indicators η h l (λ h l u h l , u h l , κ).
Step 5. Constructπ h l ∈ π h l by Marking strategy E and θ .
Step 6. Refine π h l to get a new mesh π h l+1 by procedure Refine.
Step 8. l ⇐ l + 1 and go to Step 4.
Step 6. Refine π h l to get a new mesh π h l+1 by procedure Refine.
Step 9. If h l h α l min ≥ C r , then uniformly refine the mesh π h l to get a new mesh π h l+1 and go to Step 7, else go to Step 4.

Algorithm 2M
Choose the parameter 0 < θ < 1, α, and a bound C r of h l h α l min .
Step 9. If h l h α l min ≥ C r , then uniformly refine the mesh π h l to get a new mesh π h l+1 and go to Step 7, else go to Step 4.

Algorithm 3M
Choose the parameter 0 < θ < 1, an integer i 0 , α, and a bound C r of h l h α l min .
Step 9. If h l h α l min ≥ C r , then uniformly refine the mesh π h l to get a new mesh π h l+1 and go to Step 7, else go to Step 4.
Marking strategy E Given parameter 0 < θ < 1: Step 1. Construct a minimal subset π h l of π h l by selecting some elements in π h l such that Step 2. Mark all the elements in π h l .

Marking strategy E1
To get Marking strategy E1 we only replace λ h l and u h l in Marking strategy E with λ h l and u h l , respectively.
Algorithms 1M-3M including steps with uniform refinement seem to be opposite to the adaptive concept. Indeed, the combination of adaptive algorithms and uniform refinement meets the certain mesh-grading properties, thus improving the efficiency of Algorithms 1-3 (see Tables 1-3 in Sect. 6).

Numerical experiment
In this section, we compute the smallest eigenvalue of (2.1) on the L-shaped domain (0, 1) 2 \ Our programs are compiled on MATLAB2012a under the package of Chen [28] using HP-Z230 workstation with ROM 32G and CPU 3.60 GHz.
We use the command "\" to solve (5.2) and use the sparse solver eigs(A, B, 1, sm ) to solve (2.3) for the smallest eigenvalues. Before showing the results, some symbols need to be explained: λ h l the smallest eigenvalue obtained by the lth iteration using Algorithm 1. λ R h l the smallest eigenvalue obtained by the lth iteration using Algorithm 2. λ F h l the smallest eigenvalue obtained by the lth iteration using Algorithm 3. λ M h l the smallest eigenvalue obtained by the lth iteration using Algorithm 1M. λ RM h l the smallest eigenvalue obtained by the lth iteration using Algorithm 2M. λ FM h l the smallest eigenvalue obtained by the lth iteration using Algorithm 3M. N dof the number of the degree of freedom. CPU(s) the time CPU runs from the first iteration to the current iteration. In R 2 , the initial mesh π h 0 is isosceles right triangle subdivision with mesh size √ 2 32 , and we take θ = 0.25, C r = 1.1, α = 1 2 . We fix shift from the 25th and 13th in Algorithm 3 and Algorithm 3M, respectively. The results are shown in Tables 1-3. We depict the error curves of Algorithms 1-3 and Algorithms 1M-3M in Figs. 1-3.
From Tables 1-3, we can get the conclusion that in the case the accurate are almost same, Algorithms 2-3 take about half time of Algorithm 1. In the case the accurate are almost same, Algorithm iM takes about 2 5 time of Algorithm i, i = 1, 2, 3. The smallest eigenvalue of (2.1) is unknown. Therefore, we replace it with an approximate eigenvalue λ 1 ≈ 6703.585 in R 2 with high accuracy. It is present that the relative error curves of the smallest eigenvalues derived from Algorithms 1-3 and Algorithms 1M-3M on the adaptive meshes in Figs. 1-3, whose slopes are more or less -1, which shows that all the six Morley element adaptive algorithms can get the optimal convergence rate O(h 2 ) in R 2 . In R 3 , the initial mesh π h 0 is tetrahedron subdivision with mesh size √ 3 16 , and we take θ = 0.25 and λ 1 ≈ 8290.011 with high accuracy replacing the accurate eigenvalue. It is present that the refined mesh and the relative error curves of the smallest eigenvalues derived from Algorithms 1-2 in Fig. 4, from which we see that Algorithm 2 is more efficient than Algorithm 1, but meanwhile we also see from Table 4 that the mesh size has no change.

Figure 1
The convergence rates of the smallest eigenvalue from Algorithm 1(a) and Algorithm 1M(b)

Figure 2
The convergence rates of the smallest eigenvalue from Algorithm 2(a) and Algorithm 2M(b) Figure 3 The convergence rates of the smallest eigenvalue from Algorithm 3(a) and Algorithm 3M(b) Figure 4 The refined mesh for the L-shaped domain (a) and the convergence rates of the smallest eigenvalue from Algorithm 1 and Algorithm 2(b) in R 3 Because N dof in R 3 increases very fast after uniform refinement, which leads to surpassing computer's memory, we cannot employ Algorithms 1M-3M to solve (2.1).