四阶问题有着重要的应用背景, 很多复杂的非线性问题的计算最终都归结于反复地求解一个四阶问题, 如Cahn–Hilliard方程[1-4], 声波传输问题[5-9]等.关于四阶问题的理论分析和数值计算已有很多成果, 主要包括有限元法[10-13], 谱方法及一些高阶数值方法[14-19]. 众所周知, 谱方法具有谱精度等优点, 已广泛地应用于微分方程的数值计算[20-23], 但它要求计算区域必须是方形或立方体, 为了克服这个缺陷, 通常采取谱元法去求解一般区域上的微分方程. 对于二阶问题的谱元方法, 其理论和数值计算都已经比较成熟.关于四阶问题的谱方法, 文献[20]提出了圆域上四阶Steklov方程及其特征值问题的一种有效的谱方法, 文献[21]提出了一些特殊区域上双调和特征值问题的一种谱逼近方法, 文献[22]提出了方形和立方体区域上传输特征值问题的一种有效的谱Galerkin逼近. 对于方形、立方体和一些特殊区域上四阶问题的谱方法还有很多研究成果, 就不在此一一列举了. 然而, 关于四阶问题的谱元法的文献很少, 虽然文献[24]提出了Stokes特征值问题的三角谱元方法, 文献[25]提出四阶方程的一种$ C^1 $协调四边形谱元方法, 但是对于一般区域上四阶问题的谱元法, 其基函数的构造非常复杂, 计算量也比较大. 因此, 提出一些四阶问题基于降阶格式的谱方法是非常有意义的. 据我们所知, 很少有关于简支板边界条件下四阶问题基于降阶格式的谱方法的报道.
本文的目标是提出简支板边界条件下四阶问题基于降阶格式的一种有效的谱Galerkin逼近. 首先, 通过引入一个辅助函数, 将原问题化为两个耦合的二阶问题. 其次, 通过引入适当的Sobolev空间, 建立两个耦合二阶问题的弱形式和相应的离散格式, 并利用Lax-Milgram定理, 证明弱解和逼近解的存在唯一性. 然后, 我们引入多维的非一致带权Sobolev空间, 根据其投影算子的逼近性质, 证明逼近解的误差估计. 最后, 利用Legendre多项式的正交性质构造一组适当的基函数, 推导离散格式基于张量积的矩阵形式, 并通过数值算例验证算法的有效性和理论结果的正确性. 作为一个模型, 本文考虑如下的四阶问题:
其中$ \boldsymbol{x}\in \mathbb{R}^d(d=2, 3) $, $ D $为$ \mathbb{R}^d $中的一个有界区域, $ \partial D $表示区域$ D $的边界.
本文剩余部分组织如下: 在第2节我们将推导原问题等价的二阶耦合格式及其变分形式. 在第3节, 我们将证明逼近解的误差估计. 在第4节, 我们将详细描述算法的有效实现过程. 在第5节我们将给出一些数值算例. 在第6节我们给出了一些结论性评注.
首先通过引入一个辅助变量, 将方程(1.1)化为两个耦合的二阶问题, 再引入适当的Sobolev空间及其逼近空间, 建立两个耦合二阶问题的弱形式和相应的离散格式. 令$ w(\boldsymbol{x})=-\Delta u(\boldsymbol{x}). $则方程(1.1)可化为下面两个耦合的二阶问题:
和
引入通常的Sobolev空间:
其相应的内积和范数分别为:
不失一般性, 本文仅考虑$ \psi(\boldsymbol{x})=\varphi(\boldsymbol{x})=0 $, $ D=(-1, 1)^d (d=2, 3) $的情形. 由格林公式可知方程(2.1)-(2.2)的弱形式为: 找$ (w, u) \in H_0^1(D)\times H_0^1(D) $, 使得
其中
令$ P_N $为$ N $次多项式空间, 定义逼近空间: $ X_N=(P_N)^d\cap H_{0}^1(D) $. 则(2.3)-(2.4)的离散格式为: 找$ (w_N, u_N) \in X_N\times X_N $, 使得
在这一节, 我们将证明弱解和逼近解的存在唯一性以及它们之间的误差估计.
引理1 $ a(w,v) $为定义在$ H_{0}^{1}(D)\times H_{0}^{1}(D) $上的正定连续的双线性泛函, 即
证 由Cauchy-Schwarz不等式有
另一方面, 由于$ |w|^2=| \int_{-1}^{x_i}\partial_{ x_i}wdx_i|^2\leq(x_i+1)\int_{-1}^{1}(\partial_{ x_i}w)^2dx_i, $将该不等式两边在区域$ D $上积分得
由此可得$ d\|w\|^{2}\leq 2|w|_{1}^{2}. $从而有
证毕.
引理2 若$ f\in L^2(D) $, 则$ F(v) $为$ H_0^1(D) $上的有界线性泛函. 即$ |F(v)|\lesssim\|v\|_1. $
证 对$ \forall v\in H_0^1(D) $, 由Schwarz不等式有
证毕. 由引理1, 引理2和Lax-Milgram定理可得下面的定理.
定理1 若$ f\in L^2(D) $, 则弱形式(2.3)-(2.4)存在唯一解$ (w,u) \in H_0^1(D)\times H_0^1(D) $, 离散格式(2.5)- (2.6)存在唯一解$ (w_N,u_N) \in \textbf{X}_N\times \textbf{X}_N $.
引理3 设$ w $和$ w_N $分别为弱形式(2.3)和离散格式(2.5)的解, 则有
证 由(2.3)和(2.5), 我们有
由(3.1)减去(3.2)得
由(3.3)和Schwarz不等式有
即
由$ v_N $的任意性可得
证毕. 我们用$ J_{n}^{\alpha,\beta} $表示定义在区间$ I:=(-1,1) $上关于权函数$ \omega^{\alpha,\beta}(x)=(1-t)^\alpha (1+t)^\beta $正交的Jacobi多项式, 即
其中$ t_{n}^{\alpha,\beta}=\|J_{n}^{\alpha,\beta}\|_{\omega^{\alpha,\beta}}^{2},(\alpha,\beta>-1 $). 对于整数指数对$ (l,k) $, 引入如下的广义Jacobi多项式$ ^{[26,27]} $:
其中$ n\geq n_0 $, 且
令$ \hat{J}_{n}^{l,k}(t) $为标准的广义Jacobi多项式, 即
定义$ d $维张量形式的广义Jacobi多项式和权函数:
其中$ \boldsymbol{n}=(n_1,n_2,...,n_d),\; \boldsymbol{l}=(l_1,l_2,...,l_d),\; \boldsymbol{k}=(k_1,k_2,...,k_d) $. 定义$ d $维$ N $次多项式空间:
由$ d $维广义Jacobi多项式和逼近空间$ \boldsymbol{X}_N $的定义可知: $ \boldsymbol{X}_N=\boldsymbol{Q}_{N}^{\mathit{\pmb{-1,-1}}} $. 定义正交投影算子: $ T_{N}^{\mathit{\pmb{-1,-1}}}:\; \boldsymbol{L}_{\boldsymbol\omega^{\mathit{\pmb{-1,-1}}}}^{2}(\boldsymbol{I}^d)\rightarrow \boldsymbol{Q}_{N}^{\mathit{\pmb{-1,-1}}} $, 使得
定义$ d $维带权$ Sobolev $空间, 即
其相应的半范数和范数分别为:
其中$ \; \boldsymbol{e}_\boldsymbol{i}\; $是$ \mathbb{R}^d\; $中的第$ \; \boldsymbol{i}\; $个单位向量, $ |\boldsymbol{m}|_1=\sum\limits_{i=1}^dm_i,\; \partial_{\boldsymbol{x}}^{\boldsymbol{m}}w=\partial_{x_1}^{m_1}\partial_{x_2}^{m_2}...\partial_{x_d}^{m_d}w $. 由文献[28]中定理8.1和评注8.14我们有下面的引理.
引理4 对$ \forall w \in \boldsymbol{B}_{\mathit{\pmb{-1,-1}}}^s\left(\boldsymbol{I}^d\right) $和$ 1\leq s\leq N+1 $下面的不等式成立
其中, 当$ N\gg1 $时, $ C\simeq\sqrt{2} $.
定理2 设$ w $和$ w_N $分别为弱形式(2.3)和离散格式(2.5)的解, 则当$ w\in {\boldsymbol{B}_{\mathit{\pmb{-1,-1}}}^{s}}(\boldsymbol{I}^d)\cap H_0^1(\boldsymbol{I}^d) $和$ 1\leq s\leq N+1 $时, 下面的估计式成立
证 由引理3可知
则有
由于
则由引理4和文献[28]中的(3.5.32)可得
定理3 设$ u $和$ u_N $分别为弱形式(2.4)和离散格式(2.6)的解, 则当$ u\in {\boldsymbol{B}_{\mathit{\pmb{-1,-1}}}^{s}}(\boldsymbol{I}^d)\cap H_0^1(\boldsymbol{I}^d) $和$ 1\leq s\leq N+1 $时, 下面的估计式成立
证 由(2.4)和(2.6)可知:
由(3.6)减去(3.7)可得:
对$ \forall q_N\in \boldsymbol{X}_N $, 由引理1和(3.8)有
在上式中取$ q_N=\bf{T}_{N}^{\mathit{\pmb{-1,-1}}}u $, 则有
又由于
则由Poincaré不等式, 引理4和定理2有
在这节我们将详细地描述算法的有效实现过程, 首先构造逼近空间的一组基函数, 令
其中$ L_i(t) $表示$ i $次Legendre多项式, 则逼近空间为$ X_N=\big\{u_N(\boldsymbol{x}):u_N(\boldsymbol{x})\in[P_N^0]^d\big\}. $令
则由文献[29]中的引理$ 2.1 $可知:
下面我们将推导离散格式(2.5)和(2.6)基于张量积的矩阵形式.当d=2时, 对于离散格式(2.5), 我们将寻找$ w_N=\sum_{i, j=0}^{N-2} w_{ij}\varphi_i(x_1)\varphi_j(x_2). $令
用$ \overline{W} $表示由$ W $的列构成的长度为$ (N-1)^2 $的列向量, 取$ v_N=\varphi_l(x_1)\varphi_k(x_2), (l, k=0, 1, \cdots, N-2 $), 则有:
其中$ S(l, :), M(k, :) $分别表示矩阵$ S=\{s_{ij}\}, M=\{m_{ij}\} $的第$ l $和$ k $行, $ \otimes $表示矩阵的张量积符号, 即$ S\otimes M=(s_{ij}M)_{i, j=0}^{N-2} $; 对于离散格式(2.6), 我们将寻找$ u_N=\sum_{i, j=0}^{N-2} u_{ij}\varphi_i(x_1)\varphi_j(x_2). $令
用$ \overline{U} $表示由$ U $的列构成的长度为$ (N-1)^2 $的列向量, 取$ h_N=\varphi_l(x_1)\varphi_k(x_2), (l, k=0, 1, ..., N-2) $, 类似地, 我们可以得到
则离散格式(2.5)基于张量积的矩阵形式为$ (M\otimes S+S\otimes M)\overline{W}=F, $其中
离散格式(2.6)基于张量积的矩阵形式为$ (M\otimes S+S\otimes M)\overline{U}=\hat{F}, $其中
当$ d=3 $时, 对于离散格式(2.5), 我们将寻找$ w_N=\sum_{i, j, q=0}^{N-2} w_{ij}^q\varphi_i(x_1)\varphi_j(x_2)\varphi_q(x_3). $令
用$ \overline{W}^q $表示$ W^q $的列构成的长度为$ (N-1)^2 $的列向量, 令$ W=(\overline{W}^0, \overline{W}^1, \cdots, \overline{W}^N) $, $ \overline{W} $表示由矩阵$ W $的列构成的长度为$ (N-1)^3 $的列向量. 取$ v_N=\varphi_l(x_1)\varphi_k(x_2)\varphi_n(x_3), (l, k, n=0, 1, \cdots, N-2) $, 则有:
对于离散格式(2.6), 我们将寻找$ u_N=\sum_{i, j, q=0}^{N-2} u_{ij}^q\varphi_i(x_1)\varphi_j(x_2)\varphi_q(x_3). $令
用$ \overline{U}^q $表示$ W^q $的列构成的长度为$ (N-1)^2 $的列向量, 令$ U=(\overline{U}^0, \overline{U}^1, \cdots, \overline{U}^N) $, $ \overline{U} $表示由矩阵$ U $的列构成的长度为$ (N-1)^3 $的列向量.取$ h_N=\varphi_l(x_1)\varphi_k(x_2)\varphi_n(x_3), (l, k, n=0, 1, \cdots, N-2) $, 类似地, 可以得到:
则离散格式(2.5)基于张量积的矩阵形式为$ (M\otimes M\otimes S+S\otimes M\otimes M+M\otimes S\otimes M)\overline{W}=F, $其中
离散格式(2.6)基于张量积的矩阵形式为$ (M\otimes M\otimes S+S\otimes M\otimes M+M\otimes S\otimes M)\overline{U}=\hat{F}, $其中
为了表明算法的有效性, 我们将进行一系列的数值实验. 我们在MATLAB R2016b平台上进行编程计算.
例1: 取精确解$ u(\boldsymbol{x})=\sin{\pi x_1}\sin{\pi x_2} $, 则$ w(\boldsymbol{x})=-\Delta u(\boldsymbol{x})=2\pi^2 \sin{\pi x_1}\sin{\pi x_2} $, $ u $显然满足边界条件, 将$ w(\boldsymbol{x}), u(\boldsymbol{x}) $代入方程(2.1)和(2.2)可得到$ f(\boldsymbol{x}) $. 我们分别在表格1和表格2中列出了精确解$ w(\boldsymbol{x}), u(\boldsymbol{x}) $与逼近解$ w_N(\boldsymbol{x}), u_N(\boldsymbol{x}) $在三种模下的误差. 从表格1和表格2可以观察到, 当$ N\geq20 $时, 精确解$ w(\boldsymbol{x}), u(\boldsymbol{x}) $与逼近解$ w_N(\boldsymbol{x}), u_N(\boldsymbol{x}) $在三种模下的误差都达到了大约$ 10^{-13} $的精度. 为了进一步表明我们算法的高精度, 我们在图 1、图 2和图 3、图 4中画出了当$ N $取不同值时精确解与逼近解之间的误差图像, 从图 1、图 2和图 3、图 4可以看出我们的算法是收敛的和高精度的. 另外, 作为比较, 我们在表格3中列出了基于四阶格式的谱方法在三种模下的误差结果. 从表 2和表 3可以观察到, 两种数值格式都具有谱精度.
例2: 取精确解$ u(\boldsymbol{x})=\sin{\pi x_1}\sin{\pi x_2}\sin{\pi x_3} $, 则$ w(\boldsymbol{x})=3\pi^2 \sin{\pi x_1}\sin{\pi x_2}\sin{\pi x_3} $, $ u(\boldsymbol{x}) $显然满足边界条件, 将$ w(\boldsymbol{x}), u(\boldsymbol{x}) $代入方程(2.1)和(2.2)可得到$ f(\boldsymbol{x}) $. 我们分别在表格4和表格5中列出了精确解$ w(\boldsymbol{x}), u(\boldsymbol{x}) $与逼近解$ w_N(\boldsymbol{x}), u_N(\boldsymbol{x}) $在三种模下的误差. 从表格4和表格5可以观察到, 当$ N\geq20 $时, 精确解$ w(\boldsymbol{x}), u(\boldsymbol{x}) $与逼近解$ w_N(\boldsymbol{x}), u_N(\boldsymbol{x}) $在三种模下的误差都达到了大约$ 10^{-12} $的精度. 由此可以看出我们的算法也是收敛的和高精度的.
本文提出了简支板边界条件下四阶问题基于降阶格式的一种有效的谱Galerkin逼近. 首先, 将原问题化为两个耦合的二阶问题, 建立了相应的弱形式及其离散格式, 证明了弱解和逼近解的存在唯一性及它们之间的误差估计. 其次, 构造了逼近空间中的一组适当的基函数, 推导了离散格式基于张量积的矩阵形式. 另外, 我们还给出一些数值算例, 数值结果表明了算法的有效性和理论结果的正确性. 本文提出的算法可以结合谱元法应用到更一般区域上四阶问题的数值计算.