Loading [MathJax]/extensions/TeX/boldsymbol.js
Research article

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

  • Received: 13 January 2024 Revised: 08 February 2024 Accepted: 19 February 2024 Published: 22 February 2024
  • MSC : 65N15, 65N35

  • In this paper, we developed and studied 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 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 effective 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 efficiency of the algorithm.

    Citation: Yuanqiang Chen, Jihui Zheng, Jing An. A Legendre spectral method based on a hybrid format and its error estimation for fourth-order eigenvalue problems[J]. AIMS Mathematics, 2024, 9(3): 7570-7588. doi: 10.3934/math.2024367

    Related Papers:

    [1] Mohammad Reza Foroutan, Mir Sajjad Hashemi, Leila Gholizadeh, Ali Akgül, Fahd Jarad . A new application of the Legendre reproducing kernel method. AIMS Mathematics, 2022, 7(6): 10651-10670. doi: 10.3934/math.2022594
    [2] Siqin Tang, Hong Li . A Legendre-tau-Galerkin method in time for two-dimensional Sobolev equations. AIMS Mathematics, 2023, 8(7): 16073-16093. doi: 10.3934/math.2023820
    [3] Lingling Sun, Hai Bi, Yidu Yang . A posteriori error estimates of mixed discontinuous Galerkin method for a class of Stokes eigenvalue problems. AIMS Mathematics, 2023, 8(9): 21270-21297. doi: 10.3934/math.20231084
    [4] Xiaojun Zhou, Yue Dai . A spectral collocation method for the coupled system of nonlinear fractional differential equations. AIMS Mathematics, 2022, 7(4): 5670-5689. doi: 10.3934/math.2022314
    [5] Hui-qing Liao, Ying Fu, He-ping Ma . A space-time spectral method for the 1-D Maxwell equation. AIMS Mathematics, 2021, 6(7): 7649-7668. doi: 10.3934/math.2021444
    [6] Shixian Ren, Yu Zhang, Ziqiang Wang . An efficient spectral-Galerkin method for a new Steklov eigenvalue problem in inverse scattering. AIMS Mathematics, 2022, 7(5): 7528-7551. doi: 10.3934/math.2022423
    [7] Obaid Algahtani, M. A. Abdelkawy, António M. Lopes . A pseudo-spectral scheme for variable order fractional stochastic Volterra integro-differential equations. AIMS Mathematics, 2022, 7(8): 15453-15470. doi: 10.3934/math.2022846
    [8] Chenghua Gao, Maojun Ran . Spectral properties of a fourth-order eigenvalue problem with quadratic spectral parameters in a boundary condition. AIMS Mathematics, 2020, 5(2): 904-922. doi: 10.3934/math.2020062
    [9] Tingting Jiang, Jiantao Jiang, Jing An . An efficient Fourier spectral method and error analysis for the fourth order problem with periodic boundary conditions and variable coefficients. AIMS Mathematics, 2023, 8(4): 9585-9601. doi: 10.3934/math.2023484
    [10] Chuanhua Wu, Ziqiang Wang . The spectral collocation method for solving a fractional integro-differential equation. AIMS Mathematics, 2022, 7(6): 9577-9587. doi: 10.3934/math.2022532
  • In this paper, we developed and studied 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 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 effective 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 efficiency of the algorithm.



    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.

    In this paper, we consider the fourth-order eigenvalue problem as follows:

    Δ2u(x)αΔu(x)+βu(x)=λu(x),inΩ, (2.1)
    u(x)=0,onΩ, (2.2)
    Δu(x)=0,onΩ, (2.3)

    where both α and β are non-negative constants, and ΩRd(d=2,3) is a bounded domain. Let us introduce an auxiliary variable:

    Δu(x)+w(x)=0. (2.4)

    The Eqs (2.1)–(2.3) can be restated as:

    Δw(x)+αw(x)+βu(x)=λu(x),inΩ, (2.5)
    Δu(x)=w(x),inΩ, (2.6)
    w(x)=u(x)=0,onΩ. (2.7)

    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

    Δw(x)+αw(x)+βu(x)=λu(x),inΩ, (2.8)
    βΔu(x)=βw(x),inΩ, (2.9)
    w(x)=u(x)=0,onΩ. (2.10)

    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:

    H10(Ω):=H10(Ω)×H10(Ω),H0(Ω):=L2(Ω)×L2(Ω),

    and the corresponding norms are given by

    (w,u)1,Ω=(w21+u21)12,(w,u)0,Ω=(w20+u20)12.

    Then a variational formulation of (2.8)–(2.10) is: Find λC and 0(w,u)H10(Ω) such that

    A((w,u),(v,φ))=λB((w,u),(v,φ)),(v,φ)H10(Ω), (2.11)

    where

    A((w,u),(v,φ))=Ωwˉv+αwˉv+βuˉv+βuˉφβwˉφdx,B((w,u),(v,φ))=Ωuˉvdx.

    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 λNC and 0(wN,uN)WdN such that

    A((wN,uN),(vN,φN))=λNB((wN,uN),(vN,φN)),(vN,φN)WdN. (2.12)

    For simplicity, we limit our consideration to the case where Ω=(1,1)d with d=2,3, and ab means that acb, 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

    |A((w,u),(v,φ))|(w,u)1,Ω(v,φ)1,Ω,  (3.1)
    A((w,u),(w,u))(w,u)21,Ω. (3.2)

    Proof. Using Cauchy-Schwarz inequality, we derive that

    |A((w,u),(v,φ))|=|Ωwˉv+αwˉv+βuˉv+βuˉφβwˉφdx|Ω|wv|+α|wv|+β|uv|+β|uφ|β|wφ|dx(Ω|w|2+α|w|2+β|u|2+β|u|2+β|w|2dx)12×(Ω|v|2|+α|v|2+β|v|2+β|φ|2+β|φ|2dx)12(w,u)1,Ω(v,φ)1,Ω.

    Thus, (3.1) holds. On the other hand, from Poincaré inequality, we have

    A((w,u),(w,u))=Ω|w|2+α|w|2+βuw+β|u|2βwudx=Ω|w|2+α|w|2+β|u|2dx=|w|21+αw20+β|u|21(w,u)21,Ω.

    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

    |B((w,u),(v,φ))|(w,u)0,Ω(v,φ)0,Ω.  (3.3)

    Proof. Using Cauchy-Schwarz inequality and the definition of B, we have

    |B((w,u),(v,φ))|=|Ωuvdx|(Ω|u|2dx)12×(Ω|v|2dx)12(Ω|u|2+|w|2dx)12×(Ω|v|2+|φ|2dx)12=(w,u)0,Ω(v,φ)0,Ω.

    Note that the soure problem associated with (2.11) is to find (w,u)H10(Ω) such that

    A((w,u),(v,φ))=B((f,g),(v,φ)),(f,g)H0(Ω),(v,φ)H10(Ω). (3.4)

    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

    A(T(f,g),(v,φ))=B((f,g),(v,φ)),(f,g)H0(Ω),(v,φ)H10(Ω). (3.5)

    Thus, we can obtain the equivalent operator form of (2.11) as follows

    T(w,u)=λ1(w,u). (3.6)

    Note that the corresponding adjoint problem of (2.11) is: Find λC and 0(w,u)H10(Ω) such that

    A((v,φ),(w,u))=¯λB((v,φ),(w,u)),(v,φ)H10(Ω). (3.7)

    Similarly, the solution operator T:H0(Ω)H10(Ω) can be defined by

    A((v,φ),T(f,g))=B((v,φ),(f,g)),(f,g)H0(Ω),(v,φ)H10(Ω). (3.8)

    Note that H10(Ω) is embedded in H0(Ω), together with (3.8), we can obtain the equivalent operator formulation of (3.7):

    T(w,u)=(λ)1(w,u).

    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

    A(T(f,g),T(f,g))=B((f,g),T(f,g)).

    From Lemmas 1-2 and Poincaré inequality, we have

    T(f,g)21,ΩA(T(f,g),T(f,g))=B((f,g),T(f,g))(f,g)0,ΩT(f,g)0,Ω(f,g)0,ΩT(f,g)1,Ω,

    which implies that

    T(f,g)1,Ω(f,g)0,Ω. (3.9)

    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

    A(T(f,g),(v,φ))=B((f,g),(v,φ))=A((f,g),T(v,φ)),(f,g), (v,φ)H10(Ω),

    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

    A(TN(f,g),(vN,φN))=B((f,g),(vN,φN)),(f,g),H0(Ω), (vN,φN)WdN. (3.10)

    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):

    TN(wN,uN)=λ1N(wN,uN). (3.11)

    Let us define a projection operator ΠN:H10(Ω)WdN such that

    A((w,u)ΠN(w,u),(vN,φN))=0,(w,u)H10(Ω),(vN,φN)WdN. (3.12)

    Lemma 3. Let T and TN be the operators defined by (3.5) and (3.10), respectively. Then, it holds:

    TN=ΠNT. (3.13)

    Proof. For any (w,u)H10(Ω) and (vN,φN)WdN, we derive that

    A(ΠNT(w,u)TN(w,u),(vN,φN))=A(ΠNT(w,u)T(w,u)+T(w,u)TN(w,u),(vN,φN))=A(ΠNT(w,u)T(w,u),(vN,φN))+A(T(w,u)TN(w,u),(vN,φN))=0.

    By taking (vN,φN)=ΠNT(w,u)TN(w,u), we have

    A(ΠNT(w,u)TN(w,u),ΠNT(w,u)TN(w,u))=0.

    Then, (3.13) follows from Lemma 1.

    It is clear that TN|WdN:WdNWdN is a finite rank operator. Let

    ξN=sup(w,u)H10(Ω),(w,u)1,Ω=1inf(vN,φN)WdNT(w,u)(vN,φN)1,Ω.

    With the theory of approximation, we have

    limNξN=0. (3.14)

    Let us define

    (w,u)A=[A((w,u),(w,u))]12.

    Then, Lemma 1 implies that (w,u)A is a equivalent norm in H10(Ω).

    Theorem 2. There holds:

    limNTTNL(H10(Ω),H10(Ω))=0. (3.15)

    Proof. Based on the definition of operator norm, we have

    TTNL(H10(Ω),H10(Ω))=sup(w,u)H10(Ω),(w,u)1,Ω=1T(w,u)ΠNT(w,u)1,Ω.

    For (vN,φN)WdN, we derive from Lemma 1 and (3.12) that

    T(w,u)ΠNT(w,u)21,ΩT(w,u)ΠNT(w,u)2A=A(T(w,u)ΠNT(w,u),T(w,u)ΠNT(w,u))=A(T(w,u)ΠNT(w,u),T(w,u)(vN,φN))T(w,u)ΠNT(w,u)1,ΩT(w,u)(vN,φN)1,Ω.

    That is,

    T(w,u)ΠNT(w,u)1,ΩT(w,u)(vN,φN)1,Ω.

    Thus, we have

    TTNL(H10(Ω),H10(Ω))sup(w,u)H10(Ω),(w,u)1,Ω=1inf(vN,φN)WdNT(w,u)(vN,φN)1,Ω.

    Then, (3.15) holds from (3.14).

    Note that the discrete variational form of (2.12) is: Find λNC,0(vN,φN)WdN such that

    A((v,φ),(vN,φN))=¯λNB((v,φ),(vN,φN)),(vN,φN)WdN,(v,φ)WdN. (3.16)

    Define the discrete solution operator TN:H0(Ω)WdN by

    A((v,φ),TN(f,g))=B((v,φ),(f,g)),(f,g)H0(Ω),(v,φ)WdN. (3.17)

    Then, (3.17) implies that (3.16) has the following equivalent operator form:

    TN(wN,uN)=(λN)1(wN,uN).

    Let λ and μ be the nonzero eigenvalue with algebraic multiplicity g and the ascent of (λ1T), 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:

    E=12πiΓRz(T)dz,EN=12πiΓNRz(TN)dz,

    where Rz(T)=(zT)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((λ1T)μ), where R and N denote the range and the null space, respectively. Similarly, R(EN)=gj=1N((λ1j,NTN)μj), where μj is the ascent of λ1j,NTN. For the dual problems (3.7) and (3.16), we can similarly define E,R(E),EN and R(EN).

    For two closed subspaces Q1,Q2, let

    d(Q1,Q2)=sup(w,u)Q1,(w,u)1,Ω=1inf(v,φ)Q2(w,u)(v,φ)1,Ω.

    Defining the gaps between R(E) and R(EN) in H10(Ω) as follows

    δ(R(E),R(EN))=max{d(R(E),R(EN)),d(R(EN),R(E))}.

    Denote

    εN(λ)=sup(w,u)R(E),(w,u)1,Ω=1inf(vN,φN)WdN(w,u)(vN,φN)1,Ω,εN(λ)=sup(w,u)R(E),(w,u)1,Ω=1inf(vN,φN)WdN(w,u)(vN,φN)1,Ω.

    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

    δ(R(E),R(EN))CεN(λ),|λ(1ggj=1λ1j,N)1|CεN(λ)εN(λ),|λλj,N|C[εN(λ)εN(λ)]1μ.

    Theorem 4. If limNλN=λ. Suppose for each N that (wN,uN) satisfy (wN,uN)1,Ω=1 and (λ1NTN)k(wN,uN)=0 for some positive integer kμ. Then, for any integer l with klμ, there exists a vector (w,u) such that (λ1T)l(w,u)=0 and

    (w,u)(wN,uN)1,ΩC[εN(λ)]lk+1μ.

    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:

    \begin{align} &\mathit{\boldsymbol{J^{\alpha,\beta}_{n}}}(\mathit{\boldsymbol{x}}) = \prod\limits^{d}_{j = 1}\tilde{J}^{\alpha_{j},\beta_{j}}_{n_{j}}(x_{j}), \mathit{\boldsymbol{\omega^{ \alpha,\beta}}}(\mathit{\boldsymbol{x}}) = \prod\limits^{d}_{j = 1}\omega^{\alpha_{j},\beta_{j}}(x_{j}), \forall \mathit{\boldsymbol{x}}\in I^d, \end{align} (3.18)

    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:

    \begin{align*} & B^{s}_{\mathit{\boldsymbol{\alpha}},\mathit{\boldsymbol{\beta}}}(I^d): = \{ \rho:\partial^{\mathit{\boldsymbol{k}}}_{\mathit{\boldsymbol{x}}}\rho\in L^{2}_{\mathit{\boldsymbol{\omega}}^{\mathit{\boldsymbol{\alpha+k}},\mathit{\boldsymbol{\beta+k}}}}(I^d),\; 0\leq |\mathit{\boldsymbol{k}}|_{1} \leq s\},\; \; \forall s \in \mathbb{N}, \end{align*}

    with the following norm and semi-norm

    \begin{align*} &\|\rho\|_{ B^{s}_{\mathit{\boldsymbol{\alpha}},\mathit{\boldsymbol{\beta}}}(I^d)} = \bigg(\sum\limits_{0\leq |\mathit{\boldsymbol{k}}|_{1}\leq s}\|\partial^{\mathit{\boldsymbol{k}}}_{\mathit{\boldsymbol{x}}}\rho\|^2_{\mathit{\boldsymbol{\omega}}^{\mathit{\boldsymbol{\alpha+k}},\mathit{\boldsymbol{\beta+k}}}}\bigg)^{\frac{1}{2}}, |\rho|_{B^{s}_{\mathit{\boldsymbol{\alpha}},\mathit{\boldsymbol{\beta}}}(I^d)} = \bigg(\sum^{d}_{j = 1}\|\partial^{s}_{x_{j}}\rho\|^2_{\mathit{\boldsymbol{\omega}}^{\mathit{\boldsymbol{\alpha}}+s\mathit{\boldsymbol{e}}_{j},\mathit{\boldsymbol{\beta}}+s\mathit{\boldsymbol{e}}_{j}}}\bigg)^{\frac{1}{2}}\; , \end{align*}

    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

    \begin{align*} |\mathbf{\Pi}_N^{-\mathit{\boldsymbol{1}},-\mathit{\boldsymbol{1}}}\sigma-\sigma|_{ B^{1}_{-\mathit{\boldsymbol{1}},-\mathit{\boldsymbol{1}}}(I^d)}\leq c \sqrt{\frac{(N-s)!}{(N-1)!}}\bigg(N+s\bigg)^{\frac{1-s}{2}}|\sigma|_{B^{s}_{-\mathit{\boldsymbol{1}},-\mathit{\boldsymbol{1}}}(I^d)}, \end{align*}

    where c\simeq \sqrt{2} for N\gg1 .

    Theorem 5. There exists a constant C such that

    \begin{align} \varepsilon_N(\lambda)\leq&CN^{1-m}\sup\limits_{(w,u)\in R(E), \|(w,u)\|_{1,\Omega} = 1}[|w|_{B_{\mathit{\boldsymbol{-1}},\mathit{\boldsymbol{-1}}}^m(I^d)}+|u|_{B_{\mathit{\boldsymbol{-1}},\mathit{\boldsymbol{-1}}}^m(I^d)}], \end{align} (3.19)
    \begin{align} \varepsilon_N^*(\lambda^*)\leq&CN^{1-m}\sup\limits_{(w^*,u^*)\in R(E^*), \|(w^*,u^*)\|_{1,\Omega} = 1}[|w^*|_{B_{\mathit{\boldsymbol{-1}},\mathit{\boldsymbol{-1}}}^m(I^d)}+|u^*|_{B_{\mathit{\boldsymbol{-1}},\mathit{\boldsymbol{-1}}}^m(I^d)}]. \end{align} (3.20)

    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

    \begin{align*} \varepsilon_N(\lambda) = &\sup\limits_{(w,u)\in R(E),\|(w,u)\|_{1,\Omega} = 1}\inf\limits_{(v_N,\varphi_N)\in \mathbf{W}_N^d}\|(w,u)-(v_N,\varphi_N)\|_{1,\Omega}\\ \leq&C\sup\limits_{(w,u)\in R(E),\|(w,u)\|_{1,\Omega} = 1}\inf\limits_{(v_N,\varphi_N)\in \mathbf{W}_N^d}[|w-v_N|_1^2+|u-\varphi_N|_1^2]^{\frac{1}{2}}\\ \leq&C\sup\limits_{(w,u)\in R(E),\|(w,u)\|_{1,\Omega} = 1}[|w-\mathbf{\Pi}_N^{-\mathit{\boldsymbol{1}},-\mathit{\boldsymbol{1}}}w|_1^2+|u-\mathbf{\Pi}_N^{-\mathit{\boldsymbol{1}},-\mathit{\boldsymbol{1}}}u|_1^2]^{\frac{1}{2}}. \end{align*}

    Note that

    \begin{align*} &|w-\mathbf{\Pi}_N^{-\mathit{\boldsymbol{1}},-\mathit{\boldsymbol{1}}}w|_{1}^2 = \sum\limits_{j = 1}^{d}\int_{I^d}[\partial_{x_j}(w-\mathbf{\Pi}_N^{-\mathit{\boldsymbol{1}},-\mathit{\boldsymbol{1}}}w)]^2d\mathbf{x}\\ \leq&\sum\limits_{j = 1}^{d}\int_{I^d}[\partial_{x_j}(w-\mathbf{\Pi}_N^{-\mathit{\boldsymbol{1}},-\mathit{\boldsymbol{1}}}w)]^2\prod\limits_{i = 1,i\neq j}^{d}\frac{1}{(1-x_i)(1+x_i)}d\mathbf{x}\\ = &\sum\limits_{j = 1}^{d}\int_{I^d}[\partial_{x_j}(w-\mathbf{\Pi}_N^{-\mathbf{1},-\mathbf{1}}w)]^2\mathbf{\omega}^{-\mathit{\boldsymbol{1}}+\mathit{\boldsymbol{e}}_j,-\mathit{\boldsymbol{1}}+\mathit{\boldsymbol{e}}_j}d\mathbf{x} = |w-\mathbf{\Pi}_N^{-\mathit{\boldsymbol{1}},-\mathit{\boldsymbol{1}}}w|_{ B^{1}_{-\mathit{\boldsymbol{1}},-\mathit{\boldsymbol{1}}}(I^d)}^2. \end{align*}

    Through a similar derivation, we can obtain that

    \begin{align*} &|u-\mathbf{\Pi}_N^{-\mathbf{1},-\mathbf{1}}u|_{1}^2\leq|u-\mathbf{\Pi}_N^{-\mathit{\boldsymbol{1}},-\mathit{\boldsymbol{1}}}u|_{ B^{1}_{-\mathit{\boldsymbol{1}},-\mathit{\boldsymbol{1}}}(I^d)}^2. \end{align*}

    From Lemma 4, we derive that

    \begin{align*} \varepsilon_N(\lambda)&\leq C\sup\limits_{(w,u)\in R(E),\|(w,u)\|_{1,\Omega} = 1}[|w-\mathbf{\Pi}_N^{-\mathit{\boldsymbol{1}},-\mathit{\boldsymbol{1}}}w|_1^2+|u-\mathbf{\Pi}_N^{-\mathit{\boldsymbol{1}},-\mathbf{1}}u|_1^2]^{\frac{1}{2}}\\ &\leq C\sup\limits_{(w,u)\in R(E),\|(w,u)\|_{1,\Omega} = 1}[|w-\mathbf{\Pi}_N^{-\mathbf{1},-\mathbf{1}}w|_{B_{-\mathbf{1},-\mathbf{1}}^2}(I^d)+|u-\mathbf{\Pi}_N^{-\mathbf{1},-\mathbf{1}}u|_{B_{-\mathbf{1},-\mathbf{1}}^1}(I^d)]\\ &\leq C\sqrt{\frac{(N-m)!}{(N-1)!}}(N+m)^{\frac{1-m}{2}}\sup\limits_{(w,u)\in R(E),\|(w,u)\|_{1,\Omega} = 1}[|w|_{B_{-\mathbf{1},-\mathbf{1}}^m}(I^d)+|u|_{B_{-\mathbf{1},-\mathbf{1}}^m}(I^d)]. \end{align*}

    Then, (3.19) follows from (3.5.32) in [18]. The proof is completed.

    Denote

    \begin{align*} &\mathcal{M}_m(w,u) = \sup\limits_{(w,u)\in R(E), \|(w,u)\|_{1,\Omega} = 1}[|w|_{B_{\mathbf{-1},\mathbf{-1}}^m(I^d)}+|u|_{B_{\mathbf{-1},\mathbf{-1}}^m(I^d)}],\\ &\mathcal{M}_m^*(w*,u*) = \sup\limits_{(w^*,u^*)\in R(E^*), \|(w^*,u^*)\|_{1,\Omega} = 1}[|w^*|_{B_{\mathbf{-1},\mathbf{-1}}^m(I^d)}+|u^*|_{B_{\mathbf{-1},\mathbf{-1}}^m(I^d)}]. \end{align*}

    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

    \begin{align*} &\delta(R(E),R(E_N))\leq C N^{1-m}\mathcal{M}_m(w,u),\\ &|\tau-(\frac{1}{g}\sum\limits_{j = 1}^{g}\tau_{j,N}^{-1})^{-1}| \leq CN^{2(1-m)}\mathcal{M}_m(w,u)\mathcal{M}_m^*(w*,u*),\\ &|\lambda-\lambda_{j,N}|\leq C N^{\frac{2(1-m)}{\mu}}[\mathcal{M}_m(w,u)\mathcal{M}_m^*(w*,u*)]^{\frac{1}{\mu}}. \end{align*}

    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

    \begin{align*} \|(w,u)-(w_N,u_N)\|_{1,\Omega}\leq C[N^{1-m}\mathcal{M}_m(w,u)]^{\frac{l-k+1}{\mu}}. \end{align*}

    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

    \begin{align*} &\varphi_{m}(x) = L_{m}(x)-L_{m+2}(x), (m = 0,1,\cdots,N-2), \end{align*}

    It is obvious that

    \begin{align*} &\mathbf{W}_N^d = \big{\{}\prod\limits_{k = 1}^{d}\varphi_{ik}(x_k): ik = 0,1,\cdots,N-2\big{\}}\times \big{\{}\prod\limits_{k = 1}^{d}\varphi_{jk}(x_k): jk = 0,1,\cdots,N-2\big{\}}. \end{align*}

    Denote

    \begin{align*} &a_{ij} = \int_I\varphi_{j}^{'}(x)\varphi_{i}^{'}(x)dx, b_{ij} = \int_I\varphi_{j}(x)\varphi_{i}(x)dx. \end{align*}

    \bullet Case d = 2 . We can expand the eigenfunctions as follows:

    \begin{align} &(w_N,u_N) = \big{(}\sum\limits_{i,j = 0}^{N-2}w_{ij}\varphi_{i}(x_1)\varphi_{j}(x_2),\sum\limits_{i,j = 0}^{N-2}u_{ij}\varphi_{i}(x_1)\varphi_{j}(x_2)\big{)}, \end{align} (4.1)

    where w_{ij}, u_{ij} are the expansion coefficients of the w_N and u_N , respectively.

    Denote

    \mathbf{w} = \left ( \begin{array}{llll} w_{00}&\; w_{01}\; &\cdots & w_{0,N-2}\\ w_{10}&\; w_{11}\; &\cdots & w_{1,N-2}\\ \vdots&\; \vdots\; &\cdots & \vdots\\ w_{N-2,0}&\; w_{N-2,1}\; &\cdots & w_{N-2,N-2} \end{array} \right ) ,
    \mathbf{u} = \left ( \begin{array}{llll} u_{00}&\; u_{01}\; &\cdots & u_{0,N-2}\\ u_{10}&\; u_{11}\; &\cdots & u_{1,N-2}\\ \vdots&\; \vdots\; &\cdots & \vdots\\ u_{N-2,0}&\; u_{N-2,1}\; &\cdots & u_{N-2,N-2} \end{array} \right ) .

    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:

    \begin{equation} \left[ \begin{array}{cc} \mathcal{A} & \mathcal{B} \\ \mathcal{C} &\mathcal{D} \\ \end{array} \right] \left[ \begin{array}{c} \mathbf{\bar{w }}\\ \mathbf{ \bar{ u} } \end{array} \right] = \lambda_N \left[ \begin{array}{cc} \mathbf{0} & \mathcal{E} \\ \mathbf{0} &\mathbf{0 } \\ \end{array} \right] \left[ \begin{array}{c} \mathbf{\bar{w }}\\ \mathbf{ \bar{u} } \end{array} \right], \end{equation} (4.2)

    where

    \begin{align*} &\mathcal{A} = A\otimes B+B\otimes A-\alpha B\otimes B,\ \mathcal{B} = \beta B\otimes B,\\ &\mathcal{C} = -\beta B\otimes B,\ \mathcal{D} = \beta (A\otimes B+B\otimes A),\ \mathcal{E} = B\otimes B, \end{align*}

    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:

    \begin{align} &(w_N,u_N) = \big{(}\sum\limits_{i,j,k = 0}^{N-2}w_{ijk}\varphi_{i}(x_1)\varphi_{j}(x_2)\varphi_{k}(x_3),\sum\limits_{i,j,k = 0}^{N-2}u_{ijk}\varphi_{i}(x_1)\varphi_{j}(x_2)\varphi_{k}(x_3)\big{)}. \end{align} (4.3)

    Denote

    \mathbf{w}^k = \left ( \begin{array}{llll} w_{00}^k&\; w_{01}^k\; &\cdots & w_{0,N-2}^k\\ w_{10}^k&\; w_{11}^k\; &\cdots & w_{1,N-2}^k\\ \vdots&\; \vdots\; &\cdots & \vdots\\ w_{N-2,0}^k&\; w_{N-2,1}^k\; &\cdots & w_{N-2,N-2}^k \end{array} \right ) ,
    \mathbf{u}^k = \left ( \begin{array}{llll} u_{00}^k&\; u_{01}^k\; &\cdots & u_{0,N-2}^k\\ u_{10}^k&\; u_{11}^k\; &\cdots & u_{1,N-2}^k\\ \vdots&\; \vdots\; &\cdots & \vdots\\ u_{N-2,0}^k&\; u_{N-2,1}^k\; &\cdots & u_{N-2,N-2}^k \end{array} \right ) .

    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:

    \begin{equation} \left[ \begin{array}{cc} \mathbb{A} & \mathbb{B} \\ \mathbb{C} &\mathbb{D} \\ \end{array} \right] \left[ \begin{array}{c} \mathbf{W}\\ \mathbf{ U } \end{array} \right] = \lambda_N \left[ \begin{array}{cc} \mathbf{0} & \mathbb{E} \\ \mathbf{0} &\mathbf{0 } \\ \end{array} \right] \left[ \begin{array}{c} \mathbf{W}\\ \mathbf{U} \end{array} \right], \end{equation} (4.4)

    where

    \begin{align*} \mathbb{A}& = B\otimes B\otimes A+B\otimes A\otimes B+A\otimes B\otimes B-\alpha B\otimes B\otimes B,\\ \mathbb{B}& = \beta B\otimes B\otimes B ,\ \mathbb{C} = -\beta B\otimes B\otimes B,\\ \mathbb{D}& = \beta (B\otimes B\otimes A+B\otimes A\otimes B+A\otimes B\otimes B),\ \mathbb{E} = B\otimes B\otimes B. \end{align*}

    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).

    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.

    Table 1.  Numerical results of the first four approximation eigenvalues for different N .
    N \lambda_{N}^1 \lambda_{N}^2 \lambda_{N}^3 \lambda_{N}^4
    5 30.287864207529200 163.1204205202208 402.7547327452907 813.0000000000005
    10 30.287074959045135 165.5387105473934 410.3755739012024 634.4780697093954
    15 30.287074959045280 165.5387102419903 410.3755729381884 634.4808299653813
    20 30.287074959045280 165.5387102419902 410.3755729381863 634.4808299652357
    25 30.287074959045283 165.5387102419905 410.3755729381896 634.4808299652358

     | Show Table
    DownLoad: CSV
    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.  Image (left) of the reference solution u_{40}(\mathbf{x}) and the error image (right) between reference solution and numerical solution u_{30}(\mathbf{x}) .

    We observe from Table 1 that the first four eigenvalues achieve at least 13-digit accuracy with N\geq25 . Likewise, as shown in Figures 12, 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.

    Table 2.  Numerical results of the first four approximation eigenvalues for different N .
    N \lambda_{N}^1 \lambda_{N}^2 \lambda_{N}^3 \lambda_{N}^4
    5 30.295236624303840 165.2737789038451 407.8269176502486 1247.000000000008
    10 30.287074959045206 165.5387151406006 410.3755834469236 634.4806909655549
    15 30.287074959045313 165.5387102419901 410.3755729381873 634.4808299678238
    20 30.287074959045288 165.5387102419903 410.3755729381873 634.4808299652353
    25 30.287074959045280 165.5387102419902 410.3755729381866 634.4808299652374

     | Show Table
    DownLoad: CSV

    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

    \begin{align*} &q_{(N-1)j+(n+1),(N-1)k+(m+1)} = \int_{-1}^{1}\int_{-1}^{1}\varphi_{k}^{'}\varphi_{j}^{'}\varphi_{m}\varphi_{n}+\varphi_{k}\varphi_{j}\varphi_{m}^{'}\varphi_{n}^{'}dx_1dx_2,\\ &c_{(N-1)j+(n+1),(N-1)k+(m+1)} = \int_{-1}^{1}\int_{-1}^{1}\alpha\varphi_{k}\varphi_{j}\varphi_{m}\varphi_{n}dx_1dx_2. \end{align*}

    Let \beta_1 = \partial_{x_1}\beta, \beta_2 = \partial_{x_2}\beta . We further denote

    \begin{align*} &d_{(N-1)j+(n+1),(N-1)k+(m+1)} = \int_{-1}^{1}\int_{-1}^{1}\beta_1\varphi_{k}^{'}\varphi_{j}\varphi_{m}\varphi_{n}dx_1dx_2,\\ &e_{(N-1)j+(n+1),(N-1)k+(m+1)} = \int_{-1}^{1}\int_{-1}^{1}\beta_2\varphi_{k}\varphi_{j}\varphi_{m}^{'}\varphi_{n}dx_1dx_2,\\ &f_{(N-1)j+(n+1),(N-1)k+(m+1)} = \int_{-1}^{1}\int_{-1}^{1}\beta(\varphi_{k}^{'}\varphi_{j}^{'}\varphi_{m}\varphi_{n}+\varphi_{k}\varphi_{j}\varphi_{m}^{'}\varphi_{n}^{'})dx_1dx_2,\\ &g_{(N-1)j+(n+1),(N-1)k+(m+1)} = \int_{-1}^{1}\int_{-1}^{1}\beta\varphi_{k}\varphi_{j}\varphi_{m}\varphi_{n}dx_1dx_2,\\ &h_{(N-1)j+(n+1),(N-1)k+(m+1)} = \int_{-1}^{1}\int_{-1}^{1}\varphi_{k}\varphi_{j}\varphi_{m}\varphi_{n}dx_1dx_2. \end{align*}

    Similar to the deduction of (4.2), we can obtain the equivalent matrix form of the discrete variational form (2.12) as follows:

    \begin{equation} \left[ \begin{array}{cc} \mathcal{A}_v & \mathcal{B}_v \\ \mathcal{C}_v &\mathcal{D}_v \\ \end{array} \right] \left[ \begin{array}{c} \mathbf{\bar{w }}\\ \mathbf{ \bar{ u} } \end{array} \right] = \lambda_N \left[ \begin{array}{cc} \mathbf{0} & \mathcal{E}_v \\ \mathbf{0} &\mathbf{0 } \\ \end{array} \right] \left[ \begin{array}{c} \mathbf{\bar{w }}\\ \mathbf{ \bar{u} } \end{array} \right], \end{equation} (5.1)

    where

    \begin{align*} \mathcal{A}_v = Q+C,\ \mathcal{B}_v = G,\ \mathcal{C}_v = D+E+F,\ \mathcal{D}_v = -G ,\ \mathcal{E}_v = H,\\ \end{align*}

    where

    \begin{align*} &C = (c_{(N-1)j+(n+1),(N-1)k+(m+1)}), D = (d_{(N-1)j+(n+1),(N-1)k+(m+1)}), \\ &E = (e_{(N-1)j+(n+1),(N-1)k+(m+1)}),F = (f_{(N-1)j+(n+1),(N-1)k+(m+1)}), \\ &G = (g_{(N-1)j+(n+1),(N-1)k+(m+1)}), H = (h_{(N-1)j+(n+1),(N-1)k+(m+1)}),\\ &Q = (q_{(N-1)j+(n+1),(N-1)k+(m+1)}). \end{align*}

    \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.

    Table 3.  Numerical results of the first four approximation eigenvalues for different N .
    N \lambda_{N}^1 \lambda_{N}^2 \lambda_{N}^3 \lambda_{N}^4
    5 20.516080766173257 140.7497111800699 140.9491722973433 370.2555125438320
    10 20.523346905996156 140.9328397952852 141.1180042820287 371.0947733339725
    15 20.523346901558348 140.9328457799942 141.1180110253167 371.0947913721103
    20 20.523346901558394 140.9328457799939 141.1180110253170 371.0947913721107
    25 20.523346901558362 140.9328457799946 141.1180110253167 371.0947913721118

     | Show Table
    DownLoad: CSV
    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.  Image (left) of the reference solution u_{40}(\mathbf{x}) and the error image (right) between reference solution and numerical solution u_{30}(\mathbf{x}) .

    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.

    Table 4.  Numerical results of the first four approximation eigenvalues for different N .
    N \lambda_{N}^1 \lambda_{N}^2 \lambda_{N}^3 \lambda_{N}^4
    5 48.391913926971800 202.6713333320417 463.7525086494870 831.6354398793130
    10 48.390410405809130 205.3660485651845 471.9269144900231 710.5119176894398
    15 48.390410405809410 205.3660482248708 471.9269134571854 710.5148388418000
    20 48.390410405809240 205.3660482248708 471.9269134571845 710.5148388416472
    25 48.390410405809270 205.3660482248702 471.9269134571852 710.5148388416419

     | Show Table
    DownLoad: CSV
    Figure 5.  Absolute error curves (left) between the numerical solution and the reference solution and the errors curves (right) under log-log scale.

    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.

    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.

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

    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).

    The authors declare to have no competing interests.



    [1] G. Engel, K. Garikipati, T. J. R. Hughes, M. G. Larson, L. Mazzei, R. L. Taylor, Continuous/discontinuous finite element approximations of fourth-order elliptic problems in structural and continuum mechanics with applications to thin beams and plates, and strain gradient elasticity, Comput. Method. Appl. Mech. Eng., 191 (2002), 3669–3750. https://doi.org/10.1016/S0045-7825(02)00286-4 doi: 10.1016/S0045-7825(02)00286-4
    [2] B. Li, G. Fairweather, B. Bialecki, Discrete-time orthogonal spline collocation methods for vibration problems, SIAM J. Numer. Anal., 39 (2002), 2045–2065. https://doi.org/10.1137/S0036142900348729 doi: 10.1137/S0036142900348729
    [3] J. Shen, J. Xu, J. Yang, The scalar auxiliary variable (SAV) approach for gradient flows, J. Comput. Phys., 353 (2018), 407–416. https://doi.org/10.1016/j.jcp.2017.10.021 doi: 10.1016/j.jcp.2017.10.021
    [4] J. Shen, J. Xu, J. Yang, A new class of efficient and robust energy stable schemes for gradient flows, SIAM Rev., 61 (2019), 474–506. https://doi.org/10.1137/17M1150153 doi: 10.1137/17M1150153
    [5] J. An, J. Shen, Spectral approximation to a transmission eigenvalue problem and its applications to an inverse problem, Comput. Math. Appl., 69 (2015), 1132–1143. https://doi.org/10.1016/j.camwa.2015.03.002 doi: 10.1016/j.camwa.2015.03.002
    [6] T. Tan, W. X. Cao, J. An, Spectral approximation based on a mixed scheme and its error estimates for transmission eigenvalue problems, Comput. Math. Appl., 111 (2022), 20–33. https://doi.org/10.1016/j.camwa.2022.02.009 doi: 10.1016/j.camwa.2022.02.009
    [7] S. X. Ren, T. Tan, J. An, An efficient spectral-Galerkin approximation based on dimension reduction scheme for transmission eigenvalues in polar geometries, Comput. Math. Appl., 80 (2020), 940–955. https://doi.org/10.1016/j.camwa.2020.05.018 doi: 10.1016/j.camwa.2020.05.018
    [8] N. Peng, C. Wang, J. An, An efficient finite-element method and error analysis for the fourth-order elliptic equation in a circular domain, Int. J. Comput. Math., 99 (2022), 1785–1802. https://doi.org/10.1080/00207160.2021.2007240 doi: 10.1080/00207160.2021.2007240
    [9] L. Ge, H. F. Niu, J. W. Zhou, Convergence analysis and error estimate for distributed optimal control problems governed by Stokes equations with velocity-constraint, Adv. Appl. Math. Mech., 14 (2022), 33–55. https://doi.org/10.4208/aamm.OA-2020-0302 doi: 10.4208/aamm.OA-2020-0302
    [10] H. F. Niu, D. P. Yang, J. W. Zhou, Numerical analysis of an optimal control problem governed by the stationary Navier-Stokes equations with global velocity-constrained, Commun. Comput. Phys., 24 (2018), 1477–1502. https://doi.org/10.4208/cicp.oa-2017-0045 doi: 10.4208/cicp.oa-2017-0045
    [11] J. G. Sun, Iterative methods for transmission eigenvalues, SIAM J. Numer. Anal., 49 (2011), 1860–1874. https://doi.org/10.1137/100785478 doi: 10.1137/100785478
    [12] L. Li, J. An, An efficient spectral method and rigorous error analysis based on dimension reduction scheme for fourth order problems, Numer. Methods PDE, 37 (2021), 152–171. https://doi.org/10.1002/num.22523 doi: 10.1002/num.22523
    [13] Y. P. Chen, J. W. Zhou, Error estimates of spectral Legendre-Galerkin methods for the fourth-order equation in one dimension, Appl. Math. Comput., 268 (2015), 1217–1226. https://doi.org/10.1016/j.amc.2015.06.082 doi: 10.1016/j.amc.2015.06.082
    [14] B. Bialecki, A. Karageorghis, A Legendre spectral Galerkin method for the biharmonic Dirichlet problem, SIAM J. Sci. Comput., 22 (2000), 1549–1569. https://doi.org/10.1137/S1064827598342407 doi: 10.1137/S1064827598342407
    [15] P. E. Bj\phirstad, B. P. Tj\phistheim, Efficient algorithms for solving a fourth-order equation with spectral-Galerkin method, SIAM J. Sci. Comput., 18 (1997), 621–632. https://doi.org/10.1137/S1064827596298233 doi: 10.1137/S1064827596298233
    [16] E. H. Doha, A. H. Bhrawy, Efficient spectral-Galerkin algorithms for direct solution of fourth-order differential equations using Jacobi polynomials, Appl. Numer. Math., 58 (2008), 1224–1244. https://doi.org/10.1016/j.apnum.2007.07.001 doi: 10.1016/j.apnum.2007.07.001
    [17] J. T. Jiang, J. An, J. W. Zhou, A novel numerical method based on a high order polynomial approximation of the fourth order Steklov equation and its eigenvalue problems, Discrete Cont. Dyn. B, 28 (2023), 50–69. https://doi.org/10.3934/dcdsb.2022066 doi: 10.3934/dcdsb.2022066
    [18] J. Shen, T. Tang, L. L. Wang, Spectral methods: algorithms, analysis and applications, Springer Science and Business Media, 2011. https://citations.springernature.com/book?doi = 10.1007/978-3-540-71041-7
    [19] J. Shen, T. Tang, Spectral and High-Order Methods with Applications, Science Press, 2006. https://api.semanticscholar.org/CorpusID: 117589335
    [20] E. Blåsten, X. F. Li, H. Y. Liu, Y. L. Wang, On vanishing and localizing of transmission eigenfunctions near singular points: A numerical study, Inverse Probl., 33 (2017), 105001. https://iopscience.iop.org/article/10.1088/1361-6420/aa8826/pdf
    [21] Y. Gao, H. Y. Liu, X. C. Wang, K. Zhang, On an artificial neural network for inverse scattering problems, J. Comput. Phys., 448 (2022), 110771. https://doi.org/10.1016/j.jcp.2021.110771 doi: 10.1016/j.jcp.2021.110771
    [22] Y. W. Yin, W. S. Yin, P. C. Meng, H. Y. Liu, The interior inverse scattering problem for a two-layered cavity using the Bayesian method, Inverse Probl. Imaging, 16 (2022), 673–690. https://doi.org/10.3934/ipi.2021069 doi: 10.3934/ipi.2021069
    [23] I. Babuka, J. E. Osborn, Eigenvalue problems, in: P.G. Ciarlet, J.L. Lions(Ed.), FEM (Part 1), Handb. Numer. Anal., 2 (1991), 640–787. http://refhub.elsevier.com/S0898-1221(22)00066-9/bibE5231A3DB5C6FA293BE28291A40673D9s1
  • Reader Comments
  • © 2024 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1301) PDF downloads(39) Cited by(0)

Figures and Tables

Figures(5)  /  Tables(4)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog