四阶特征值问题基于降阶格式的一种有效的Legendre-Galerkin逼近
An Efficient Legendre-Galerkin Approximation Based on Reduced Order Scheme for Fourth Order Eigenvalue Problems
摘要: 本文提出了四阶特征值问题基于降阶格式的一种有效的Legendre-Galerkin逼近。首先,我们引入了一个辅助函数,将原问题转化为一个二阶混合格式。通过引入一些适当的Sobolev空间,其相应的变分形式被建立,并在解足够光滑条件下证明了其等价性。其次,基于Legendre多项式的正交性质,两组紧凑的基函数被构造,并导出具有稀疏系数矩阵的线性特征系统。最后,我们给出了两个数值例子,数值结果表明了算法的收敛性与高精度。
Abstract: In this paper, an efficient Legendre-Galerkin approximation based on reduced order scheme for fourth order eigenvalue problems is presented. First, we introduce an auxiliary function to trans-form the original problem into a second order mixed format. By introducing some suitable Sobolev Spaces, the corresponding variational form is established, and its equivalence is proved if the solu-tion is sufficiently smooth. Secondly, based on the orthogonal property of Legendre polynomials, two groups of compact basis functions are constructed, and a linear characteristic system with sparse coefficients matrix is derived. Finally, we give two numerical examples, and the numerical results show the convergence and high precision of the algorithm.
文章引用:魏涛. 四阶特征值问题基于降阶格式的一种有效的Legendre-Galerkin逼近[J]. 应用数学进展, 2023, 12(4): 1981-1988. https://doi.org/10.12677/AAM.2023.124203

1. 引言

特征值问题广泛地应用于科学和工程领域,如量子力学,流体力学和随机过程 [1] [2] [3] [4] [5] 等。其数值计算方法一直备受很多学者关注。因此,提出一些有效求解特征值问题的高精度数值方法是非常有意义的。

到目前为止,已经有一些有限元法 [6] [7] [8] ,有限差分法 [9] [10] [11] 相继被提出并应用于求解特征值问题。然而,这些数值方法要获得高精度的数值解需要花费大量计算时间和内存容量。众所周知,目前很少有将谱方法应用于求解基于混合格式的四阶特征值问题。因此,本文提出了四阶特征值问题基于降阶格式的一种有效的Legendre-Galerkin逼近。首先,我们引入了一个辅助函数,将原问题转化为一个二阶混合格式。通过引入一些适当的Sobolev空间,其相应的变分形式被建立,并在解足够光滑条件下证明了其等价性。其次,基于Legendre多项式的正交性质,两组紧凑的基函数被构造,并导出具有稀疏系数矩阵的线性特征系统。最后,我们给出了两个数值例子,数值结果表明了算法的收敛性与高精度。

本文其余部分组织如下:在第2节,我们推导了四阶特征值问题等价的二阶耦合格式及其变分形式。在第3节,我们给出了算法的有效实现。第4节,我们给出了一些数值例子。第5节,我们给出了结论性注记。

2. 耦合的降阶格式及其变分形式

本文考虑如下的四阶特征值问题:

{ Δ 2 ψ α Δ ψ + β ψ = λ ψ , Ω , ψ = ψ n = 0 , Ω , (1)

其中 α , β 为非负常数, Ω d 中的有界区域, n 表示单位外法向量。

ω = Δ ψ ,则方程组(1)等价于如下的二阶混合格式:

{ Δ ψ = ω , Ω , Δ ω + α ω + β ψ = λ ψ , Ω , ψ = ψ n = 0 , Ω . (2)

L 2 ( Ω ) H m ( Ω ) H 0 m ( Ω ) 为通常的Sobolev空间, m 分别为 L 2 ( Ω ) H m ( Ω ) H 0 m ( Ω ) 中的范数, | | m H m ( Ω ) H 0 m ( Ω ) 中的半范数。记

V = H 1 ( Ω ) , M = H 0 1 ( Ω ) , H = L 2 ( Ω ) ,

a ( ω , v ) = Ω ω v d x , b ( ω , ψ ) = Ω ω ψ d x .

则 的弱形式为:求 ( λ , ω , ψ ) × V × M 使得:

{ a ( ω , v ) + b ( v , ψ ) = 0 , v V , b ( ω , ξ ) + ( α ω , ξ ) + ( β ψ , ξ ) = ( λ ψ , ξ ) , ξ M . (3)

不失一般性,我们仅考虑 Ω = [ 1 , 1 ] 2 的情况。设 P N 为区间 [ 1 , 1 ] 上的N次多项式空间,定义逼近空间 V N = ( P N × P N ) H 1 ( Ω ) M N = ( P N × P N ) H 0 1 ( Ω ) ,则(3)式的离散格式为:求 ( λ N , ω N , ψ N ) × V N × M N 使得:

{ a ( ω N , v N ) + b ( v N , ψ N ) = 0 , v N V N , b ( ω N , ξ N ) + ( α ω N , ξ N ) + ( β ψ N , ξ N ) = ( λ N ψ N , ξ N ) , ξ N M N . (4)

定理1:四阶特征值问题(1)与混合变分问题(3)等价。

证明:若 ψ 为(1)的解,由格林公式可知 ( ω = Δ ψ , ψ ) 为(3)的解。反之,若 ( ω , ψ ) 为(3)的解,则对 ξ H 0 1 ( Ω ) ,由格林公式可得:

b ( ω , ξ ) = Ω ω ξ d x = Ω ω n ξ d S Ω Δ ω ξ d x = Ω Δ ω ξ d x .

将上式代入(3)的第二个等式可得:

Ω ( Δ ω + α ω + β ψ λ ψ ) ξ d x = 0.

由变分法基本引理可得:

Δ ω + α ω + β ψ = λ ψ . (5)

类似地,由格林公式,变分法基本引理和(3)的第一个等式可得:

ω = Δ ψ . (6)

另一方面,对于 v H 1 ( Ω ) ,由(3)的第一个等式可得:

Ω ( ω + Δ ψ ) v d x Ω ψ n v d S = 0. (7)

在(7)中取 v = ψ n ,结合(6)可得: ψ | Ω = ψ n | Ω = 0 。再将(6)代入(5)可知 ψ 为(1)的解。证毕。

3. 算法的有效实现

在这节,我们将详细描述算法的实现过程。首先构造逼近空间中的一组基函数,我们用 L n ( x 1 ) 表示n次Legendre多项式。令

η i ( x 1 ) = L i ( x 1 ) L i + 2 ( x 1 ) , ( i = 0 , 1 , , N 2 ) ;

ξ i ( x 1 ) = η i ( x 1 ) , ( i = 0 , 1 , , N 2 ) , ξ N 1 ( x 1 ) = 1 , ξ N ( x 1 ) = 1 + x 1 2 .

则有

M N = span { η i ( x 1 ) η j ( x 2 ) , i , j = 0 , 1 , , N 2 } ;

V N = span { ξ i ( x 1 ) ξ j ( x 2 ) , i , j = 0 , 1 , , N } .

s k j = 1 1 ξ j ξ k d x 1 , a k j = 1 1 ξ j ξ k d x 1 , b k j = 1 1 η j ξ k d x 1 ,

c k j = 1 1 η j ξ k d x 1 , d k j = 1 1 ξ j η k d x 1 , e k j = 1 1 ξ j η k d x 1 ,

h k j = 1 1 η j η k d x 1

ω N = i , j = 0 N ω i j ξ i ( x 1 ) ξ j ( x 2 ) ,

ψ N = i , j = 0 N 2 ψ i j η i ( x 1 ) η j ( x 2 ) ,

W = ( ω 00 ω 01 ω 0 N ω 10 ω 11 ω 1 N ω N 0 ω N 1 ω N N ) ,

Ψ = ( ψ 00 ψ 01 ψ 0 , N 2 ψ 10 ψ 11 ψ 1 , N 2 ψ N 2 , 0 ψ N 2 , 1 ψ N 2 , N 2 ) .

W ¯ 表示由W的列构成的长度为 ( N + 1 ) 2 的列向量。用 Ψ ¯ 表示由 Ψ 的列构成的长度为 ( N 1 ) 2 的列向量,则对

v N = ξ m ( x 1 ) ξ n ( x 2 ) , m , n = 0 , 1 , 2 , , N , ξ N = η k ( x 1 ) η l ( x 2 ) , k , l = 0 , 1 , 2 , , N 2 ,

我们有

a ( ω N , v N ) = Ω ω N v N d x = i , j = 0 N 1 1 1 1 ξ i ( x 1 ) ξ j ( x 2 ) ξ m ( x 1 ) ξ n ( x 2 ) d x ω i j = i , j = 0 N 1 1 1 1 ξ i ( x 1 ) ξ m ( x 1 ) ξ j ( x 2 ) ξ n ( x 2 ) d x ω i j = i , j = 0 N 1 1 ξ i ( x 1 ) ξ m ( x 1 ) d x 1 1 1 ξ j ( x 2 ) ξ n ( x 2 ) d x 2 ω i j = i , j = 0 N a m i a n j ω i j = A ( m , : ) W A ( n , : ) T = A ( n , : ) A ( m , : ) W ¯ .

b ( v N , ψ N ) = Ω v N ψ N d x = i , j = 0 N 2 1 1 1 1 ( ξ m ( x 1 ) ξ n ( x 2 ) , ξ m ( x 1 ) ξ n ( x 2 ) ) ( ( η i ( x 1 ) η j ( x 2 ) , η i ( x 1 ) η j ( x 2 ) ) d x ψ i j = i , j = 0 N 2 1 1 1 1 ( ξ m ( x 1 ) ξ n ( x 2 ) η i ( x 1 ) η j ( x 2 ) + ξ m ( x 1 ) ξ n ( x 2 ) η i ( x 1 ) η j ( x 2 ) ) d x ψ i j = i , j = 0 N 2 1 1 1 1 ( ξ m ( x 1 ) η i ( x 1 ) ξ n ( x 2 ) η j ( x 2 ) + ξ m ( x 1 ) η i ( x 1 ) ξ n ( x 2 ) η j ( x 2 ) ) d x ψ i j = i , j = 0 N 2 ( b m i c n j + c m i b n j ) ψ i j = ( B ( m , : ) Ψ C ( n , : ) T + C ( m , : ) Ψ B ( n , : ) T ) = ( C ( n , : ) B ( m , : ) + B ( n , : ) C ( m , : ) ) Ψ ¯ .

b ( ω N , ξ N ) = Ω ω N ξ N d x = i , j = 0 N 1 1 1 1 ( ξ i ( x 1 ) ξ j ( x 2 ) , ξ i ( x 1 ) ξ j ( x 2 ) ) ( ( η k ( x 1 ) η l ( x 2 ) , η k ( x 1 ) η l ( x 2 ) ) d x ω i j = i , j = 0 N 1 1 1 1 ( ξ i ( x 1 ) ξ j ( x 2 ) η k ( x 1 ) η l ( x 2 ) + ξ i ( x 1 ) ξ j ( x 2 ) η k ( x 1 ) η l ( x 2 ) ) d x ω i j = i , j = 0 N 1 1 1 1 ( ξ i ( x 1 ) η k ( x 1 ) ξ j ( x 2 ) η l ( x 2 ) + ξ i ( x 1 ) η k ( x 1 ) ξ j ( x 2 ) η l ( x 2 ) ) d x ω i j = i , j = 0 N ( d k i e l j + e k i d l j ) ω i j = ( D ( k , : ) W E ( l , : ) T + E ( k , : ) W D ( l , : ) T ) = ( E ( l , : ) D ( k , : ) + D ( l , : ) E ( k , : ) ) W ¯ .

( ω N , ξ N ) = Ω ω N ξ N d x = i , j = 0 N 1 1 1 1 ξ i ( x 1 ) ξ j ( x 2 ) η k ( x 1 ) η l ( x 2 ) d x ω i j = i , j = 0 N 1 1 1 1 ξ i ( x 1 ) η k ( x 1 ) ξ j ( x 2 ) η l ( x 2 ) d x ω i j = i , j = 0 N 1 1 ξ i ( x 1 ) η k ( x 1 ) d x 1 1 1 ξ j ( x 2 ) η l ( x 2 ) d x 2 ω i j = i , j = 0 N e k i e l j ω i j = E ( k , : ) W E ( l , : ) T = E ( l , : ) E ( k , : ) W ¯ .

( ψ N , ξ N ) = Ω ψ N ξ N d x = i , j = 0 N 2 1 1 1 1 η i ( x 1 ) η j ( x 2 ) η k ( x 1 ) η l ( x 2 ) d x ψ i j = i , j = 0 N 2 1 1 1 1 η i ( x 1 ) η k ( x 1 ) η j ( x 2 ) η l ( x 2 ) d x ψ i j = i , j = 0 N 2 1 1 η i ( x 1 ) η k ( x 1 ) d x 1 1 1 η j ( x 2 ) η l ( x 2 ) d x 2 ψ i j = i , j = 0 N 2 h k i h l j ψ i j = H ( k , : ) Ψ H ( l , : ) T = H ( l , : ) H ( k , : ) Ψ ¯ .

其中 S = ( s i j ) , A = ( a i j ) , B = ( b i j ) , C = ( c i j ) , D = ( d i j ) , E = ( e i j ) , H = ( h i j ) A ( m , : ) 表示矩阵 A = ( a i j ) 的第m行, 表示矩阵的张量积,即 A B = ( a i j B ) 。记

T 1 = A ( n , : ) A ( m , : ) , T 2 = ( C ( n , : ) B ( m , : ) + B ( n , : ) C ( m , : ) ) , T 3 = E ( l , : ) D ( k , : ) + D ( l , : ) E ( k , : ) , T 4 = α E ( l , : ) E ( k , : ) , T 5 = β H ( l , : ) H ( k , : ) , T 6 = H ( l , : ) H ( k , : ) .

因此(4)等价于下面的矩阵形式:

( T 1 T 2 T 3 + T 4 T 5 ) ( W ¯ Ψ ¯ ) = ( O O O λ N T 6 ) ( W ¯ Ψ ¯ )

4. 数值例子

在这一节,我们将在MATLAB R2017a平台上进行两个数值例子,通过数值结果来验证我们提出的算法是一种高精度数值方法。

例1:我们的第一个算例取 α = 1 β = 1 ,用第3节提出的算法求解问题(1),对于不同的N,我们在表1中列出了前4个特征值的数值结果( λ N 2 , λ N 3 为重根)。另外,我们以 N = 60 时的数值解作为一个参考解,并在图1中列出了 λ N i ( i = 1 , 2 , 3 , 4 ) 关于不同的N的误差曲线。

Table 1. For different N, the first four approximate eigenvalue numerical results

表1. 对于不同的N,前4个逼近特征值的数值结果

Figure 1. The error between the approximation eigenvalue λ N i ( i = 1 , 2 , 3 , 4 ) and the reference solution

图1. 逼近特征值 λ N i ( i = 1 , 2 , 3 , 4 ) 与参考解之间的误差

表1图1中可以观察到我们的算法是收敛的和高精度的。

例2:我们取 α = β = 0 ,用第3节提出的算法求解问题(1),对于不同的N,我们在表2中列出了前4个特征值的数值结果( λ N 2 , λ N 3 为重根)。我们还以 N = 60 时的数值解作为一个参考解,并在图2中列出了 λ N i ( i = 1 , 2 , 3 , 4 ) 关于不同的N的误差曲线。

Table 2. For different N, the first four approximate eigenvalue numerical results

表2. 对于不同的N,前4个逼近特征值的数值结果

Figure 2. The error between the approximation eigenvalue λ N i ( i = 1 , 2 , 3 , 4 ) and the reference solution

图2. 逼近特征值 λ N i ( i = 1 , 2 , 3 , 4 ) 与参考解之间的误差

表2图2中再次观察到我们的算法是收敛的和高精度的。

5. 结论性注记

本文提出了四阶特征值问题基于降阶格式的一种有效的谱Galerkin逼近。首先,我们建立了原问题的弱形式及其离散格式。其次,我们给出了算法的实现过程。最后通过具体的数值算例来验证算法的有效性和理论结果的正确性。

本文只给出了二维区域 [ 1 , 1 ] 2 上基于降阶格式的四阶特征值问题的数值试验结果,但文中的方法也可以直接推广到一般的矩形区域以及长方体区域。

本文通过将一个四阶特征值问题转化为一个二阶混合特征值问题,不仅降低了问题的复杂度,还可以结合谱元法应用于一般区域上的四阶特征值问题的计算。

参考文献

[1] 刘诗焕, 朱先阳. 高维空间中阻尼Boussinesq方程初值问题的整体解[J]. 西南师范大学学报(自然科学版), 2018, 43(9): 1-5.
[2] 张琪慧, 尚月强. Navier-Stokes方程的亚格子模型后处理混合有限元方法[J]. 西南大学学报(自然科学版), 2019, 41(3): 67-74.
[3] Hu, J., Huang, Y.Q. and Shen, H.M. (2004) The Lower Approximation of Eigen-value by Lumped Mass Finite Element Method. Journal of Computational Mathematics, 22, 545-556.
[4] Grebenkov, D.S. and Nguyen, B.T. (2013) Geometrical Structure of Laplacian Eigenfunctions. SIAM Review, 55, 601-667.
https://doi.org/10.1137/120880173
[5] Boffi, D. (2010) Finite Element Approximation of Eigenvalue Problems. Acta Numerica, 19, 1-120.
https://doi.org/10.1017/S0962492910000012
[6] Zhou, J.-W., Zhang, J. and Xing, X.-Q. (2016) Galerkin Spec-tral Approximations for Optimal Contral Problems Governed by the Four Order Equation with an Integral Constraint on State. Computers & Mathematics with Applications, 72, 2549-2561.
https://doi.org/10.1016/j.camwa.2016.08.009
[7] Davis, C.B. (2014) A Partition of Unity Method with Penalty for Fourth Order Problems. Journal of Scientific Computing, 60, 228-248.
https://doi.org/10.1007/s10915-013-9795-8
[8] Sun, J.-G. (2012) A New Family of High Regularity Elements. Numerical Methods for Partial Differential Equations, 28, 1-16.
https://doi.org/10.1002/num.20601
[9] Mercier, B., Osborn, J., Rappaz, J., et al. (1981) Eigenvalue Approximation by Mixed and Hybrid Methods. Mathematics of Computation, 36, 427-453.
https://doi.org/10.1090/S0025-5718-1981-0606505-9
[10] Bjorstad, P.E. and Tjos-theim, B.P. (1997) Timely Communication: Efficient Algorithms for Solving a Fourth Order Equation with the Spec-tral-Galerkin Method. SIAM Journal on Scientific Computing, 18, 621-632.
https://doi.org/10.1137/S1064827596298233
[11] Lv, T., Lin, Z.-B. and Shi, J.-M. (1988)A Fourth Order Finite Difference Approximation to the Eigenvalues of a Clamped Plate. Journal of Computational Mathematics, 6, 267-271.