众所周知, 非线性发展方程的解通常无法直接用解析式写出来, 或是写出来的表达式非常复杂, 所以利用数值方法给出其近似解就显得尤为重要.而对于有限元方法这一主流方向, 我们常见的线性化BE (Backward-Euler)方法和CN (Crank-Nicolson)方法凭借可以避免在每一个时间层都要求解非线性方程的劣势且不降低计算精度的优势, 成为了该方向的研究热点之一.事实上, 研究一个非线性发展方程的线性化有限元方法总会涉及到一个有限元解$ U_h^n $关于某种模的有界性问题, 由于这些模的先验估计不容易直接得到, 通常的处理技巧就是利用逆不等式.比如在二维的情况下, 考虑有限元解$ \|U_h^n\|_{0, \infty} $有界时的经典做法是

$ \|U_h^n\|_{0, \infty}\leq\|U_h^n-I_hu^n\|_{0, \infty}+\|I_hu^n\|_{0, \infty}\leq Ch^{-1}\|U_h^n-I_hu^n\|_{0}+\|I_hu^n\|_{0, \infty}, $

其中$ u^n $是原始问题的解, $ I_h $是某个插值算子或者投影算子.由于通常有误差估计$ \|U_h^n-I_hu^n\|_{0} = O(h^{m_1}+\tau^{m_2}) $ ($ m_1, m_2 $为某些正数), 要使$ \|U_h^n\|_{0, \infty} $有界就不可避免的要对时间步长$ \tau $有一个限制, 从而导致空间网格参数$ h $与时间步长$ \tau $需要满足某个比值关系(即网格比).在实际计算中, 这样的网格比经常会导致时间步长变的非常小, 从而引起很长的耗时.因此, 怎样甩掉这些限制就成了备受关注的课题.最近, 为了克服这一严重缺陷, 孙伟伟、李步扬、王冀鲁、高华东等学者都在此方面做出了许多有价值的工作.其主要思想(见文献[1])是通过引入一个时间离散方程系统, 并利用其解$ U^n $把误差$ u^n-U_h^n $分裂成两部分—时间误差$ u^n-U^n $和空间误差$ U^n-U_h^n $, 利用时间误差的结果得到关于时间离散方程解的正则性, 再利用空间误差得到有限元解的无网格比有界性

$ \|U_h^n\|_{0, \infty}\leq\|U_h^n-I_hU^n\|_{0, \infty}+\|I_hU^n\|_{0, \infty}\leq Ch^{-1}\|U_h^n-I_hU^n\|_{0}+\|I_hU^n\|_{0, \infty}. $

事实上, 由于空间误差的分析过程中甩掉了经典误差估计中的截断误差项, 只要空间上的误差能写成$ \|U_h^n-I_hU^n\|_{0} = O(h(h^{m_3}+\tau^{m_4})) $ ($ m_3, m_4 $为某些正数)的形式, 网格比即可去掉.随后, 王冀鲁、高华东、司志勇等又把该思想应用于非线性多孔介质流问题^{[2, 3]}, 非线性的Joule Heating方程^[4], 非线性Thermistor方程^{[5, 6]}, 非线性Schrödinger方程^{[7, 8]}和非线性Navier-Stokes方程^[9]等.以上研究都考虑了这些非线性发展方程在协调元下关于无网格比的收敛性, 有许多问题需要进行更深层次的研究.

首先, 为了提高有限元解的逼近精度, 超收敛的思想已成为了一个重要的研究途径.事实上, 在理论分析和实际计算中, 若有好的网格, 有限元解与有限元插值的误差在某种范数的意义下比有限元解与真解的误差要小得多, 即超逼近现象.从上个世纪80年代开始, 以林群院士为代表的众多学者专家都在此方面取得了许多有出色的成果, 所以如何将无网格比收敛的结果推广到无网格比超收敛上去是我们感兴趣的话题.但是, 为了达到超收敛的结果, 如果我们把文献[1-9]中所考虑的区域换成更具一般性的矩形区域(不再满足$ C^2 $的条件), 则由椭圆的正则性可以看到, 引入的时间离散方程解的有界性就很难达到$ H^3 $ -模.因此, 如何在时间离散方程解的有界性较弱的前提下, 探讨无网格比的超逼近结果就显得尤为重要.

其次, 由于非协调元方法在大多数情况下对方程解的正则性要求比较低, 因此人们对非协调元的研究一直保持着较高的热度(见以石钟慈院士为代表的众多学者专家所得到的具有特色的工作).这样就很有必要研究怎样利用自由度少、精度高的低阶非协调单元研究非线性发展方程的无网格比的超收敛性.

再次, 传统的有限元方法对解的光滑度要求都比较高, 这会给实际计算造成很多困难, 因此混合有限元方法受到了高度的关注.事实上, 混合有限元方法的关键性问题在于如何构造出合适的空间对, 使其满足LBB条件, 这其实是不容易做到的.因此构造特别的格式来降低LBB条件的难度成为了一个热点, 比如:文献[10-13]对二阶椭圆问题提出了一种混合元格式, 它具有当两个逼近空间满足一个简单的包含关系时, LBB条件容易满足且能避勉因涉及散度算子带来的麻烦等优点.另一方面, 直接绕开LBB稳定性条件(如最小二乘法、稳定化有限元方法等)也成为了大家另一个关注的方面.事实上, 1998年, Pani在文献[14]中提出了一种称之为$ H^{1} $-Galerkin混合有限元方法.这种方法不需要所选取的混合元空间满足LBB相容性条件, 并被广泛应用在各种方程上.例如, 长波方程^[15], 双曲方程^[16], 带有记忆项的方程^[17], 积分微分方程^[18-20], 抛物方程^[21].因此, 怎样利用$ H^1 $-Galerkin方法得到非线性发展方程无网格比的超收敛结果是值得深思的.

最后, 对于线性化的全离散格式来说, 由于当时时刻的时间层分析需要用到上一时刻时间层的结论, 我们往往会选择数学归纳法进行证明.但是对于每一个时间层的结果到最后都应该由一个统一的系数来控制这个问题显然就不是一件容易的事了.更进一步地, 根据不同非线性问题的具体特点, 针对不同方程的逼近格式, 设计新的高效的有限元数值算法来验证理论分析的正确性也是必须且困难的.

最近, 我们在文献[22-34]中, 在分裂思想的基础上, 博采众家之长, 创新性的把无网格比、高精度分析与非协调单元、新的线性化离散格式等特色和优势有机结合起来, 形成非线性发展方程在全离散格式上无网格比约束的有限元超收敛分析的一套新理论体系.与此同时, 更是尝试着探索、研究一些特殊的方程, 考虑绕过分裂的方法也达到无网格比的超收敛结果.

近期所做的工作, 主要的创新点集中表现在以下几个方面:

(1) 超收敛结果对方程解的光滑性要求比较高, 但构造时间离散辅助问题(即时间离散方程系统)时, 在多边形区域(例如矩形)下, 就无法保证其解较强模的有界性.因此我们利用了一些特殊的、不同以往的技巧, 在其解空间较弱的条件下得到无网格比超收敛的结论; 巧用Taylor展开式对非线性项进行处理, 以保证对时间步长$ \tau $的阶不丢失.

(2) 由于选择的全离散格式是线性化的形式, 在利用数学归纳法分析第$ n $层的结果时需要用到第$ n-1 $层的结论, 我们用一个统一的系数来控制每一个时间层的结果, 这也是其数学归纳法成立的关键所在.

(3) 构造了非线性双曲方程新的二阶格式, 以此得到无网格比超收敛结果.而以往对非线性双曲方程的无网格比研究甚至连收敛性也没有见到报道.

(4) 对一些特殊的非线性发展方程, 抛弃分裂误差思想, 采用一些新的技巧也证明了其无网格比超收敛性.

本文的目的是在前期我们所做的工作的基础上, 挑拣出有特色的创新点给予说明, 以期窥探出对非线性发展方程无网格比超收敛分析的重要方法和思路, 起到抛砖引玉的作用.

非线性抛物方程有着深刻的物理背景, 它的有限元方法也越来越受人们关注.例如:文献[35]针对一般的非线性抛物方程建立了两种线性化的格式, 当$ \tau\leq h $时, 利用线性三角形元得到了$ L^2 $ -模意义下的收敛结果.文献[36]在$ \tau = O(h^{\frac{q}2}), q\geq 4 $限制下利用两层时间离散的方法讨论了非线性抛物方程的最优误差估计.文献[37]采用了一个非线性$ H^1 $投影, 当$ \tau^4 = O(h^{q}), q\leq3 $时, 得到了其解在$ L^2 $ -模和$ H^1 $ -模意义下的最优误差估计.文献[1]利用分裂技巧摆脱了此类限制, 给出了一类称之为Joule Heating的非线性抛物型方程的协调元无网格比收敛性分析.我们看到, 一方面, 一般的非线性抛物方程中的非线性项$ \nabla\cdot(a (u)\nabla u) $的处理对于超收敛的分析是很有挑战的, 特别是在分析空间误差的时候, 怎样处理$ a(u) $才可以在不降阶的情况下使其结果能提出空间网格参数$ h $, 才能在使用逆不等式的时候不产生网格比, 但同时还得能维持数学归纳法所需要要的系数统一性?另一方面, 当非线性抛物方程右端的非线性项是局部Lipschitz连续时, 对有限元解的正则性要求可能会更苛刻, 如何把这些限制考虑进去且得到无网格比超收敛结果是我们想要研究的方向之一.

考虑如下非线性抛物方程

$ \begin{equation} \begin{cases} u_t-\nabla\cdot(a(u)\nabla u) = f(u), & (X, t)\in\Omega\times(0, T], \\ u = 0, & (X, t)\in\ \partial \Omega\times(0, T], \\ u(X, 0) = u_{0}(X), & X\in\Omega, \end{cases} \end{equation} $

(2.1)

其中$ \Omega\subset \mathbb{R}^{2} $是一个矩形, 其边界为$ \partial \Omega $, $ 0<T<\infty $, $ X = (x, y) $, $ u_{0}(X) $是已知光滑函数. $ a(u) $关于$ u $是二阶可导连续的函数, 其中$ 0<a_0\leq a(u)\leq a_1 $, $ a_0, a_1 $是某些正数.

利用低阶非协调单元$ EQ_1^{rot} $(参见文献[21, 38-40]), 对(2.1)式开展了无网格比超收敛性质的研究.在矩形区域下, 引进了一个时间离散方程, 证明了时间离散方程解的$ H^2 $ -模有界.绕过时间离散方程解较弱的正则性, 利用Taylor展开, 在不降低时间方向阶的前提下, 得到其CN格式的无网格比超逼近结果.另一方面, 限制(2.1)式右端项为仅满足局部Lipschitzt连续条件, 利用数学归纳法, 巧妙的使用几个不等式, 在每一层都得到数值解$ L^{\infty} $ -模有界的前提下, 保证了结果系数的统一性, 采用双线性元(参见文献[13, 41]), 得到其在BE格式下无网格比超逼近结果.

首先, 设$ f(u) $是$ \Omega $上整体Lipschitz连续的函数, $ \Omega $是一个四条边都平行于坐标轴的矩形, $ \Gamma_{h} $是一个拟一致正则矩形剖分.对于给定的$ K\in\Gamma_{h} $, 令其四个顶点和四条边分别为$ a_{i}, i = 1\sim4 $和$ l_{i} = \overline{a_{i}a_{i+1}}, i = 1\sim4\pmod{4} $.记$ h = \max\limits_{K\in \Gamma_h}{\mathrm{diam} K} $.定义非协调有限元空间$ EQ_1^{rot} $:

$ V_{h} = \{v_h;v_h|_{K}\in \operatorname{span}\{1, x, y, x^2, y^2\}, \int_F[v_h]ds = 0, F\subset\partial K, \forall K\in \Gamma_{h}\}, $

其中$ [v_h] $表示$ v_h $跨过单元边界$ F $的跳度, 而当$ F\subset\partial\Omega $时, $ [v_h] = v_h $.令$ I_{h}:H^1(\Omega)\rightarrow V_h $为相对应的插值算子, 且$ I_{h} = I_h|_{K} $满足

$ \begin{align} \label{chapter3.1.1} \int_{l_i}(I_{h}u-u)ds = 0, i = 1\sim4, \int_{K}(I_{h}u-u)dX = 0. \end{align} $

则文献[21, 38-40]证明了下面重要引理.

引理1   若$ u\in H^{2}(\Omega)\cap H^{1}_0(\Omega) $, 则对于任意的$ v_h\in V_h $, 有

$ \begin{eqnarray} &&(\nabla_h(u-I_hu), \nabla_h v_h)_h = 0, \end{eqnarray} $

(2.2)

$ \begin{eqnarray} &&\left|\sum\limits_{K} \int_{\partial K}\frac{\partial u}{\partial n}v_hds\right| = O(h)|u|_2\|v_h\|_h. \end{eqnarray} $

(2.3)

进一步地, 若$ u\in H^{3}(\Omega)\cap H^{1}_0(\Omega) $, 则有

$ \begin{equation} \; \left|\sum\limits_{K} \int_{\partial K}\frac{\partial u}{\partial n}v_hds\right| = O(h^2)|u|_3\|v_h\|_h. \end{equation} $

(2.4)

这里$ \nabla_h $表示分片梯度, $ (\star, \star)_h = \sum\limits_{K} \int_{K}\star\cdot\star dX $且$ \|\cdot\|_h = (\sum\limits_{K}|\cdot|_{1, K}^2)^{\frac12} $是$ V_h $上的一个能量模.文献[40]证明了对于任意正整数$ m, v_h\in V_h $,

$ \begin{equation} \|v_h\|_{0, 2m}\leq \|v_h\|_{h}. \end{equation} $

(2.5)

设$ \{t_n:t_n = n\tau;0\leq n\leq N\} $是$ [0, T] $上的一个等距剖分, 时间步长为$ \tau = T/N $, 设$ t^{n-\frac12} = (n-\frac12)\tau $, 且$ u(X, t_n) = u^n $, 若$ \{\sigma^n\}_{n = 0}^{N} $为一列函数.记

$ \begin{eqnarray*} &&\bar{\partial}_t\sigma^n = \frac{\sigma^n-\sigma^{n-1}}{\tau}, \; \; \bar{\sigma}^n = \frac{\sigma^n+\sigma^{n-1}}{2}, n = 1, 2, \cdots, N, \\ &&\hat{\sigma}^n = \frac{3}{2}\sigma^{n-1}-\frac{1}{2}\sigma^{n-2}, n = 2, \cdots, N. \end{eqnarray*} $

利用这些记号, 考虑(2.1)式的线性化Galerkin有限元逼近:寻找$ U_h^n\in V_h $, 使得对于任意的$ v_h\in V_{h0} $,

$ \begin{equation} (\bar{\partial}_tU^n_h, v_h)+(a(\hat{U}_h^n)\nabla_h \bar{U}_h^n, \nabla_h v_h)_h = (f(\hat{U}_h^n), v_h), n = 2, \cdots, N. \end{equation} $

(2.6)

类似于文献[35], 由以下方程求解$ U_h^1 $:

$ \begin{eqnarray} &&(\frac{U^{1, 0}_h-U^0_h}{\tau}, v_h)+(a(U^0_h)\frac{\nabla_h U^{1, 0}_h+\nabla_h U^0_h}{2}, \nabla_h v_h)_h = (f(U_h^0), v_h), \end{eqnarray} $

(2.7)

$ \begin{eqnarray} &&(\frac{U^{1}_h-U^0_h}{\tau}, v_h)+(a(\frac{U^{1, 0}_h+U^0_h}{2})\frac{\nabla_h U^{1}_h+\nabla_h U^0_h}{2}, \nabla_h v_h)_h = (f(\frac{U^{1, 0}_h+U^0_h}{2}), v_h), \end{eqnarray} $

(2.8)

其中$ v_h\in V_h $, $ U_h^0 = I_hu_0 $.可以看到先利用(2.7)式计算出$ U_h^{1, 0} $, 再利用(2.8)式得到$ U_h^1 $.由于线性化后, (2.6)–(2.8)式是一个线性系统, 其解的存在唯一性是显然的.下面分步骤地阐述有创新性的重要过程.

第一步  建立一个时间离散系统, 当$ n> 1 $时求$ U^n $满足

$ \begin{equation} \begin{cases} \bar{\partial}_tU^n-\nabla\cdot(a(\hat{U}^n)\nabla \bar{U}^n) = f(\hat{U}^n), & X\in\Omega, \\ U^n = 0, & X\in\partial\Omega, \\ U^0(X) = u_0(X), &X\in\Omega. \end{cases} \end{equation} $

(2.9)

当$ n = 1 $时, 利用以下式子计算$ U^1 $:

$ \begin{equation} \frac{U^{1, 0}-U^0}{\tau}-\nabla\cdot(a(U^0)\frac{\nabla U^{1, 0}+\nabla U^0}{2}) = f(U^0) \end{equation} $

(2.10)

和

$ \begin{equation} \frac{U^{1}-U^0}{\tau}-\nabla\cdot(a(\frac{U^{1, 0}+U^0}{2})\frac{\nabla U^{1}+\nabla U^0}{2}) = f(\frac{U^{1, 0}+U^0}{2}), \end{equation} $

(2.11)

其中$ U^{1, 0}(X)|_{\partial\Omega} = 0, U^{1}(X)|_{\partial\Omega} = 0. $

令$ e^{1, 0}\triangleq u^1-U^{1, 0}, e^n\triangleq u^n-U^n (n = 0, 1, 2, \cdots, N) $.通过分析时间误差, 给出$ U^{1, 0}, U^n(n = 0, 1, 2, \cdots, N) $的正则性.设$ u $和$ U^m (m = 0, 1, 2, \cdots, N) $分别为(2.1)和(2.9)–(2.11)式的解, $ u\in L^{\infty}(0, T;H^{3}(\Omega)), u_t, u_{tt}\in L^{\infty}(0, T;H^{2}(\Omega)), u_{ttt}\in L^{\infty}(0, T;L^{2}(\Omega)) $, 则对于$ m = 1, \cdots, N, $存在$ \tau_0>0 $, 使得当$ \tau\leq\tau_0 $时, 有

$ \begin{equation} \|e^{1, 0}\|_0+\tau^{\frac12}\|e^{1, 0}\|_1+\tau\|e^{1, 0}\|_2+\|e^m\|_1+\tau\|e^{m}\|_2\leq C_0\tau^2 \end{equation} $

(2.12)

和

$ \begin{equation} \|\frac{U^{1, 0}-U^0}{\tau}\|_2+\|\bar{\partial}_t U^m\|_2\leq C_0, \end{equation} $

(2.13)

其中$ C_0 $是一个与$ m, h $和$ \tau $无关的正数.

此时, 注意到由于$ \Omega $是矩形, 其边界不属于$ C^{1} $, 那么就不容易得到$ U^n $的$ H^3 $ -模有界性.因此随后的无网格比超收敛分析需要利用新的方法得到.

第二步   讨论空间误差$ \tau\sum\limits_{i = 1}^n\|\bar{\partial}_t(I_hU^i-U_h^i)\|_0^2 $, 也为最终无网格比超逼近结果$ \|I_hu^n-U_h^n\|_{h} = O(h^2+\tau^2) $做好准备.给出记号

$ \begin{eqnarray} &&U^{1, 0}-U^{1, 0}_h = (U^{1, 0}-I_hU^{1, 0})+(I_hU^{1, 0}-U^{1, 0}_h)\triangleq r^{1, 0}+\theta^{1, 0}, \end{eqnarray} $

(2.14)

$ \begin{eqnarray} &&U^i-U^i_h = (U^i-I_hU^i)+(I_hU^i-U^i_h)\triangleq r^i+\theta^i, i = 0, 1, 2, \cdots, n. \end{eqnarray} $

(2.15)

令$ U^m $和$ U_h^m $分别为(2.9)–(2.11)和(2.6)–(2.8)式的解, 其中$ m = 1, 2, \cdots, N $, 在前面所做工作的前提下, 当$ \tau $充分小时, 有

$ \begin{equation} \tau\sum\limits_{i = 1}^m\|\bar{\partial}_t\theta^i\|_0^2 = O(h^2). \end{equation} $

(2.16)

第三步   令

$ \begin{equation} I_hu^{1}-U^{1, 0}_h\triangleq\xi^{1, 0}, u^i-U^i_h = (u^i-I_hu^i)+(I_hu^i-U^i_h)\triangleq\eta^i+\xi^i, i = 0, 1, 2, \cdots, n. \end{equation} $

(2.17)

设$ u^n $和$ U_h^n $分别为(2.1)和(2.6)–(2.8)式的解, 则对$ n = 1, 2, \cdots, N $, 有

$ \begin{equation} \|\xi^n\|_h = O(h^2+\tau^2). \end{equation} $

(2.18)

在这个过程中注意到, 将$ \tau $从内积的一端转向另一端的恒等变化

$ (a^n, \frac{b^n-b^{n-1}}{\tau}) = \frac{1}{\tau}((a^n, b^n)-(a^{n-1}, b^{n-1}))-(\frac{a^n-a^{n-1}}{\tau}, b^{n-1}), $

化简过之后需要估计误差

$ \|\frac{(a(\hat{u}^n)-a(\hat{u}^{n-1}))-(a(\hat{U}_h^n)-a(\hat{U}_h^{n-1}))}{\tau}\|_0. $

如果按照传统的方法, 则有

$ \|\frac{(a(\hat{u}^n)-a(\hat{u}^{n-1}))-(a(\hat{U}_h^n)-a(\hat{U}_h^{n-1}))}{\tau}\|_0\leq C\|\bar{\partial}_t\hat{u}^n\|_0 +C\|\bar{\partial}_t\hat{U}_h^n\|_0, $

这样最终的结果会降一阶.但是若利用Taylor展开则有

$ \begin{align} \ &\|\frac{(a(\hat{u}^n)-a(\hat{u}^{n-1}))-(a(\hat{U}_h^n)-a(\hat{U}_h^{n-1}))}{\tau}\|_0\\ = \ &\|\frac{(a_u(\hat{u}^{n-1})(\hat{u}^n-\hat{u}^{n-1})+\frac12a_{uu}(\mu_1^n)(\hat{u}^n-\hat{u}^{n-1})^2)}{\tau} \\ \ &-\frac{(a_u(\hat{U}_h^{n-1})(\hat{U}_h^n-\hat{U}_h^{n-1})+\frac12a_{uu}(\mu_2^n)(\hat{U}_h^n-\hat{U}_h^{n-1})^2)}{\tau}\|_0\\ \leq \ & Ch^2+C\tau^2+C\|\bar{\partial}_t(\hat{U}_h^n-I_h\hat{U}^{n})\|_0+C\|\hat{U}_h^n-I_h\hat{U}^{n}\|_{0}+C\|\hat{U}_h^{n-1}-I_h\hat{U}^{n-1}\|_{0}, \end{align} $

其中

$ \mu_1^n = \hat{u}^{n-1}+\lambda_1(\hat{u}^n-\hat{u}^{n-1}), \mu_2^n = \hat{U}_h^{n-1}+\lambda_2(\hat{U}_h^n-\hat{U}_h^{n-1}), \; \; 0<\lambda_1, \lambda_2<1. $

这样就可以保持到想要的结果.

在这里我们还指出本节的结果对于正方形网格上的非协调$ Q^{rot}_1 $单元(参见文献[42]), 矩形的带约束的非协调$ Q^{rot}_1 $单元(参见文献[43])和$ P_1 $ -非协调四边形单元(参见文献[44])也成立, 因为这些单元都满足引理1.

限制$ f(u) $为局部Lipschitz连续的, 利用双线性协调单元可研究(2.1)式的无网格比超逼近性质(区域及剖分如前面一样).定义其有限元空间$ V_{h0} $为

$ V_{h0} = \{v_h;v_h|_{K_{h}}\in \operatorname{span}\{1, x, y, xy\}, v_{h}|_{\partial\Omega} = 0\}, $

其中$ I_{h}:H^2(\Omega)\rightarrow V_{h0} $是相对应的插值算子, 且对于以上双线性元有以下高精度结果^[41].

引理2   若$ u\in H^{3}(\Omega)\cap H^{1}_0(\Omega) $, 则对于任意的$ v_h\in V_{h0} $, 有

$ \begin{equation} \label{chapter3.1.332} (\nabla(u-I_hu), \nabla v_h) = O(h^2)|u|_3\|v_h\|_1.\nonumber \end{equation} $

考虑(2.1)式的线性化有限元逼近:寻找$ U_h^n\in V_{h0} $, 使得

$ \begin{equation} (\bar{\partial}_tU^n_h, v_h)+(a(U_h^{n-1})\nabla U_h^n, \nabla v_h) = (f(U_h^{n-1}), v_h), n = 1, 2, \cdots, N, \; \forall\; v_h\in V_{h0}, \end{equation} $

(2.19)

其中$ U_h^0 = I_hu_0 $, 显然(2.19)式在每一个时间层只需要解决一个线性问题.

引入时间离散方程:当$ n\geq 1 $时求$ U^n $满足

$ \begin{equation} \begin{cases} \bar{\partial}_tU^n-\nabla\cdot(a(U^{n-1})\nabla U^n) = f(U^{n-1}), & X\in\Omega, \\ U^n = 0, & X\in\partial\Omega, \\ U^0(X) = u_0(X), &X\in\Omega. \end{cases} \end{equation} $

(2.20)

接下来, 我们分别就时间误差和空间误差中的新技巧予以说明.

时间误差   记$ e^n\triangleq u^n-U^n (n = 0, 1, 2, \cdots, N) $, 则有

$ \begin{equation} (\tau\sum\limits_{i = 1}^{m}\|e^{i}\|_2^2)^{\frac{1}{2}}+\tau^{\frac{1}{2}}\|e^{m}\|_2\leq C_0\tau. \end{equation} $

(2.21)

注意到在第$ n-1 $层中归纳假设(2.21)式成立后, 进一步地得到$ \|e^m\|_2\leq\tau^{\frac{1}{4}} $$ (m\leq n-1) $是必需的.主要表现在以下两个方面.

1.在第$ n $层估计中, 误差方程左端有$ a_0\|\Delta e^{n}\|_0^2 $, 右端部分在$ \|e^m\|_2\leq\tau^{\frac{1}{4}} $$ (m\leq n-1) $的前提下有估计项

$ \begin{align} (a_u(U^{n-1})\nabla e^{n-1}\nabla e^n, \Delta e^{n})\leq C\|e^{n-1}\|_{2}\|e^{n}\|_{2}\|\Delta e^{n}\|_0 \leq C\tau^{\frac{1}{4}}\|\Delta e^{n}\|_0^2 \end{align} $

和

$ \begin{align} ((a(u^{n-1})-a(U^{n-1}))\Delta e^n, \Delta e^{n}) \leq C\|e^{n-1}\|_2\|\Delta e^n\|_0^2 \leq C\tau^{\frac{1}{4}}\|\Delta e^{n}\|_0^2. \end{align} $

可以看到, 此时当$ \tau $充分小时, 在第$ n $层误差方程的右端才可以去掉$ \|\Delta e^n\|_0^2 $.

2.由于$ f $的局部Lipschitz连续, 要想估计误差方程右端项$ (f(u^{n-1})-f(U^{n-1}), \Delta e^{n}), $则必须得有$ \|e^{n-1}\|_{0, \infty}\leq C\|e^{n-1}\|_{2}\leq C\tau^{\frac14} $做前提, 这样就保证了

$ (f(u^{n-1})-f(U^{n-1}), \Delta e^{n})\leq C\|\nabla e^{n-1}\|_0^2+\frac{a_0}6\|\Delta e^n\|_0^2. $

空间误差   注意到, 由于$ f $的局部Lipschitz连续, 在估计下面误差时,

$ (f(U^{n-1})-f(U_h^{n-1}), \theta^n), $

利用前面的结论以及协调元的性质, 有

$ \begin{align} \|U^{n-1}-U_h^{n-1}\|_{0, \infty} \ &\leq\|e^{n-1}\|_{2}+\|u^{n-1}-I_h u^{n-1}\|_{0, \infty}+Ch^{-1}\|\theta^{n-1}\|_{0}\\ \ & \leq Ch+C\tau^{\frac12}. \end{align} $

(2.22)

这样可以估计得到

$ \begin{align} (f(U^{n-1})-f(U_h^{n-1}), \theta^n)\leq C h^4\|U^{n-1}\|_2+C\|\theta^{n-1}\|_0^2. \end{align} $

注1   更进一步地, 我们在文献[34]中将文献[10-12]中的混合元和无网格比的思想有机的结合起来, 利用分裂内积等思想, 得到了(2.1)式关于原始变量$ u $的$ H^1 $ -模和$ \vec{p} = \nabla u $的$ L^2 $ -模的无网格比超收敛结果.

对于非线性Schrödinger方程来说, 文献[7, 8]利用分裂技巧得到了协调元的无网格比收敛结果.但是怎样能够使得其无网格比的超收敛结果成立呢？首先, 要使得时间误差$ H^2 $ -模的阶比文献[7, 8]有所提高, 这样才能为下面的无网格比超逼近结果做铺垫.其次, 可以看到, 文献[7, 8]是利用投影算子来得到的无网格比的收敛结果的.事实上, 对我们而言, 仅用投影算子是不能直接利用插值后处理技巧得到整体超收敛的.所以从这个角度考虑, 利用插值算子更有优势.但是, 若仅仅考虑插值算子得到超逼近的结果, 则在$ t_n $时刻对时间离散方程解$ U^n $的正则性要求就比较高了, 而在矩形网格下, 要想得到$ U^n $的$ H^3 $ -模以上的有界性并不是那么容易的事情.既然单独利用插值算子或者投影算子都不能得到令人满意的效果, 那把二者结合在一起是否会有更好的结果?

考虑如下Schrödinger方程

$ \begin{equation} \begin{cases} iu_t+\Delta u+ f(|u|^2)u = 0, & (X, t)\in\Omega\times(0, T], \\ u = 0, & (X, t)\in\ \partial \Omega\times(0, T], \\ u(X, 0) = u_{0}(X), & X\in\Omega, \end{cases} \end{equation} $

(3.1)

其中$ \Omega $同(2.1)式, $ i $是虚数单位, $ u_0(X) $是已知复值函数.另外, $ f(s) $是一个实值函数, 且关于$ s $是二阶可导连续的.

在克服非线性Schrödinger方程由于虚数单位$ i $所带来困难的同时, 需要采用新的技巧, 得到每一层数值解的$ L^{\infty} $ -模有界性质, 保证其每一层解的存在唯一性, 并导出其无网格比超逼近性质.一方面, 我们利用新的技巧得到了比文献[8]更高阶的时间误差, 从而导出了结果更好的时间离散方程解的正则性, 这也为下面无网格比的超逼近奠定了基础.另一方面, 在证明过程中我们引入了经典的Ritz投影算子, 避开了对引入的时间离散方程解的正则性要求过高的麻烦, 得到了合适的空间误差结果.同时, 结合插值算子和Ritz投影算子相结合的思想达到了超逼近的结果.最后, 利用文献[41]中的插值后处理算子得到了整体超收敛性质.

选用上一部分的区域, 剖分和双线性元单元, 且仍定义其有限元空间为$ V_{h0} $.令$ R_h: H_0^1\rightarrow V_{h0} $是定义在$ V_{h0} $上的相对应的Ritz投影算子

$ \begin{align} (\nabla(u-R_hu), \nabla v_h) = 0, \forall v_h\in V_{h0}, \end{align} $

(3.2)

有

$ \begin{align} \|\nabla R_hu\|_0\leq C\|\nabla u\|_{0} \end{align} $

(3.3)

且

$ \begin{align} \|u-R_hu\|_0\leq Ch^s\|u\|_{s}, s = 1, 2, \forall u\in H^1_0(\Omega)\cap H^2(\Omega). \end{align} $

(3.4)

更进一步地, 当$ u\in H^{3}(\Omega) $, 又由文献[13]可知

$ \begin{align} \|I_hu-R_hu\|_1 = O(h^2)\|u\|_3, \end{align} $

(3.5)

其中$ I_h $是定义在$ V_{h0} $上相对应的插值算子.下面我们仍采用了分裂技巧, 分别给出时间误差和空间误差上分析时所遇到的困难和解决方法.

时间误差 对于CN格式   将误差方程相邻两层相减, 则

$ \begin{align} i\tilde{\partial}_te^n+\Delta \tilde{e}^n+P_1^n\triangleq R_1^n+R_2^n+R_3^n \end{align} $

(3.6)

就变成

$ \begin{align} i\tilde{\partial}_{tt}e^n+\tilde{\partial}_{t}\Delta \bar{e}^n+\tilde{\partial}_{t}P_1^n = \tilde{\partial}_{t}R_1^n+\tilde{\partial}_{t}R_2^n+\tilde{\partial}_{t}R_3^n. \end{align} $

(3.7)

至此, 得到的结果$ \|e^n\|_2 = O(\tau^2) $相关估计, 可以比文献[8]的结果高二分之一阶, 也就是这样的结果导出了$ \|\tilde{\partial}_{tt}U^m\|_2\leq C_0 $, 为后面的无网格比超逼近结果奠定了基础.

对于BE格式  得到结果$ \|e^n\|_2 = O(\tau) $比文献[7]中的结果阶高二分之一阶, 这也导出了$ \|\tilde{\partial}_{t}U^m\|_2\leq C_0 $, 也为下面的空间误差做出了铺垫.

空间误差    1.对误差方程相邻两层相减, 对比文献[45]中的结果, 可以看到不需要$ \|\tilde{\partial}_tU_h^n\|_{0, \infty} $的有界性也得到了无网格比高精度的结果, 这就改进了已有结论.

2.使用插值算子和投影算子相结合的思想:若仅使用插值算子, 为了得到高精度结果, 避免不了利用高精度结果$ (\nabla (u^n-I_hu^n), \nabla v_h) = O(h^2)\|u^n\|_3\|v_h\|_1 $或者$ (\nabla (u^n-I_hu^n), \nabla v_h) = O(h^2)\|u^n\|_4\|v_h\|_0 $, 则对$ U^n $和$ u^n $的正则性要求过于苛刻.然而, 在$ \Omega $为一个矩形的前提下, 目前只能得到$ \|U^n\|_2 $的有界性, 此时选用投影算子$ R_h $是合适的.另一方面, 若仅仅使用投影算子$ R_h $, 则不能构造相应于$ R_h $的插值后处理算子, 也就不能得到整体超收敛结果了.

Sobolev方程起源于流体通过裂隙岩石的流动、二阶流体的热力学剪切和粘土的固结等物理现象.到目前为止, 已有很多文献研究了它的数值方法.例如:文献[46]得到了当$ \tau = O(h^{d/3})\; (d\leq3) $时, 在三种情况下关于$ H^1 $ -模最优误差估计结果.而文献[47]利用混合有限元方法, 在条件$ \tau = O(h) $下得到了最优误差估计.文献[38]控制条件为$ \tau = O(h^{1+\varepsilon})(\varepsilon>0) $时分别利用协调有限元和非协调有限元讨论了其特征有限元方法, 也得到了$ H^1(\Omega) $ -模和$ L^2(\Omega) $ -模的最优误差估计.

大家都知道, $ H^1 $-Galerkin方法是一个不需要满足LBB条件的混合有限元方法, 加上一些技巧的应用, 还可得到关于流量$ \vec{p} = \nabla u $散度模的误差估计.但是由于$ H^1 $-Galerkin混合有限元方法需要的时间离散方程解的正则性较高, 在矩形区域下不容易得到, 所以对于一般的诸如非线性抛物方程利用上述分裂技巧直接处理暂时还不能去掉$ h $和$ \tau $的比值.非线性Sobolev方程有着其自身的特点, 它比非线性抛物方程多了一个非线性的导数项, 正是多了这一项, 使得我们考虑在分析的时候可以不用以上的分裂法就得到无网格比超收敛结果.因此我们通过与前面不同的分析, 不引进时间离散方程, 即在不必考虑所谓的时间误差的前提下, 避免由于时间离散方程解的正则性达不到相应的要求而带来的麻烦, 给出了无网格比的超逼近结论.

考虑如下非线性Sobolev方程:

$ \begin{equation} \begin{cases} u_t-\nabla\cdot(a(u)\nabla u_t+b(u)\nabla u) = f(X, t), & (X, t)\in\Omega\times(0, T], \\ u = 0, & (X, t)\in\ \partial \Omega\times(0, T], \\ u(X, 0) = u_{0}(X), & X\in\Omega, \end{cases} \end{equation} $

(4.1)

这里对于正数$ b_{1} $, 有$ |b(u)|\leq b_{1} $, 其余同(2.1)式中的假设.

首次尝试选择非协调元单元对($ EQ_{1}^{rot} $及零阶Raviart-Thomas (RT)单元)构造$ H^{1} $-Galerkin混合有限元格式去解决(4.1)式的无网格比超逼近问题.给出一个重要引理, 完全不同于文献[1-9]的思路, 先估计$ \|\bar{\partial}_t\nabla\xi^n\|_0 $和$ \|\nabla\xi^1\|_0 $的结果, 再通过不等式$ \|\nabla \xi^n\|_0\leq C\tau\sum\limits_{i = 1}^n\|\bar{\partial}_t\nabla\xi^i\|_0 $给出$ \|\nabla\xi^n\|_0 $超收敛性, 其中$ \xi^n = I_hu^n-U_h^n $, $ I_h $是相对应的插值算子.

采用非协调$ EQ_1^{rot} $元, 假定所有符号同第二节.设零阶RT单元的空间$ \vec{W}_{h} $定义为

$ \vec{W}_{h} = \{\vec{w}_h = (w_h^{1}, w_h^{2})\in H(\operatorname{div};\Omega);\vec{w}_h|_{K}\in \{1, x\}\times \{1, y\}, \forall K\in \Gamma_{h}\}. $

对于$ \vec{q} = (\vec{q}_{1}, \vec{q}_{2})\in(H^{1}(\Omega))^{2} $, 定义相对应的插值算子$ \Pi_{h} $为

$ \begin{align} \ &\Pi_{h}:\vec{q}\in(H^{1}(\Omega))^{2}\rightarrow\Pi_{h} \vec{q}\in\vec{W}_{h}, \\ \ &\Pi_{h}|_{K} = \Pi_{K} = (\Pi^1_{K}, \Pi^2_{K}), \int_{l_{i}}(\vec{q}-\Pi_{K}\vec{q})\cdot\vec{n}_ids = 0, \end{align} $

(4.2)

其中$ \vec{n}_i $是单元边上$ l_{i} $的外法线.

令$ \vec{q} = a(u)\nabla u_t+b(u)\nabla u $, 相对应的弱形式是:寻找$ \{u, \vec{q}\}:[0, T]\to H_{0}^{1}(\Omega)\times H(\operatorname{div};\Omega) $, 使得

$ \begin{equation} \begin{cases} (\vec{q}, \nabla v) = (a(u)\nabla u_t, \nabla v)+(b(u)\nabla u, \nabla v), & (X, t)\in\Omega\times(0, T], \\ (\alpha(u)\vec{q}, \vec{w})+(\nabla\cdot \vec{q}, \nabla\cdot\vec{w})-(\beta(u)\nabla u, \vec{w}) = -(f(X, t), \nabla\cdot\vec{w}), & (X, t)\in\Omega\times(0, T], \\ u(X, 0) = u_{0}(X), & X\in\Omega, \end{cases} \end{equation} $

(4.3)

其中$ \alpha(u) = \frac{1}{a(u)}, \; \beta(u) = \frac{b(u)}{a(u)}. $

给出线性化的CN全离散格式:寻找$ \{U^n_{h}, \vec{Q}^n_{h}\}:[0, T]\to V_{h}\times\vec W_{h} $, 使得$ n\ge 2 $时,

$ \begin{equation} \begin{cases} (a(\hat{U}^n_h)\bar{\partial}_t \nabla_h U^n_h, \nabla_h v_h)_h+(b(\hat{U}_h^n)\nabla_h \bar{U}_h^n, \nabla_h v_h)_h\\ = (\bar{\vec{Q}}^n_h, \nabla_h v_h)_h, & \forall v_h\in V_{h}, \\ (\alpha(\hat{U}^n_h)\bar{\vec{Q}}_h^n, \vec{w}_h)+(\nabla\cdot\bar{\vec{Q}}_h^n, \nabla\cdot\vec{w}_h)-(\beta(\hat{U}^n_h)\nabla_h \bar{U}_h^n, \vec{w}_h)_h\\ = -(f(X, t^{n-\frac12}), \nabla\cdot\vec{w}_h), & \forall \vec{w}_h\in \vec{W}_h. \end{cases} \end{equation} $

(4.4)

当$ n = 1 $时,

$ \begin{align} \begin{cases} \ &(a(U^0_h)\frac{\nabla_h U^{1, 0}_h-\nabla_h U^{0}_h}{\tau}, \nabla_h v_h)+(b(U_h^0)\frac{\nabla_h U_h^{1, 0}+\nabla_h U_h^{0}}{2}, \nabla_h v_h)\nonumber\\ \ & = (\frac{\vec{Q}_h^{1, 0}+\vec{Q}_h^{0}}{2}, \nabla_h v_h), \forall v_h\in V_{h}, \nonumber\\ \ & (\alpha(U^0_h)\frac{\vec{Q}_h^{1, 0}+\vec{Q}_h^{0}}{2}, \vec{w}_h)+(\frac{\nabla\cdot\vec{Q}_h^{1, 0}+\nabla\cdot\vec{Q}_h^0}{2}, \nabla\cdot\vec{w}_h) -(\beta(U^0_h) \frac{\nabla_h U_h^{1, 0}+\nabla_h U_h^{0}}{2}, \vec{w}_h)\nonumber\\ \ & = (f(X, t^{n-\frac12}), \nabla\cdot\vec{w}_h), \forall \vec{w}_h\in \vec{W}_h \nonumber \end{cases} \end{align} $

和

$ \begin{align} \begin{cases} \ &(a(\frac{U^{1, 0}_h+U^{0}_h}{2})\frac{\nabla_h U^1_h-\nabla_h U^0_h}{\tau}, \nabla_h v_h)+(b(\frac{U^{1, 0}_h+U^{0}_h}{2})\frac{\nabla_h U^1_h+\nabla_h U^0_h}{2}, \nabla_h v_h)\nonumber\\ \ & = (\frac{\vec{Q}_h^{1}+\vec{Q}_h^0}{2}, \nabla_h v_h), \forall v_h\in V_{h}, \nonumber\\ \ & (\alpha(\frac{U^{1, 0}_h+U^{0}_h}{2})\frac{\vec{Q}_h^{1}+\vec{Q}_h^0}{2}, \vec{w}_h)+(\frac{\nabla\cdot\vec{Q}_h^{1}+\nabla\cdot\vec{Q}_h^0}{2} , \nabla\cdot\vec{w}_h) -(\beta(\frac{U^{1, 0}_h+U^{0}_h}{2}) \frac{\nabla_h U_h^1+\nabla_h U_h^0}{2}, \vec{w}_h)\nonumber\\ \ & = (f(X, t^{n-\frac12}), \nabla\cdot\vec{w}_h), \forall \vec{w}_h\in \vec{W}_h, \nonumber \end{cases} \end{align} $

其中$ U_h^0 = I_hu_0 $和$ \vec{Q}_h^0 = \Pi_h\vec{q}^0 $, $ \vec{q}^0 = a(u_0)\nabla u^0_t+b(u_0)\nabla u_0 $.

下面先给出一个新的引理.

引理3   对于任意的$ u\in H^{1}_0(\Omega) $, 则有

$ \begin{equation} (\nabla_h (u-I_hu), \vec{w}_h)_h = 0, \; \; \forall \vec{w}_{h}\in\vec{W}_{h}, \end{equation} $

(4.5)

这里仅给出该引理成立的关键点.事实上, 由$ \vec{W}_h $的定义可知, $ \nabla\cdot \vec{w}_h $在$ K $上为常数, 利用$ EQ_{1}^{rot} $的插值定义以及$ \vec{W}_h $空间性质,

其中$ l_1, l_3 $分别为$ K $的下边和上边, $ l_2, l_4 $分别为$ K $右边和左边, 则有

$ \begin{align} \ (\nabla_h(u-I_hu), \vec{w}_h)_h = &\sum\limits_{K} \int_{\partial K}(u-I_{K}u)\cdot\vec{n}\cdot\vec{w}_hds-\sum\limits_{K} \int_{K}(u-I_{K}u)\cdot(\nabla\cdot\vec{w}_h)dxdy\\ \ = &\sum\limits_{K}(w_h^2\cdot \int_{l_3}(u-I_{K}u)dx)-\sum\limits_{K} (w_h^2\cdot \int_{l_1}(u-I_{K}u)dx)\\ \ & +\sum\limits_{K} (w_h^1\cdot \int_{l_2}(u-I_{K}u)dy)-\sum\limits_{K} (w_h^1\cdot \int_{l_4}(u-I_{K}u)dy)\\ \ & -\sum\limits_{K}((\nabla\cdot\vec{w}_h)\cdot \int_{K}(u-I_{K_{h}}u)dxdy) = 0. \end{align} $

(4.6)

注意到, 由于该引理的证明利用了$ EQ_1^{rot} $的插值定义, 若此时将单元对换成协调单元对$ Q_{11}\times Q_{10}\times Q_{01} $, 此结果将不再成立.令

$ \begin{eqnarray*} &&u^i-U^i_h = u^i-I_hu^i+I_hu^i-U^i_h\triangleq\eta^i+\xi^i, \\ &&\vec{q}^i-\vec{Q}^i_h = \vec{q}^i-\Pi_h\vec{q}^i+\Pi_h\vec{q}^i-\vec{Q}_h^i\triangleq\vec{r}^i+\vec{\theta}^i, i = 0, 1, 2, \cdots, n. \end{eqnarray*} $

由引理3, 并利用数学归纳法, 可分以下几步分析说明非线性Sobolev方程的无网格比超收敛结果

$ \begin{equation} \|\xi^{n}\|_{h}+\|\vec{\theta}^{n}\|_{H(\operatorname{div};\Omega)} = O(h^{2}+\tau^{2}). \end{equation} $

(4.7)

第一步    $ \|\bar{\partial}_t\xi^n\|_h^2\leq Ch^4+C\tau^4+C\|\xi^n\|_h^2+C\|\vec{\theta}^n\|_0^2. $

第二步   利用$ \|\xi^n\|_h^2\leq C\tau\sum\limits_{i = 1}^n\|\bar{\partial}_t\xi^i\|_h^2 $, 得到

$ \|\bar{\partial}_t\xi^n\|_h^2\leq Ch^4+C\tau^4+C\tau\sum\limits_{i = 1}^n\|\bar{\partial}_t\xi^i\|_h^2+C\|\vec{\theta}^n\|_0^2. $

利用Gronwall引理, 当$ \tau $充分小时有

$ \begin{equation} \|\bar{\partial}_t\xi^n\|_h^2\leq Ch^4+C\tau^4+C\|\vec{\theta}^n\|_0^2, \end{equation} $

(E4.8)

所以

$ \begin{equation} \|\xi^n\|_h^2\leq C\tau\sum\limits_{i = 1}^n\|\bar{\partial}_t\xi^i\|_h^2\leq Ch^4+C\tau^4+C\tau\sum\limits_{i = 1}^n\|\vec{\theta}^i\|_0^2. \end{equation} $

(4.9)

第三步   进一步地,

$ \begin{equation} \|\vec{\theta}^n\|_0^2+\|\nabla\cdot\vec{\theta}^n\|_0^2\leq Ch^4+C\tau^4+C\|\xi^n\|_h^2. \end{equation} $

(4.10)

将(4.9)式代入(4.10)式, 当$ \tau $充分小, 利用Gronwalls引理有$ \|\vec{\theta}^n\|_0+\|\nabla\cdot\vec{\theta}^n\|_0\leq Ch^2+C\tau^2, $再利用(4.9)式得到$ \|\xi^n\|_h\leq Ch^2+C\tau^2. $

这里强调以下三点

(1) 如果直接估计$ \|\xi^n\|_h $, 对$ \tau $和$ h $的比例限制将不可避免;

(2) $ \|\xi^n\|_h^2\leq C\tau\sum\limits_{i = 1}^n\|\bar{\partial}_t\xi^i\|_h^2 $在分析中起到重要的作用;

(3) 在通常的估计中, $ \|\vec{\theta}^n\|_{H({\rm div;\Omega})} $比$ \|\xi^n\|_h $的阶低一阶, 而在这里保持了一样的阶.

在物理上, 双曲方程是一类一直很受关注的偏微分方程, 它可以用来描述声波和电磁波的传播等现象, 其中也有很多文献关注其非线性问题的有限元方法.例如, 文献[48]和[49]分析了非线性双曲方程的全离散格式, 其中文献[48]讨论了混合有限元方法, 达到了最优误差估计.文献[49]利用Galerkin交替方向法讨论了一类三维非线性双曲方程, 利用先验估计的结果得到了误差的$ H^1(\Omega) $ -模和$ L^2(\Omega) $ -模.但上述结果也都没有摆脱$ h $和$ \tau $的比值限制, 在文献[48]和[49]中分别需要假设条件$ \tau = O(h), h^{r} = O(\tau)(1\leq r\leq k+1, k\geq0) $和$ \tau = O(h^2) $.因此, 如何有效的对非线性双曲方程展开无网格比的研究具有相当重要的科学价值.另一方面, 在现有的参考文献中对非线性双曲方程的二阶线性化格式讨论的非常少, 怎样构造新的非线性双曲方程的线性化全离散格式, 使得其有更好的稳定性以及超收敛结果也值得进行深入的探讨.

考虑如下非线性双曲方程

$ \begin{equation} \begin{cases} u_{tt}-\nabla\cdot(a(u)\nabla u) = f(u), & (X, t)\in\Omega\times(0, T], \\ u = 0, & (X, t)\in\ \partial \Omega\times(0, T], \\ u(X, 0) = u_{0}(X), \; \; \; u_t(X, 0) = u_{1}(X), & X\in\Omega. \end{cases} \end{equation} $

(5.1)

关于方程的一些基本假设如同第二节.

我们将对(5.1)式创造性地构造一个新的线性化二阶格式, 技巧性地证明其截断误差的二阶性质, 给出其相对应的时间离散方程解的正则性, 并由此得到在非协调单元下无网格比的超逼近结果.

本节仍采用第二节中非协调$ EQ_1^{rot} $元的空间.记

$ \begin{align} \ & \bar{\partial}_{tt}\sigma^{n} = \frac1{\tau^2}(\sigma^{n+1}-2\sigma^{n}+\sigma^{n-1}), \bar{\partial}_{t}\sigma^{n} = \frac1{2\tau}(\sigma^{n+1}-\sigma^{n-1}) , \tilde{\partial}_{t}\sigma^{n} = \frac1{\tau}(\sigma^{n}-\sigma^{n-1}) , \\ \ &\sigma^{n, \frac14} = \frac14(\sigma^{n+1}+2\sigma^{n}+\sigma^{n-1}), \hat{\sigma}^{n} = \frac14(3\sigma^{n}+2\sigma^{n-1}-\sigma^{n-2}), \bar{\sigma}^{n} = \frac1{2}(\sigma^{n}+\sigma^{n-1}), \end{align} $

则有

$ \begin{align} \ & \bar{\partial}_{t}\sigma^{n} = \tilde{\partial}_{t}\bar{\sigma}^{n+1} , \sigma^{n, \frac14} = \frac12(\bar{\sigma}^{n+1}+\bar{\sigma}^{n}), \hat{\sigma}^{n} = \frac12(3\bar{\sigma}^{n}-\bar{\sigma}^{n-1}), \\ \ & \tilde{\partial}_{t}\hat{\sigma}^{n} = \frac32\bar{\partial}_{t}\sigma^{n-1} -\frac12\bar{\partial}_{t}\sigma^{n-2}, (a^n, \bar{\partial}_{t}b^n) = \tilde{\partial}_{t}(a^n, \bar{b}^{n+1}) -(\tilde{\partial}_{t}a^n, \bar{b}^{n}). \end{align} $

(5.2)

利用这些记号, 考虑(5.1)式线性化的逼近方程:寻找$ U_h^n\in V_{h} $, 使得当$ n\geq 2 $, 有

$ \begin{equation} (\bar{\partial}_{tt}U_h^n, v_h)+(a(\hat{U}_h^n)\nabla_h U_h^{n, \frac{1}{4}}, \nabla_h v_h)_h = (f(\hat{U}_h^n), v_h), \forall v_h\in V_h. \end{equation} $

(5.3)

利用如下方程求解$ U_h^2 $:

$ \begin{equation} (\frac{U_h^{2, 0}-2U_h^1+U_h^0}{\tau^2}, v_h)+(a(\bar{U}_h^1)\nabla_h \frac{U_h^{2, 0}+2U_h^{1}+U_h^{0}}4, \nabla_h v_h)_h = (f(\bar{U}_h^1), v_h), \forall v_h\in V_h \end{equation} $

(5.4)

和

$ \begin{equation} (\bar{\partial}_{tt}U_h^1, v_h)+(a(\frac{U_h^{2, 0}+U_h^{0}}2)\nabla_h U_h^{1, \frac{1}{4}}, \nabla_h v_h)_h = (f(\frac{U_h^{2, 0}+U_h^{0}}2), v_h), \forall v_h\in V_h, \end{equation} $

(5.5)

其中$ U_h^0 = I_hu^0, U_h^1 = I_h(u_0+u_1\tau+\frac{1}{2}u_{tt}(0)\tau^2) $, 且$ u_{tt}(0) $可以利用$ u_{tt}(0) = \nabla\cdot(a(u_0)\nabla u^0)+f(u_0) $得到, $ u_0 $是已知函数.

时间误差   引入时间离散方程, 利用其解$ U^n $分裂误差, 通过估计时间误差得到$ U^n $的正则性.令$ e^n\triangleq u^n-U^n $, 则有

$ \begin{equation} \|e^{2, 0}\|_1+\tau\|e^{2, 0}\|_2+\|\tilde{\partial}_te^n\|_1+\|\bar{e}^{n}\|_2 +\tau\|e^{n}\|_2+\tau\|\tilde{\partial}_t\bar{e}^{n}\|_2\leq C_0\tau^2 \end{equation} $

(5.6)

和

$ \begin{equation} \|\frac{U^{2, 0}-U^{1}}{\tau}\|_2+\|\tilde{\partial}_t U^n\|_2\leq C_0. \end{equation} $

(5.7)

注意到, 在这一节里针对非线性双曲方程构造了一个新的线性化的二阶格式, 可以看到要证明其截断误差为$ O(\tau^2) $是非常不容易的.另一方面, 在得到误差结果时, 由于$ C_0 $在第$ n $层和第$ n+1 $层必须统一, 则在估计误差时, 我们需要利用$ \bar{e}^{m+1}\leq \tau\; (m\leq n-1) $, 而不是$ \bar{e}^{m+1}\leq C_0\tau^2\; (m\leq n-1) $来得到结果.

空间误差   利用以上结果得到超逼近结果, 令$ \theta^n\triangleq I_hU^n-U_h^n $, 则有

$ \begin{align} \|\bar{\theta}^n\|_h\leq C_1 h(h+\tau). \end{align} $

(5.8)

在此过程中, 有以下几点需要特别关注

1.如果将本节使用的$ EQ_1^{rot} $单元换成协调单元(比如$ Q_{11} $单元), 由于没有$ \bar{\partial}_{tt}U^n $的和$ \bar{\partial}_{tt}e^n $的$ H^2 $ -模的有界性, 则结果(5.8)式将不成立.

2. (5.8)式的右端不能直接被估计成$ C_1(h^2+\tau^2) $, 否则$ h $和$ \tau $的比值将无法避免.事实上, 也不能将其估计为$ C_1h $, 因为利用数学归纳法, 我们需要统一$ n $层和$ n+1 $层的, 所以利用了$ \|\bar{\theta}^m\|_h\leq h $.

进一步地, 在误差估计过程中会出现以下项

$ (a(\hat{U}^n)-a(\hat{U}_h^n))\nabla u^{n, \frac{1}{4}}, \nabla_h\bar{\partial}_t\theta^n)_h. $

如果利用类似前面的估计方法, 将$ \tau $从内积的一端转向另一端, 则有结果$ \|\bar{\theta}^n\|_h\leq C_1 (h^2+\tau^2) $, 这时将不可避免的出现网格比.

为了克服关键问题, 重新分裂内积为

$ a(\hat{U}^n)-a(\hat{U}_h^n) = a(\hat{U}^n)-a(I_h\hat{U}^n)+a(I_h\hat{U}^n)-a(\hat{U}_h^n), $

则有

$ \begin{eqnarray*} &&a(I_h\hat{U}^n)-a(\hat{U}^n) = -a_u(\hat{U}^{n})\hat{r}^{n-1}+\frac{1}{2}a_{uu}(\mu_1^n)(\hat{r}^n)^2, \\ &&a(I_h\hat{U}^n)-a(\hat{U}_h^n) = a_u(\hat{U}_h^{n})\hat{\theta}^{n}+\frac{1}{2}a_{uu}(\mu_{2}^n)(\hat{\theta}^n)^2, \end{eqnarray*} $

其中$ \mu_1^n = \hat{U}^{n}+\lambda^n_{1}(I_h\hat{U}^n-\hat{U}^{n}) $, $ \mu_{2}^n = \hat{U}_h^{n}+\lambda^n_{2}(I_h\hat{U}_h^n-\hat{U}_h^{n}) $以及$ r^n = U^n-I_hU^n $.这样可得到

$ \begin{align} \|\tilde{\partial}_{t}\theta^{n+1}\|_0+\|\bar{\theta}^{n+1}\|_h\leq C_{13}h(h+\tau). \end{align} $

最后可导出以下结果

$ \begin{eqnarray*} \|I_h\bar{u}^n-\bar{U}_h^n\|_h&\leq& C\|\bar{e}^n-I_h\bar{e}^n\|_h +C\|\bar{e}^n\|_1+C\|\bar{\theta}^n\|_h\\ &\leq& Ch\|\bar{e}^n\|_2+C\tau^2+Ch(h+\tau)\leq Ch^2+C\tau^2. \end{eqnarray*} $

本文在矩形区域下, 对非线性抛物方程、非线性Schrödinger方程、非线性Sobolev方程、非线性双曲方程, 讨论了它们线性化的全离散格式的构造以及无网格比超收敛分析.

首先, 把原来文献中已有的对非线性发展方程无网格比收敛的结果, 进一步延伸到对非线性发展方程无网格比的超逼近和超收敛的研究中.由于要想得到高精度的结果, 对原始方程解正则性的要求往往都会比较高.那么怎样绕过在矩形区域下, 引入的时间离散方程解的有界性达不到$ H^3(\Omega) $时, 有技巧的得到无网格比超逼近结果是之前文献所不曾见过的.

其次, 对于非线性发展方程, 其非线性项中$ a(u) $的处理也是之前很少有的.我们创新性的利用Taylor展开式, 保持了最后时间方向上的误差不丢失阶.

再次, 在分析过程中发现, 对于一些方程, 若仅仅利用插值算子将会提高对方程解的正则性要求; 仅仅利用投影算子, 没有办法进行插值后处理, 进而没有办法得到整体超收敛的性质.所以在对非线性Schrödinger方程的分析当中, 利用了插值算子和投影算子相结合的处理方式得到了其无网格比超逼近和超收敛的结果.

而后, 由于分裂误差的做法使得引入时间离散方程解在矩形区域下正则性结果不是很好, 而此时想要利用$ H^1 $-Galerkin有限元方法得到无网格比超逼近结果还是非常困难的.对于非线Sobolev方程来说, 我们发现了此类方程自身的特点, 创新性地抛弃了分裂误差的方法, 利用一些特殊的技巧也得到了到了其无网格比超收敛结果.

最后, 由于还没有文献对非线性双曲方程的无网格比讨论过, 特别是对其二阶线性化的全离散格式研究还很少有文献进行报道.我们尝试着构造一个新的线性化二阶格式, 证明了其截断误差的二阶性质, 并给出了无网格比超逼近性质.

对未来的进一步工作, 我们提出以下几点供参考

(1) 各向异性有限元的研究是目前该领域独具特色和挑战性的又一前沿热点和难点(参见文献[50-53]).怎样把原来已有的拟一致的正则剖分下的无网格比超收敛结果拓展到各向异性剖分上去.另一方面, 如何将目前只对矩形区域进行的无网格比超收敛研究拓展到一般区域.

(2) 二重网格算法是处理非线性问题的有效方法(参见文献[54, 55]), 但在误差估计中粗网格的网格参数或者细网格的网格参数都会不可避免的和时间步长产生一个比值.事实上, 我们在文献[22]中已经讨论了半线性抛物方程二重网格算法的无网格比收敛性.那么怎样才可以得到非线性发展方程二重网格算法的无网格比高精度分析呢?由于二重网格算法在粗网格上往往是给出一个非线性的逼近方程, 要去掉网格比约束条件, 便要每一层都引入一个非线性的时间离散方程.这样, 非线性系数的有界性在当时时间层是未知的, 得到它的有界性便是非常困难的事了.另一方面, 二重网格算法需要利用粗网格的结果得到在细网格上的逼近方程, 使得它每一层都是一个线性问题, 那么此时使用光滑边界来提高此处引入的另一组时间离散方程系统解的正则性也许是可行的办法, 但这样再处理超收敛分析就是一个挑战了.

(3) 很多数学物理问题都需要利用积分微分方程进行求解.数学、工程技术和自然科学的许多领域, 例如在系统识别、流体力学、信号重构与图像恢复、电磁场理论与静电学、地球物理勘探等方面都将归结为求解积分微分方程的问题.积分微分方程已经成为众多学者研究的一份重要方向.更值得一提的是, 对于我们比较难以处理的双曲方程, 也可以利用$ u_t = q $使其变成一个抛物型的积分微分方程.对于积分微分方程的无网格比收敛性研究, 现在还没有涉及到, 尤其是该如何解决其积分项的估计都还是待解决的问题, 因此这方面的工作还有待进一步地研究.

(4) 怎样将无网格比超收敛高精度结果进一步推广应用到更为复杂的二阶问题上去, 比如说非定常不可压缩Navier-Stokes方程?事实上, 人们对于该方程的有限元方法做了大量研究, 如协调有限元方法、全离散加罚有限元方法、多尺度方法、稳定化有限元方法等^[56-58].关于无网格比, 文献[59]利用特征有限元方法绕过了Navier-Stokes方程的非线性项的繁琐, 再利用分裂技巧得到了无网格比收敛的结果.那么如果换成一般的有限元方法, 处理它的非线性项使得无网格比结果不掉阶是一个麻烦的事情.关于非线性四阶问题^[59-61] (如EFK方程、Cahn-Hilliard方程、四阶非线性双曲方程等)的无网格比超收敛研究也是我们所关心的话题.事实上, 对于非线性的Navier-Stokes或者是四阶问题, 混合有限元方法是常用的一个方法, 那么怎样利用它得到这些问题的无网格比超收敛结果也是我们未来研究的方向之一.

以上这些问题都是全新的、悬而未决的关键性问题, 作为延伸和后续的工作, 亟待建立一套新的框架性理论体系和分析模式去深入展开对这些问题的研究和探索.可以说, 任何实质性的创新和突破都将进一步丰富有限元方法的内涵, 提升其品味.