1.
Introduction
Fourth-order equations and eigenvalue problems have been extensively utilized in many scientific and engineering fields [1,2], and the numerical computation of numerous intricate and nonlinear fourth-order equations and eigenvalue problems, such as the Cahn-Hilliard equation and the transmission eigenvalue problem [3,4,5,6,7], stem from repeatedly solving a linear fourth-order equation and eigenvalue problem. Numerous theoretical and numerical computing results have been obtained on fourth-order equations and eigenvalue problems, primarily involving the finite element methods [8,9,10,11] and spectral methods [12,13,14,15,16,17].
For the finite element method applied to fourth-order equations and eigenvalue problems, a C1 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.
2.
Hybrid format and its Legendre spectral approximation
In this paper, we consider the fourth-order eigenvalue problem as follows:
where both α and β are non-negative constants, and Ω⊂Rd(d=2,3) is a bounded domain. Let us introduce an auxiliary variable:
The Eqs (2.1)–(2.3) can be restated as:
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
Let Hm(Ω) and Hm0(Ω) be the usual Sobolev spaces of order m, and their norms and seminorms are denoted by ‖⋅‖m and |⋅|m, respectively. Especially, we denote L2(Ω) by H0(Ω), 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)∈H10(Ω) such that
where
Define the finite element space WdN=H10(Ω)∩(PdN×PdN), here PN denotes the space of Nth-order polynomials. Then the corresponding discrete variational form of (2.11) reads: Find λN∈C and 0≠(wN,uN)∈WdN such that
3.
Error estimation of the approximate solutions
For simplicity, we limit our consideration to the case where Ω=(−1,1)d with d=2,3, and a≲b means that a≤cb, where c is a positive constant independent of N.
Lemma 1. A((w,u),(v,φ)) is a continuous and coercive bilinear functionals defined on H10(Ω)×H10(Ω). That is, for any given ((w,u),(v,φ))∈H10(Ω)×H10(Ω), there hold
Proof. Using Cauchy-Schwarz inequality, we derive that
Thus, (3.1) holds. On the other hand, from Poincaré inequality, we have
This proof is completed.
Lemma 2. B((w,u),(v,φ)) is a continuous bilinear functional defined on H0(Ω)×H0(Ω). That is, for any (w,u)∈H0(Ω) and (v,φ)∈H0(Ω), 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)∈H10(Ω) 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:H0(Ω)→H10(Ω) 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 λ∗∈C and 0≠(w∗,u∗)∈H10(Ω) such that
Similarly, the solution operator T∗:H0(Ω)→H10(Ω) can be defined by
Note that H10(Ω) is embedded in H0(Ω), together with (3.8), we can obtain the equivalent operator formulation of (3.7):
Theorem 1. Both T:H10(Ω)→H10(Ω) and T∗:H10(Ω)→H10(Ω) are complete continuous operators.
Proof. We can obtain by taking (v,φ)=T(f,g) in (3.5) that
From Lemmas 1-2 and Poincaré inequality, we have
which implies that
Assuming S is a bounded subset in H10(Ω), since H10(Ω) is compactly embeded in H0(Ω), then S is the sequentially compact set in H0(Ω). It follows from (3.9) that TS is a sequentially compact set in H10(Ω). Therefore, T:H10(Ω)→H10(Ω) 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∗:H10(Ω)→H10(Ω) is also a complete continuous operator.
We can similarly define the discrete solution operator TN by
It is obvious that both TN:H10(Ω)→WdN and H0(Ω)→WdN are all finite rank operators. Using (3.10), we can obtain the equivalent operator formulation of (2.12):
Let us define a projection operator ΠN:H10(Ω)→WdN such that
Lemma 3. Let T and TN be the operators defined by (3.5) and (3.10), respectively. Then, it holds:
Proof. For any (w,u)∈H10(Ω) and (vN,φN)∈WdN, we derive that
By taking (vN,φN)=ΠNT(w,u)−TN(w,u), we have
Then, (3.13) follows from Lemma 1.
It is clear that TN|WdN:WdN→WdN is a finite rank operator. Let
With the theory of approximation, we have
Let us define
Then, Lemma 1 implies that ‖(w,u)‖A is a equivalent norm in H10(Ω).
Theorem 2. There holds:
Proof. Based on the definition of operator norm, we have
For ∀(vN,φN)∈WdN, 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 λ∗N∈C,0≠(v∗N,φ∗N)∈WdN such that
Define the discrete solution operator T∗N:H0(Ω)→WdN by
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 Rz(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. Similarly, R(EN)=g∑j=1N((λ−1j,N−TN)μj), where μj is the ascent of λ−1j,N−TN. For the dual problems (3.7) and (3.16), we can similarly define E∗,R(E∗),E∗N and R(E∗N).
For two closed subspaces Q1,Q2, let
Defining the gaps between R(E) and R(EN) in H10(Ω) as follows
Denote
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
Theorem 4. If limN→∞λN=λ. Suppose for each N that (wN,uN) satisfy ‖(wN,uN)‖1,Ω=1 and (λ−1N−TN)k(wN,uN)=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 \mathit{\boldsymbol{n}} = (n_1, n_2, \cdots, n_d)\in \mathbb{N}^d , \mathit{\boldsymbol{\alpha}} = (\alpha_1, \alpha_2, \cdots, \alpha_d) , \mathit{\boldsymbol{\beta}} = (\beta_1, \beta_2, \cdots, \beta_d), I = (-1, 1) . Define the non-uniformly weighted Sobolev space:
with the following norm and semi-norm
where \mathit{\boldsymbol{e}}_{j} is the j th unit vector in \mathbb{R}^{d} , \mathit{\boldsymbol{k}} = (k_1, k_2, \cdots, k_d), |\mathit{\boldsymbol{k}}|_{1} = k_1+k_2+\cdots+k_d . From the theorem 8.1 and remark 8.14 in [18], we have following lemma.
Lemma 4. There exists a projection operator \mathbf{\Pi}_N^{-\mathit{\boldsymbol{1}}, -\mathit{\boldsymbol{1}}}: \mathbf{H}^1_0(\Omega)\rightarrow \mathbf{W}^d_N , such that for any \sigma\in B^{s}_{-\mathit{\boldsymbol{1}}, -\mathit{\boldsymbol{1}}}(\mathcal{D}) , it holds for 1\leq s\leq N+1 that
where c\simeq \sqrt{2} for N\gg1 .
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
From Lemma 4, we derive 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
Theorem 7. If \lim\limits_{N\rightarrow \infty}\lambda_{N} = \lambda . Suppose for each N that (w_N, u_N) satisfy \|(w_N, u_N)\|_{1, \Omega} = 1 and (\lambda_{N}^{-1}-T_N)^k(w_N, u_N) = 0 for some positive integer k\leq \mu . Then, for any integer l with k\leq l \leq \mu , there exists a vector (w, u) such that (\lambda^{-1}-T)^l(w, u) = 0 and
4.
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 \mathbf{W}_N^d . Denote by L_m(x) the Legendre polynomial of degree m . Let
It is obvious that
Denote
\bullet Case d = 2 . We can expand the eigenfunctions as follows:
where w_{ij}, u_{ij} are the expansion coefficients of the w_N and u_N , respectively.
Denote
We use \mathbf{\bar{u}} to denote the vector formed by the columns of \mathbf{u} . Now, plugging the expressions of (4.1) in (2.12), and taking (v_N, \varphi_N) through all the basis functions in \mathbf{W}_N^2 , the discrete variational form (2.12) is equivalent to the following matrix form:
where
where A = (a_{ij}), B = (b_{ij}) , \otimes represents the tensor product symbol of the matrix.
\bullet Case d = 3 . Here, we can expand the eigenfunctions as follows:
Denote
Denote by \bar{\mathbf{w}}^k and \bar{\mathbf{u}}^k the vectors formed by the columns of \mathbf{w}^k and \mathbf{u}^k , respectively. Let \mathbf{W} = (\bar{\mathbf{w}}^0, \bar{\mathbf{w}}^1, \cdots, \bar{\mathbf{w}}^{N-2}) , \mathbf{U} = (\bar{\mathbf{u}}^0, \bar{\mathbf{u}}^1, \cdots, \bar{\mathbf{u}}^{N-2}) . Denote by \mathbf{\bar{W}} and \mathbf{\bar{U}} the vectors formed by the columns of \mathbf{W} and \mathbf{U} . Now, plugging the expressions of (4.3) in (2.12), and taking (v_N, \varphi_N) through all the basis functions in \mathbf{W}_N^3 , 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).
5.
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.
\mathbf{Example\; \; 1} We take \Omega = (-1, 1)^2 , \alpha = \beta = 1 . The numerical results of the first fourth eigenvalues \lambda^{1}_{N}, \lambda^{2}_{N} (double eigenvalue), \lambda^{3}_{N} , \lambda^{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.
We observe from Table 1 that the first four eigenvalues achieve at least 13-digit accuracy with N\geq25 . 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 \alpha and \beta are all positive constants, our algorithm is applicable to the case where \alpha and \beta 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 \alpha and \beta are all variable coefficients. Let us denote
Let \beta_1 = \partial_{x_1}\beta, \beta_2 = \partial_{x_2}\beta . We further 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
where
\mathbf{Example\; \; 2} We take \Omega = (-1, 1)^2 , \alpha = 1, \beta = e^{\sin(x_1+x_2)} . The numerical results of the first fourth eigenvalues \lambda^{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.
From Table 3, it is observed once again that the first four numerical eigenvalues arrive at an accuracy of approximately 14 -digits when N \geq 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.
\mathbf{Example\; \; 3} We take \Omega = (-1, 1)^3 , \alpha = \beta = 1 . The numerical results of the first fourth eigenvalues \lambda^{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\geq25 . Furthermore, from Figure 5, we can see that our algorithm is also convergent and spectral accurate.
6.
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.
Acknowledgments
The authors would like to thank the anonymous referees for their valuable comments and suggestions, which led to considerable improvement of the article. The research is supported by National Natural Science Foundation of China (Grant No. 11661022).
Conflict of interest
The authors declare to have no competing interests.