考虑如下的四阶非线性双曲方程[1-2]
其中$ {\Omega}\subset {\Bbb R}^{2} $是具有Lipschitz连续边界的有界凸多边形区域, $\partial{\Omega} $为$ \Omega $的边界, $T\in(0, +\infty)$, $\gamma $为一正的定值, $X=(x, y)$, $u_{0}\triangleq u_{0}(X), u_{1}\triangleq u_{1}(X) $是充分光滑的已知函数, $f(u) $关于$ u $是全局Lipschitz连续的, 即存在常数$ C>0 $使得
针对四阶方程有限元方法的研究一直受到学者们的关注:一种途径是Galerkin有限元方法[3-11].对协调的有限元而言, 有限元空间的分片多项式具有整体的$ C^{1} $光滑性.为满足这种光滑性要求多项式具有较高次数, 如文献[3]分别利用三角形网格下的Argyris元, Bell元及矩形网格下的Bogner-Fox-Schmit(BFS)元分析四阶椭圆问题时, 其形函数空间分别是五次多项式, 不完全五次多项式以及双三次多项式.文献[4]在矩形网格下利用双三次Hermite型有限元研究了Cahn-Hilliard方程在$ L^{\infty} $和$ W^{1, \infty} $模意义下的最优误差估计和导数逼近的超逼近结果.为了降低形函数空间多项式的次数及减少单元自由度的个数, 人们自然地想到非协调元, 并构造出一些非协调元且已应用到实际中去, 如Morley元[5-7]和Veubeke元[8]是非$ C^{0} $连续的三角形元. Adini元[9-10]是$ C^{0} $连续的非协调元, 但它的法向导数的平均值在跨过单元边界时是并不连续, 其收敛性依赖于矩形网格特殊的几何性质.文献[11]进一步地构造了四阶问题的三个$ C^{0} $非协调元, 证明了其中一个三角形单元是一阶收敛的, 另外两个矩形单元是二阶收敛的.从上述文献可以看出, 无论是协调元还是非协调元都要求广义解具有较高的光滑度以及较高次数的插值多项式(三次及以上)并且一些非协调元对网格剖分限制较苛刻.为克服以上缺点, 另一种途径是利用混合有限元方法, 其思想是将高阶方程进行降阶, 化为低阶方程组成的系统.由于混合元方法具有对空间要求光滑度较低, 并能同时得到原始变量和中间变量数值解等优势, 该方法被广泛应用于高阶发展方程的数值求解[12-16].其中, 文献[12]讨论四阶障碍间题的稳定化混合有限元逼近, 借助网格依赖范数技巧和正则化方法, 给出了基于$ C^{0} $协调有限元空间$ (\vec{W}_{h}, V_{h}) $的混合有限元逼近和相关变量的误差估计.文献[13]研究了四阶抛物型积分-微分方程的混合间断时空有限元格式, 借助椭圆投影算子证明了离散解的稳定性, 存在唯一性和收敛性.文献[14]采用混合有限体积元方法对文献[13]的问题进行了研究, 同样得到了有限元解的稳定性, 存在唯一性和收敛性.文献[15]对四阶抛物方程提出了一个H1-Galerkin混合元格式, 通过引入三个适当的辅助变量, 形成四个一阶的方程组系统, 借助椭圆投影算子, 得到了一维和二维情形下的半离散和全离散格式下的最优误差估计及多维情形下半离散格式最优误差估计, 最后通过数值算例验证了算法的可行性.文献[16]利用双线性元及零阶Raviart-Thomas元对文献[15]的抛物方程建立了半离散和向后欧拉全离散H1-Galerkin混合有限元格式.借助积分恒等式技巧和单元特殊构造, 证明了两个新的重要性质, 从而在半离散及全离散格式下导出了相关变量的最优误差估计以及分别具有$ O(h^{2}) $阶及$ O(h^{2}+\Delta t) $阶的超逼近结果.
文献[17-18]在传统混合元基础上, 对椭圆问题提出了扩展的混合元方法, 该方法能同时逼近未知函数, 未知函数的梯度及扩散通量.这种方法在处理具有复杂边界和小扩散系数问题方面具有优势, 但其混合元空间必须满足所谓离散的LBB条件, 这使得混合元的构造或自由度的选取变得相当复杂.为了避免这个问题, 文献[19-20]对二阶椭圆问题提出了一个新的混合元格式, 与传统的混合元相比具有空间对匹配更容易满足离散的LBB条件, 自由度少且可避免对矢量有限元空间的试探函数进行散度运算等优点.因此, 该格式已被广泛应用在二阶椭圆方程的高精度分析[21]; 线性抛物方程[22], 线性Sobolev方程[23-24]的超收敛分析以及线弹性方程[25]的最优误差估计等.
对四阶抛物问题而言, 文献[26]通过引入两个中间变量$ v=-\Delta u $和$ \vec{p}=-\nabla u$, 将四阶方程转化为三个二阶的方程组成的系统并利用线性三角形元构造了一个扩展的混合元格式.借助投影算子, 在半离散和全离散格式下分别得到了原始变量$ u $和扩散项$ v=-\Delta u $在$ H^{1} $模意义下及流量$ \vec{p}=-\nabla u $在$ L^{2} $模意义的最优误差估计.
本文的主要目的是将文献[19-20]构造混合元空间对的优势与文献[26]扩展的混合元格式相结合, 对一类非线性四阶双曲方程$ (1.1) $利用双线性元空间$ Q_{11} $及其梯度空间$ Q_{01}\times Q_{10} $提出一个低阶的矩形混合元格式.该格式即保持了文献[19-20]混合元空间的优势, 又具有能同时逼近原始变量及多个中间变量的特点.采用文献[28-29]插值算子和投影算子相结合的估计技巧, 利用上述两个单元所诱导的插值算子的高精度分析结果以及对时间$ t $的导数转移技巧, 在半离散格式下导出了原始变量$ u $和中间变量$ v=-\Delta u $在$ H^{1} $模及流量$ \vec{p}=-\nabla u $在$ L^{2} $模意义下具有$ O(h^{2}) $阶的超逼近和超收敛结果.另一方面, 在全离散格式下, 证明了单独利用插值或投影所不能得到的$ u $和$ v $在$ H^{1} $模意义下及$ \vec{p} $在$ L^{2} $模意义下具有$ O(h^{2}+(\Delta t)^{2}) $阶的超逼近和超收敛结果.
本文用$ W^{m, p} $表示通常的Sobolev空间, 其范数和半范分别记为$ \|\cdot\|_{m, p} $和$ |\cdot|_{m, p}$, 特别地, 当$ p=2 $时, $W^{m, p} $简记为$ H^{m}(\Omega)$, 相应的范数和半范简记为$ \|\cdot\|_{m} $和$ |\cdot|_{m}$,
文中出现的$ C $均表示与剖分尺寸无关的常数, 在不同的地方取值可以不同.
为了方便, 不妨设$ \Omega $的边界$ \partial{\Omega} $分别平行于$ x $轴与$ y $轴, $\Gamma_{h} $是$ \Omega $的矩形单元剖分族, 满足正则性假设.对任何$ K\in\Gamma_{h}$, 设其中心为$ (x_K, y_K)$, 四个顶点坐标分别为$ a_{i}(x, y)$(i=1, 2, 3, 4)$.单元$ K $平行于$ x $轴和$ y $轴的边分别为$ l_1, l_3 $及$ l_2, l_4$.边$ l_1=\overline{a_1a_2}$, $l_3=\overline{a_3a_4} $及$ l_2=\overline{a_2a_3}$, $l_4=\overline{a_4a_1} $的边长分别为$ 2h_{x, K}$, $2h_{y, K}$, 记
定义混合有限元空间如下
其中$Q_{ij}={\rm span}\{x^ry^s, 0\leq r\leq i, 0\leq s\leq j\}$.显然, 有限元空间$ M_h $和$ \vec{W}_h $满足$ \nabla M_h\subset \vec{W}_h$.
对于$ u\in H^2(\Omega), \vec{w}\in (H^1(\Omega))^2$, 设$ I_h:H^2(\Omega)\rightarrow M_h $和$ \Pi_h:(H^1(\Omega))^2\rightarrow \vec{W}_h $分别为由$ M_h $和$ \vec{W}_h $上诱导的插值算子, 满足$I_h|_K=I_K, \Pi_h|_K=\Pi_K $及
其中$ \vec{\nu} $是对应边$ l_i$(i=1, 2, 3, 4) $的单位切向量.
文献[27]利用积分恒等式技巧已证明如下结论.
引理2.1 若$u\in H^3(\Omega), \vec{p}\in (H^2(\Omega))^2 $有
若$ u\in H^{3}(\Omega)$, 类似于文献[27]的证明可得如下估计式
其中$ (u, v) $表示$ L^{2} $上的内积.
定义$ R_{h}$: $H_{0}^{1}(\Omega)\rightarrow M_{h} $的Ritz投影, 即对$ \forall u\in H_{0}^{1}(\Omega)$, 满足
文献[28]给出了精确解与Ritz投影之间的误差估计.
引理2.2 若$u\in H^1_{0}(\Omega)\cap H^2(\Omega) $有
根据文献[20, 28]可以得到如下估计结果.
引理2.3 若$(\hat{R}_{h}v, \hat{R}_{h}\vec{p}):[0, T]\rightarrow M_{h}\times\vec{W}_{h} $的定义如下
则有
下面, 我们给出插值与投影之间的高精度结果, 它将在后面的误差分析中起重要作用.
引理2.4 若$u, v\in H^1_{0}(\Omega)\cap H^3(\Omega), \vec{p}\in (H^{2}(\Omega))^{2} $有
证 令
根据(2.6)式可得
在$ (2.11) $式中分别取$ \phi_{h}=\tau, \vec{w}_{h}=\nabla\tau$, 然后两式相减, 根据(2.1)式则有
从而有
在(2.11)式第二式中取$ \vec{w}_{h}=\vec{\theta}$, 根据(2.2)式易得
利用(2.12)式可得
类似于文献[28]的证明方法, 可以得到(2.9)式.引理得证.
令$ \vec{p}=-\nabla u$, $v=\nabla\cdot\vec{p}$, 则问题(1.1)等价于下面问题
问题(3.1)的变分问题为:求$\{u, v, \vec{p}\}:[0, T]\rightarrow H^1_0(\Omega)\times H^1_0(\Omega)\times (L^2(\Omega))^2$, 使得
考虑(3.2)式的半离散格式:求$\{u_{h}, v_{h}, \vec{p}_{h}\}:[0, T]\rightarrow M_h\times M_h\times \vec{W}_h$, 使得
容易验证问题(3.3)存在唯一解.下面, 将在半离散格式下给出混合元解的超逼近和超收敛结果.
定理3.1 设$ \{u, v, \vec{p}\} $和 $ \{u_{h}, v_{h}, \vec{p}_h\} $分别是问题$(3.1) $和$(3.3) $的解, 当$ u, v\in H^3(\Omega)$, $u_{tt}, v_{t}, v_{tt}\in H^2(\Omega)$, $\vec{p}\in (H^2(\Omega))^{2}$, 则有如下超逼近性质
其中
由(3.1)和(3.3)式可得下面的误差方程
在方程(3.7)第一式式中取$ \phi_h=\tau_{t}$, 同时对第二、三式先对$ t $求导数, 然后分别取$ \chi_{h}=\xi_{tt}, $ \vec{w}_{h}=\nabla\xi_{tt}$, 则有
根据(3.8)式可得
下面依次估计$ A_i$(i=1, 2, \cdot\cdot\cdot, 5)$.
利用Schwarz不等式, Young不等式及插值理论得
由于$ f $满足Lipschitz条件, 则有
利用导数转移技巧, 则有
将上述关于$ A_i$(i=1, 2, \cdot\cdot\cdot, 5) $的估计代入(3.9)式可得
将上式两端都乘以$ \frac{2}{\min\{1, \gamma\}}$, 然后从$ 0 $到$ t $积分, 并注意到$ \tau(0)=\xi_{t}(0)=\xi(0)=0$, $ \|\xi\|_{0}^{2}\leq C\int_{0}^{t}\|\xi_{t}\|_{0}^{2}{\rm d}s$, 再根据不等式$ \int_{0}^{t}\int_{0}^{s}\varphi^{2}{\rm d}s{\rm d}\tau\leq C\int_{0}^{t}\varphi^{2}{\rm d}s $ (见文献[30]), 则有
利用Gronwall不等式得
另一方面, 注意到$ \xi(0)=\nabla\xi(0)=0$, 则由$ \|\xi\|_{1}^{2}\leq C\int_{0}^{t}\|\xi_{t}\|_{1}^{2}{\rm d}s $可得
即
在方程(3.7)第三式中取$ \vec{w}_{h}=\vec{\theta}$, 由Schwarz不等式可得
根据(3.11)式则有
根据(2.9)和(3.11)式, 利用三角不等式得
同理, 由(2.10)-(2.11)式及(3.10)和(3.12)式, 利用三角不等式得
定理得证.
为得到整体超收敛, 将相邻的四个单元合并为一个大单元, 采用文献[23]中构造的插值后处理算子$ I^2_{2h} $及 $ \Pi^2_{2h}$
及
结合定理3.1和上述插值后处理算子的性质可得如下整体超收敛结果.
定理3.2 在定理$3.1 $的条件下, 则有
证 根据(3.13)及(3.4)式, 利用三角不等式得
根据(3.13)及(3.5)式, 利用三角不等式, 类似上述证明可得
根据(3.14)及(3.6)式, 利用三角不等式得
注3.1 若直接利用插值, 利用双线性元$ Q_{11} $及$ Q_{01}\times Q_{10} $元的高精度结果, 借助对时间$ t $的导数转移技巧, 当$ u, u_{t}, u_{tt}, v, v_{t}\in H^{3}(\Omega)$, $v_{tt}\in H^{2}(\Omega)$, $\vec{p}\in (H^{2}(\Omega))^{2}$, 则有如下超逼近结果
与定理$ 3.1 $比较, 可以看出这里采用插值和投影相结合的方法降低了$ u_{t}, u_{tt}, v_{t} $的光滑度.
注3.2 可以验证很多著名的非协调元, 如矩形网格下的$ EQ_{1}^{rot} $元[31-36], 正方形网格下的$ Q_{1}^{rot} $元[37]或带约束的旋转$Q_{1}$元[38](等价于矩形网格上的P1非协调元[39]), 由于其相容误差只能估计为
其中$ \|\cdot\|_{h}=\Big(\sum\limits_{K}|\cdot|_{1, K}^{2}\Big)^{1/2}$是$ M_{h} $上的模.所以, 到目前为止还无法得到定理$3.1 $的结果.但是, 在本文定理$3.1 $条件下, 若增加条件$ v_{t}\in H^{3}(\Omega)$, $\vec{p}_{t}\in (H^{2}(\Omega))^{2}$, 利用导数转移技巧也可以得到半离散格式下具有$ O(h^{2}) $阶的超逼近结果.而这里给出的混合元格式的总体自由度只有$ 4NP$ (这里$ NP $是关于$ \Omega $的剖分中所有节点的个数).
为了使时间步长$ \Delta t $和空间剖分参数$ h $都具有二阶精度, 在这一部分, 我们将给出问题的全离散逼近格式, 并进行误差分析.
设$ 0=t_0 < t_1 < \cdots <t_{N-1}<t_N=T $是$ [0, T] $上步长为$ \Delta t=T/N $的剖分, $t_n=n\Delta t, n=0, 1, 2, \cdots, N$, $U^n $表示$ t=t_n=n\Delta t $时$ u(t_n) $在$ M_{h} $中的逼近.为了方便起见, 我们引入下面一些记号
定义(3.1)的全离散逼近格式:求$ \{U^{n}, V^{n}, \vec{P}^{n}\}\in M_{h}\times M_{h}\times \vec{W}_{h}$, 对$\forall\phi_{h}, \chi_{h}\in M_{h}, \vec{w}_{h}\in \vec{W}_{h}$, 使得
其中$ u_{tt}(0)=-(\gamma\Delta^{2}u_{0}-\Delta u_{0}-\Delta u_{1}+f(u_{0})). $
为了进行误差估计, 我们引入如下记号
由(3.1)和(4.1)式, 可得误差方程
其中$ R_{1}^{n}={{\bar{\partial }}_{tt}}{{u}^{n}}-u_{tt}^{n, \frac{1}{4}}=O({{(\Delta t)}^{2}}{{u}_{tttt}}), R_{2}^{n}={{\bar{\partial }}_{t}}{{v}^{n}}-v_{t}^{n, \frac{1}{4}}=O({{(\Delta t)}^{2}}{{v}_{ttt}}).$
定理4.1 设$ \{u(t_{n}), v(t_{n}), \vec{p}(t_{n})\} $和$ \{U^{n}, V^{n}, \vec{P}^{n}\} $分别是问题(3.1)和(4.1)的解, 假设$ u, v\in{L^{\infty}}(0, T;H^{3}(\Omega)), u_{tt}, v_{t}, v_{tt}\in {L^{\infty}(0, T;H^{2}(\Omega))}, u_{tttt}, v_{ttt}\in{L^{\infty}}(0, T;L^{2}(\Omega))$, $ \vec{p}\in{L^{\infty}}(0, T;(H^{2}(\Omega))^{2})$, 则有
证 根据方程(4.2)可得
在(4.6)式第一式中取$ \phi_{h}=\bar{\partial}_{t}\tau^{n}$, 同时在(4.6)式第二、三式中分别取$ \chi_{h}=\bar{\partial}_{tt}\xi^{n}$, $\vec{w}_{h}=\nabla\bar{\partial}_{tt}\xi^{n} $可得
注意到
则(4.7)式的左端可表示为
下面对(4.7)式右端各项进行估计.
根据Schwarz不等式得
由于$ f(\cdot) $满足Lipschitz条件, 则有
根据$ R_{i}^{n}\ (i=1, 2) $的定义可得
将(4.8)式及上述对$ E_{i}\ (i=1, \cdots, 7) $的估计代入到(4.7)式中去, 两端都乘以$ 2\Delta t$, 然后两端从$ 1 $到$ J-1(J\leq N) $求和, 则有
下面, 我们来估计(4.9)式中的$ \|\bar{\partial}_{t}\xi^{\frac{1}{2}}\|_{1}^{2}, \|\nabla\tau^{\frac{1}{2}}\|_{0}^{2}, \|\tau^{\frac{1}{2}}\|_{0}^{2} $这三项.注意到
故
将(4.10)-(4.11)式代入到(4.9)式, 则
在(4.12)式中取适当小的$ \Delta t$, 使得$ \frac{1}{2}-C\Delta t>0$, 然后利用离散的Gronwall引理, 则有
根据(4.2)式的第二、三式, 当$ t=t_{J-\frac{1}{2}} $时分别取$ \chi_{h}=\xi^{J-\frac{1}{2}} $和$ \vec{w}_{h}=\nabla\xi^{J-\frac{1}{2}}$, 然后两式相减得
利用(4.13)式可得
根据(4.2(c))式, 当$ t=t_{J-\frac{1}{2}} $时取$ \vec{w}_{h}=\vec{\theta}^{J-\frac{1}{2}}$, 并利用(4.14)式可得
利用引理2.4及(4.13)-(4.15)式可得
类似于定理3.2的证明过程, 对全离散混合元格式可得下面的超收敛结果.
定理4.2 在定理4.1的条件下有
注4.1 在全离散格式下, 之所以能得到原始变量$ u $和中间变量$ v $及$ \vec{p} $的超逼近和超收敛结果, 是因为采用了插值和投影相结合的技巧:一方面, 如果只利用插值, 则在$(4.7) $式的右端需要估计$(\nabla\bar{\partial}_{t}\eta^{n}, \nabla\bar{\partial}_{tt}\xi^{n})$.根据引理2.1的(2.1)式及逆不等式知, 当$ \Delta t=O(h) $时
这样只能得到具有$ O(h) $阶结果.再者, 若利用非协调元, 为了得到全离散格式下的超逼近超收敛结果, 需要相容误差比插值误差至少高二阶, 尽管矩形网格下的类Wilson元[40], 广义矩形网格下改进的Wilson元[41]和三角形剖分下的$ P_{1}^{\rm mod} $元[42]的相容误差具有$ O(h^{3}) $阶, 但都要求$ u, v\in H^{4}(\Omega)$, 比本文光滑度度高.同时要构造满足匹配关系$ \nabla M_{h}\subset\vec{W}_{h} $混合元空间对也是很困难的, 这也是今后我们进一步要讨论的课题.另一方面, 如果只利用投影算子, 如何构造与投影算子相对应的后处理算子, 也是我们要进一步研究的课题.