A Legendre spectral method based on a hybrid format and its error estimation for fourth-order eigenvalue problems

: In this paper, we developed and studied an e ffi cient Legendre spectral method for fourth order eigenvalue problems with the boundary conditions of a simply supported plate. Initially, a new variational formulation based on a hybrid format and its discrete variational form were established. We then employed the spectral theory of complete continuous operators to establish the prior error estimates of the approximate solutions. By integrating approximation results of some orthogonal projection operators in weighted Sobolev spaces, we further gave the error estimation for the approximating eigenvalues and eigenfunctions. In addition, we developed an e ff ective set of basis functions by utilizing the orthogonal properties of Legendre polynomials, and subsequently derived the matrix eigenvalue system of the discrete variational form for both two-dimensional and three-dimensional cases, based on a tensor product. Finally, numerical examples were provided to demonstrate the exponential convergence and e ffi ciency of the algorithm.

For the finite element method applied to fourth-order equations and eigenvalue problems, a C 1 continuous finite element space is typically required.This not only complicates the construction of basis functions but also leads to a significant number of degrees of freedom.The spectral method, which is known to all, is a high-order numerical method with a spectral accuracy and plays a vital role in finding numerical solutions of many differential equations [18][19][20][21][22].However, it is necessary for the computational domain to be square or cubic.To address this limitation, some spectral element methods are commonly employed to solve differential equations on general domains.For the spectral element methods applied to second-order problems, both their theoretical foundation and numerical calculations are well-established.However, for the spectral element methods applied to fourth-order problems on general domains, the construction of basis functions is also intricate and the computational load is substantial.Therefore, it is highly significant to introduce spectral methods based on hybrid format for the fourth-order equations and eigenvalue problems.To our knowledge, there are few reports on Legendre spectral methods based on a hybrid format for fourth-order eigenvalue problems with the boundary conditions of simply supported plates.
Thus, our aim of the current paper is to propose an efficient Legendre spectral method for fourth order eigenvalue problems with the boundary conditions of a simply supported plate.Initially, a new variational formulation based on a hybrid format and its discrete variational form are established.We then employ the spectral theory of complete continuous operators to establish the prior error estimates of the approximate solutions.By integrating approximation results of some orthogonal projection operators in weighted Sobolev spaces, we further give the error estimation for approximating eigenvalues and eigenfunctions.In addition, we developed an effective set of basis functions by utilizing the orthogonal properties of Legendre polynomials, and subsequently derive the matrix eigenvalue system of the discrete variational form for both two-dimensional and three-dimensional cases, based on tensor product.Finally, numerical examples are provided to demonstrate the exponential convergence and efficiency of the algorithm.
The rest of this article is arranged as follows.In Section 2, we derive the equivalent hybrid format and its Legendre spectral approximation.In Section 3, we provide the error estimation of the approximating eigenvalues and eigenfunctions.In Section 4, we carry on a detailed description for the efficient implementation of the discrete variational form.In Section 5, we present some numerical examples to validate the theoretical findings and the effectiveness of the algorithm.Finally, we give in Section 6 a concluding remark.

Hybrid format and its Legendre spectral approximation
In this paper, we consider the fourth-order eigenvalue problem as follows: in Ω, (2.1) where both α and β are non-negative constants, and ) is a bounded domain.Let us introduce an auxiliary variable: The Eqs (2.1)-(2.3)can be restated as: in Ω, (2.6) Without loss of generality, we assume β > 0. If β = 0, we can add u(x) to both sides of the equation (2.5).At this time, only the corresponding eigenvalue becomes λ + 1, and the structure of the equation remains unchanged.By multiplying both sides of equation (2.6) with β, equations (2.5)-(2.7)can be rewritten as in Ω, (2.9) Let H m (Ω) and H m 0 (Ω) be the usual Sobolev spaces of order m, and their norms and seminorms are denoted by ∥ • ∥ m and | • | m , respectively.Especially, we denote L 2 (Ω) by H 0 (Ω), equipped by the inner product ⟨σ, ϱ⟩ := Ω σ ρdx and norm ∥σ∥ 0 = √ ⟨σ, σ⟩, here ρ denotes the complex conjugate of ϱ.Define the product Sobolev spaces as follows: and the corresponding norms are given by Then a variational formulation of (2.8)-(2.10)is: Find λ ∈ C and 0 (w, u) ∈ H 1 0 (Ω) such that where Define the finite element space , here P N denotes the space of Nth-order polynomials.Then the corresponding discrete variational form of (2.11) reads: Find λ N ∈ C and (2.12)
Lemma 2. B((w, u), (v, φ)) is a continuous bilinear functional defined on H 0 (Ω) × H 0 (Ω).That is, for any (w, u) ∈ H 0 (Ω) and (v, φ) ∈ H 0 (Ω), it holds Proof.Using Cauchy-Schwarz inequality and the definition of B, we have Note that the soure problem associated with (2.11) is to find (w, u) ∈ H 1 0 (Ω) such that Using Lemmas 1-2 and Lax-Milgram Theorem, the soure problem (3.4) exists an unique solution.Thus, we can define a solution operator T : H 0 (Ω) → H 1 0 (Ω) such that Thus, we can obtain the equivalent operator form of (2.11) as follows Note that the corresponding adjoint problem of (2.11) is: Find Similarly, the solution operator T * : H 0 (Ω) → H 1 0 (Ω) can be defined by Note that H 1 0 (Ω) is embedded in H 0 (Ω), together with (3.8), we can obtain the equivalent operator formulation of (3.7): From Lemmas 1-2 and Poincaré inequality, we have Assuming S is a bounded subset in H 1 0 (Ω), since H 1 0 (Ω) is compactly embeded in H 0 (Ω), then S is the sequentially compact set in H 0 (Ω).It follows from (3.9) that T S is a sequentially compact set in H 1 0 (Ω).Therefore, T : H 1 0 (Ω) → H 1 0 (Ω) is a complete continuous operator.We derive from (3.4) and (3.7) that which means that in the sense of inner product A(•, •), T * is the adjoint operator of T .Thus, T * : We can similarly define the discrete solution operator T N by N are all finite rank operators.Using (3.10), we can obtain the equivalent operator formulation of (2.12): Let us define a projection operator (3.12) Lemma 3. Let T and T N be the operators defined by (3.5) and (3.10), respectively.Then, it holds: By taking (v N , φ N ) = Π N T (w, u) − T N (w, u), we have Then, (3.13) follows from Lemma 1.
It is clear that With the theory of approximation, we have Let us define Then, Lemma 1 implies that ∥(w, u)∥ A is a equivalent norm in H 1 0 (Ω).Theorem 2. There holds: Proof.Based on the definition of operator norm, we have For ∀(v N , φ N ) ∈ W d N , we derive from Lemma 1 and (3.12) that That is, Thus, we have Then, (3.15) holds from (3.14).Note that the discrete variational form of (2.12) is: Find Define the discrete solution operator T * N : (3.17) Then, (3.17) implies that (3.16) has the following equivalent operator form: Let λ and µ be the nonzero eigenvalue with algebraic multiplicity g and the ascent of (λ −1 − T ), respectively.It follows from (3.15) that g eigenvalues λ j,N ( j = 1, 2, • • • , g) will converge to λ.Let ρ(T ) and σ(T ) be the resolvent set and the spectrum set, respectively.Define the spectral projection operators: where R z (T ) = (z − T ) −1 , and Γ lies in ρ(T ) and is a circle centered at λ −1 that does not enclose any other points within σ(T ).
According to [23], E is a projection onto the generalized eigenvectors space corresponding to λ −1 and T , that is, R(E) = N((λ −1 − T ) µ ), where R and N denote the range and the null space, respectively.
, where µ j is the ascent of λ −1 j,N − T N .For the dual problems (3.7) and (3.16), we can similarly define Defining the gaps between R(E) and R(E N ) in Based on the Theorems 8.1-8.4 in [23], we have the following prior error estimates.
Theorem 3.There exists a constant C such that Suppose for each N that (w N , u N ) satisfy ∥(w N , u N )∥ 1,Ω = 1 and (λ −1 N − T N ) k (w N , u N ) = 0 for some positive integer k ≤ µ.Then, for any integer l with k ≤ l ≤ µ, there exists a vector (w, u) such that (λ −1 − T ) l (w, u) = 0 and In order to offer the error estimates for the approximation of eigenvalues and eigenfunctions, we begin by introducing the d-dimensional Jacobian polynomial and weight function: where ). Define the non-uniformly weighted Sobolev space: with the following norm and semi-norm , where e j is the jth unit vector in From the theorem 8.1 and remark 8.14 in [18], we have following lemma.
Theorem 5.There exists a constant C such that Proof.We only give the proof for (3.19), and the same argument can be applied to (3.20).Using Poincaré inequality, we derive that Note that Through a similar derivation, we can obtain that Then, (3.19) follows from (3.5.32) in [18].The proof is completed.

Denote
We can obtain from Theorems 3-5 the error estimation of the approximating eigenvalues and eigenfunctions.Theorem 6.There exists a constant C such that Then, for any integer l with k ≤ l ≤ µ, there exists a vector (w, u) such that (λ −1 − T ) l (w, u) = 0 and

Efficient implementation of the discrete variational form.
In order to effectively solve problem (2.12), we first construct a set of basis functions for the approximation space W d N .Denote by L m (x) the Legendre polynomial of degree m.Let It is obvious that Denote • Case d = 2.We can expand the eigenfunctions as follows: where w i j , u i j are the expansion coefficients of the w N and u N , respectively. Denote .
We use ū to denote the vector formed by the columns of u.Now, plugging the expressions of (4.1) in (2.12), and taking (v N , φ N ) through all the basis functions in W 2 N , the discrete variational form (2.12) is equivalent to the following matrix form: where where A = (a i j ), B = (b i j ), ⊗ represents the tensor product symbol of the matrix.
• Case d = 3.Here, we can expand the eigenfunctions as follows: w i jk φ i (x 1 )φ j (x 2 )φ k (x 3 ), Denote by wk and ūk the vectors formed by the columns of w k and u k , respectively.Let W = ( w0 , w1 . Denote by W and Ū the vectors formed by the columns of W and U. Now, plugging the expressions of (4.3) in (2.12), and taking (v N , φ N ) through all the basis functions in W 3 N , we obtain the matrix form of the discrete variational form (2.12) as follows: where It is noted that each block matrix in (4.2) and (4.4) is sparse, and each non-zero element in them can be precisely calculated by utilizing the orthogonal properties of Legendre polynomials [18].Therefore, we can employ the sparse solver eigs(A, B, k, ′ sm ′ ) to effectively solve (4.2) and (4.4).

Numerical experiment
In this section, a series of numerical experiments will be presented to confirm the theoretical findings and demonstrate the efficiency of our algorithm.Our program is compiled and executed in MATLAB R2019a.
Example 1 We take Ω = (−1, 1) 2 , α = β = 1.The numerical results of the first fourth eigenvalues λ 1 N , λ 2 N (double eigenvalue), λ 3 N , λ 4 N for different N are listed in Table 1.To intuitively demonstrate the spectral accuracy of our algorithm, we employ the numerical solution with N = 40 as a reference solution and plot the absolute error curves of approximate eigenvalues as well as corresponding error curves under a log-log scale in Figure 1.Additionally, we also give an image of the reference solution for the eigenfunction and an error image between the reference solution and the approximate solution with N = 30 in Figure 2.   Table 1.Numerical results of the first four approximation eigenvalues for different N. We observe from Table 1 that the first four eigenvalues achieve at least 13-digit accuracy with N ≥ 25.Likewise, as shown in Figures 1-2, our algorithm is both convergent and spectral accurate.
As a comparison, we list in Table 2 the numerical results of the first four approximate eigenvalues obtained by directly solving the fourth-order eigenvalue problem using the classical Legendre spectral method.From Tables 1 and Table 2, we can observe that the convergence orders of the two numerical methods are almost the same.Howerver, for the fourth-order eigenvalue problem in general domain, if the spectral element method is directly applied to solve it, not only does the construction of the basis function become complex, but the computational load is also significant.On the contrary, the spectral element method for second-order problems is relatively mature in theoretical analysis and numerical calculation.Therefore, by transforming a fourth-order eigenvalue problem into a second-order coupled system, not only are the difficulties of theoretical analysis overcome, but the construction of the basis function is also relatively simple, facilitating efficient programming.
Although our theoretical analysis is based on the case where α and β are all positive constants, our algorithm is applicable to the case where α and β are variable coefficients.Thus, we provide a two-dimensional numerical example with variable coefficients in Example 2. Initially, we derive the equivalent matrix form for the case where α and β are all variable coefficients.Let us denote Similar to the deduction of (4.2), we can obtain the equivalent matrix form of the discrete variational form (2.12) as follows: where Example 2 We take Ω = (−1, 1) 2 , α = 1, β = e sin(x 1 +x 2 ) .The numerical results of the first fourth eigenvalues λ j N ( j = 1, 2, 3, 4) for different N are listed in Table 3.Similarly, in order to intuitively demonstrate the spectral accuracy of our algorithm, we employ the numerical solution with N = 40 as a reference solution and plot the absolute error curves of approximate eigenvalues as well as corresponding error curves under a log-log scale in Figure 3. Additionally, we also give an image of the reference solution for the eigenfunction and an error image between the reference solution and the approximate solution with N = 30 in Figure 4.
Table 3. Numerical results of the first four approximation eigenvalues for different N.   From Table 3, it is observed once again that the first four numerical eigenvalues arrive at an accuracy of approximately 14-digits when N ≥ 25.Additionally, it is observed from Figures 3-4 that our algorithm is convergent and spectral accurate.
Next, we provide a three-dimensional numerical example in Example 3.
Example 3 We take Ω = (−1, 1) 3 , α = β = 1.The numerical results of the first fourth eigenvalues λ j N ( j = 1, 2, 3, 4) for different N are listed in Table 4. Again, in order to intuitively demonstrate the spectral accuracy of our algorithm, we employ the numerical solution with N = 40 as a reference solution and plot absolute error curves of approximate eigenvalues as well as corresponding error curves under a log-log scale in Figure 5.As shown in Table 4, the first four eigenvalues achieve at least 14-digit accuracy when N ≥ 25.Furthermore, from Figure 5, we can see that our algorithm is also convergent and spectral accurate.

Conclusions
In this paper, an efficient Legendre spectral method is proposed and studied for fourth order eigenvalue problems with the boundary conditions of a simply supported plate.By introducing an auxiliary variable, the fourth order eigenvalue problem is transformed into a coupled second-order eigenvalue problem.By utilizing the hybrid format, a fresh weak formulation and its corresponding discrete variational form are formulated.Error estimates for the eigenvalues and eigenfunction approximations are also derived.In addition, numerical results validate the effectiveness of the algorithm and the correctness of theoretical results.
The algorithm proposed in this paper can be combined with the spectral element method to be applied to the numerical computation of fourth-order problems on more general domains.

Use of AI tools declaration
The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

N 4 Figure 1 .
Figure 1.Absolute error curves (left) between the numerical solution and the reference solution and the errors curves (right) under log-log scale.

Figure 2 .
Figure 2. Image (left) of the reference solution u 40 (x) and the error image (right) between reference solution and numerical solution u 30 (x).

4 Figure 3 .
Figure 3. Absolute error curves (left) between the numerical solution and the reference solution and the errors curves (right) under log-log scale.

Figure 4 .
Figure 4. Image (left) of the reference solution u 40 (x) the error image (right) between reference solution and numerical solution u 30 (x).

4 Figure 5 .
Figure 5. Absolute error curves (left) between the numerical solution and the reference solution and the errors curves (right) under log-log scale.

Table 2 .
Numerical results of the first four approximation eigenvalues for different N.

Table 4 .
Numerical results of the first four approximation eigenvalues for different N.