考虑扩散方程边值问题

$ \begin{align} -\mathrm{div}( {{\Lambda}{{\nabla}{u}}}) = f, \qquad & \hbox{在 }\Omega\hbox{ 内}, \end{align} $

(1.1)

$ \begin{align} u = g_D, \qquad & \hbox{在 }\Gamma_{D}\hbox{ 上}, \end{align} $

(1.2)

$ \begin{align} - {{\Lambda}{{\nabla}{u}}} \cdot \mathit{\boldsymbol{n}} = g_N, \qquad & \hbox{在 }\Gamma_{N}\hbox{ 上}, \end{align} $

(1.3)

其中$ \Omega \subset \mathbb{R}^{2} $是平面上的多边形区域, ${\Lambda}{}$是扩散张量, $ f $是源项, $ \partial\Omega = \bar{\Gamma}_{D} \cup \bar{\Gamma}_{N} $是$ \Omega $的边界, $ \mathit{\boldsymbol{n}} $为单位外法向量, $ g_{D} $和$ g_{N} $分别是Dirichlet和Neumann边值.

上述扩散问题有着广泛的实际应用背景, 例如在辐射流体力学、辐射磁流体力学、油藏模拟等应用中, 扩散过程与其他物理过程相耦合.我们常常需要在网格强扭曲、强间断、强各项异性、强非线性等极端条件下求解扩散方程, 这是一件具有挑战性的工作.有限体积方法是求解扩散问题最常用的方法之一, 它具有局部守恒、简单容易实现等良好性质.多年来, 许多科研工作者致力于扩散方程的有限体积方法研究, 提出了众多的离散格式, 如九点格式(NPS^[1]、多点通量逼近格式(MPFA)^[2]、支撑算子格式^[3]、广义差分格式^[4]等.按照未知量的类型, 这些格式可大致分为单元中心格式、节点格式、杂交格式、混合格式等等, 详细的最新研究进展可参见文献^[5-9].

我国学者李德元教授等人在二维网格上基于积分插值法构造了一个单元中心型有限体积格式^[1], 由于该格式在结构四边形网格上有九点模板, 故常称其为九点格式, 它是若干辐射流体力学程序的基本格式^{[10, 11]}.由于单元中心格式在一个单元上只有一个未知量, 在构造离散格式时需要引入辅助未知量来提高精度.九点格式的辅助未知量定义在网格节点处, 如何用单元中心未知量逼近节点辅助未知量是九点格式研究中的一个重要内容.一个理想的二维节点插值算法应当具有简单、保正、不依赖于网格拓扑、不依赖于间断线的位置、有二阶精度、易于向三维推广等性质.目前已知的节点插值算法, 如算术平均插值^[1]、逆距离加权插值^[12]、最小二乘插值^[13]、显式加权算法^[14]等, 均只能满足部分要求.由于缺乏容易向三维推广的二阶插值算法, 严重制约了九点格式在三维问题中的应用.

本文在多点通量逼近的边辅助未知量插值算法的基础上, 通过一个特殊的极限技巧获得了一个新的节点插值算法, 并在给定假设下证明了该算法中局部线性系统的可解性.本文插值算法的一个重要之处在于容易向三维推广, 且满足除保正性以外的其他要求.全文的推导基于所谓的线性精确准则^[15], 即当扩散张量关于网格是分片常数、解析解关于网格是分片线性函数时, 算法推导的每一步都是精确成立的.为了和通常的等式相区别, 用符号$ \simeq $表示相关等式是在线性精确意义下成立的.

最初的九点格式是在扩散系数为标量的情形下推导的, 但我们容易将其推广到扩散系数为张量的情形^[14], 为保持本文内容的完整性, 我们简要介绍一下扩散张量情形下九点格式的一种构造算法.首先将$ \Omega $剖分为互不重叠的多边形网格, K表示其中的一个单元, σ表示单元边,K的所有边组成的集合记为ε_K, 单元K关于边σ的单位外法向量记为n_K,σ.记流向量为$ {F} = -\Lambda \nabla u $, 这里假定${\Lambda}{}$关于网格是分片常数, 并用$ {\Lambda}{}_K $表示$ \Lambda $在单元$ K $上的限制.单元K通过边σ的流记为F_K,σ, 即

$ \begin{equation} {{{F}}_{{ {K}}, { {\sigma}}}} = \int_{ {\sigma}}{F} \cdot {{\mathit{\boldsymbol{n}}}_{{ {K}}, { {\sigma}}}}\, ds = -\int_{ {\sigma}}\nabla u \cdot (\Lambda^T_K {{\mathit{\boldsymbol{n}}}_{{ {K}}, { {\sigma}}}})\, ds. \end{equation} $

(2.1)

在单元K上对方程(1)两端积分并利用散度定理, 有

$ \begin{equation} \sum\limits_{ {\sigma}\in {\mathcal{E}_{{ {K}}}}} {{{F}}_{{ {K}}, { {\sigma}}}} = \int_{ {K}}fdxdy. \end{equation} $

(2.2)

于是有限体积法解扩散问题就归结为构造流$ {{{F}}_{{ {K}}, { {\sigma}}}} $的离散表达式.

在图 1中, $ K $和$ L $表示两个相邻单元, 它们的单元中心分别是$ O_K $和$ O_L $, $ \sigma $是它们的公共边, 其顶点为$ A $和$ B $, $ \mathit{\boldsymbol{n}}_{K, \sigma}, \mathit{\boldsymbol{n}}_{L, \sigma} $分别是$ K, L $关于$ \sigma $的单位外法向量. $ d_{K, \sigma}, d_{L, \sigma} $分别表示$ O_K, O_L $到$ \sigma $的距离, $ u_{A}, u_{B}, u_{K}, u_{L} $表示$ u $在点$ A, B, O_K, O_L $处的近似值.对于向量$ v = (a, b)^T $, 记$ v^{\circ} = (b, -a)^T $. 此外, 为叙述简单起见, 引入以下记号(见图 1)

$ \begin{equation} \begin{array}{lll} \mathit{\boldsymbol{t}} = \overrightarrow{BA}, \quad & \mathit{\boldsymbol{a}}_K = \overrightarrow{O_KA}, \quad & \mathit{\boldsymbol{a}}_L = \overrightarrow{O_LA}, \\[0.3cm] \mathit{\boldsymbol{s}} = \overrightarrow{O_KO_L}, \quad & \mathit{\boldsymbol{b}}_K = \overrightarrow{O_KB}, \quad & \mathit{\boldsymbol{b}}_L = \overrightarrow{O_LB}. \end{array} \end{equation} $

(2.3)

首先, 有如下的向量分解

$ \begin{equation} {\Lambda}_{\alpha}^T\mathit{\boldsymbol{t}}^{\circ} = \lambda_{\alpha}^{t}\mathit{\boldsymbol{t}}+\lambda_{\alpha}^{n}\mathit{\boldsymbol{t}}^{\circ}, \quad \alpha = K, L, \end{equation} $

(2.4)

其中

$ \begin{equation*} \lambda_{\alpha}^{n} = \frac{(\mathit{\boldsymbol{t}}^{\circ})^T\Lambda_{\alpha}\mathit{\boldsymbol{t}}^{\circ}}{|\sigma|^2}, \quad \lambda_{\alpha}^{t} = \frac{(\mathit{\boldsymbol{t}}^{\circ})^T\Lambda_{\alpha}\mathit{\boldsymbol{t}}}{|\sigma|^2}, \end{equation*} $

$ |\sigma| $表示$ \sigma $的长度.其次, 在单元$ K $上, 根据文献[14]的引理2.1,

$ \begin{equation} \nabla u\simeq \frac{1}{|\sigma|d_{K, \sigma}}[(u_B-u_K)\mathit{\boldsymbol{a}}_K^{\circ}-(u_A-u_K)\mathit{\boldsymbol{b}}_K^{\circ}]. \end{equation} $

(2.5)

注意到$ |\sigma|\mathit{\boldsymbol{n}}_{K, \sigma} = \mathit{\boldsymbol{t}}^{\circ}, \ t = \mathit{\boldsymbol{a}}_K-\mathit{\boldsymbol{b}}_K, \ {a}_K^{\circ}\cdot {t} = -|\sigma|d_{K, \sigma} $, 将(2.4)和(2.5)式代入(2.1)式, 经过简单的计算后就可以得到

$ \begin{equation} {F}_{K, \sigma}\simeq -\frac{\lambda_{K}^{n}}{|\sigma|d_{K, \sigma}}[(\mathit{\boldsymbol{a}}_K \cdot \mathit{\boldsymbol{t}})(u_{B}-u_K)+(-\mathit{\boldsymbol{b}}_K\cdot\mathit{\boldsymbol{t}})(u_{A}-u_{K})] -(u_{A}-u_{B})\lambda_{K}^{t}. \end{equation} $

(2.6)

类似地, 在单元$ L $上, 有

$ \begin{align} {F}_{L, \sigma}\simeq -\frac{\lambda_{L}^{n}}{|\sigma|d_{L, \sigma}}[(\mathit{\boldsymbol{a}}_L \cdot\mathit{\boldsymbol{t}})(u_{B}-u_L)+(-\mathit{\boldsymbol{b}}_L \cdot\mathit{\boldsymbol{t}})(u_{A}-u_{L})] -(u_{B}-u_{A})\lambda_{L}^{t}. \end{align} $

(2.7)

通常要求通过边$ \sigma $的流满足连续条件, 即

$ \begin{equation} {F}_{K, \sigma} + {F}_{L, \sigma} = 0. \end{equation} $

(2.8)

由此得到

$ \begin{equation} \left(\frac{d_{K, \sigma}}{\lambda_{K}^{n}}+ \frac{d_{L, \sigma}}{\lambda_{L}^{n}}\right) {F}_{K, \sigma} = \frac{d_{K, \sigma}}{\lambda_{K}^{n}}{F}_{K, \sigma} -\frac{d_{L, \sigma}}{\lambda_{L}^{n}} {F}_{L, \sigma}. \end{equation} $

(2.9)

将(2.6)和(2.7)式代入(2.9)式的右端, 利用$ a_K-b_K = a_L-b_L = t, a_K-a_L = b_K-b_L = s $, 并经过简单化简就得到内部边$ \sigma $上流的离散表达式

$ \begin{equation} {F}_{K, \sigma} \simeq -\mathcal{K}_{\sigma}|\sigma|[u_L-u_K-\mathcal{D}_{\sigma}(u_{A}-u_{B})], \end{equation} $

(2.10)

其中

$ \begin{equation*} \mathcal{K}_{\sigma} = \frac{\lambda_{K}^{n}\lambda_{L}^{n}}{\lambda_{L}^{n}d_{K, \sigma}+\lambda_{K}^{n}d_{L, \sigma}}, \quad \mathcal{D}_{\sigma} = \frac{\mathit{\boldsymbol{t}}\cdot\mathit{\boldsymbol{s}}}{\, |\sigma|^{2}}-\frac{1}{|\sigma|}\left(\frac{\lambda_{K}^{t}}{\lambda_{K}^{n}}d_{K, \sigma}+\frac{\lambda_{L}^{t}}{\lambda_{L}^{n}}d_{L, \sigma}\right). \end{equation*} $

容易看出, 上述离散流涉及单元中心未知量$ u_K $和$ u_L $以及节点辅助未知量$ u_A $和$ u_B $, 要得到纯单元中心格式, 需要用单元中心量对节点辅助未知量进行插值.

如图 2(a), 假设$ Q_0 $是一内部节点, 它周围的单元$ \Omega_{k}\; (1\leq k\leq n) $按逆时针顺序排列. $ Q_0P_k $和$ Q_0P_{k+1} $是$ \Omega_k $的以$ Q_0 $为顶点的边, $ O_k $为其中心.假设$ T_k $是边$ Q_0P_k $上的一点, 由下式确定

$ \begin{equation} T_k = (1-\tau) Q_0+\tau P_k, \ \ k = 1, \cdots, n, \end{equation} $

(3.1)

其中$ \tau\in (0, 1) $.令$ u_0 $, $ u_k $, $ \bar{u}_k $分别表示$ Q_0 $, $ O_k $, $ T_k $处的未知量, $ \Lambda_k $表示$ \Lambda $在$ \Omega_k $上所取的分片常数值, $ \mathit{\boldsymbol{n}}_{k} $表示$ \Omega_{k} $关于$ Q_0P_k $的单位外法向, 见图 2(b).本文中$ k $将被视为以$ n $为周期的指标, 例如当$ k = n+1 $和$ 0 $时, 有$ P_{n+1} = P_1 $, $ T_{0} = T_n $.

在多点通量逼近格式中, 边未知量$ \bar{u}_k $可以通过单元中心未知量$ {u}_k $表示.为叙述简洁起见, 引入以下记号

$ \begin{equation} \begin{array}{l} {\boldsymbol p}_{k} = \overrightarrow{Q_{0}P_{k}}, \quad {\boldsymbol s}_{k} = \overrightarrow{Q_{0}O_{k}}, \quad {\boldsymbol t}_{k, 1} = \overrightarrow{O_{k}T_{k}}, \quad {\boldsymbol t}_{k, 2} = \overrightarrow{O_{k}T_{k+1}}, \\[0.3cm] \overline{\bf U} = (\bar{u}_1, \cdots, \bar{u}_{n})^T, \quad {\bf U} = (u_1, \cdots, u_{n})^T, \quad {\mathbb T}_k = \left({\boldsymbol e_k}, \boldsymbol e_{k+1}\right), \end{array} \end{equation} $

(3.2)

其中$ e_k $表示$ n $阶单位矩阵的第$ k $个列向量.显然, $ {\mathbb T}_k $是一个$ n\times 2 $的矩阵, 并且

$ \begin{equation} {\mathbb T}_k^T \, \overline{\bf U} = \left(\begin{array}{l} \bar{u}_k\\[0.2cm] \bar{u}_{k+1} \end{array}\right), \quad {\mathbb T}_k^T\, {\bf U} = \left(\begin{array}{l} {u}_k\\[0.2cm] {u}_{k+1} \end{array}\right). \end{equation} $

(3.3)

在单元$ \Omega_k $的边$ Q_0P_k $和$ Q_0P_{k+1} $上, 定义如下离散流

$ \begin{equation*} f_{k, 1} = \int_{Q_0P_k}\left(- {\Lambda}_{k} \nabla u_{h}^{(k)}\right) \cdot \mathit{\boldsymbol{n}}_k\, ds, \quad f_{k, 2} = \int_{Q_0P_{k+1}}\left(- {\Lambda}_{k} \nabla u_{h}^{(k)}\right) \cdot (-\mathit{\boldsymbol{n}}_{k+1})\, ds, \end{equation*} $

其中$ u_{h}^{(k)} $为$ \triangle O_kT_{k+1}T_k $上标准的$ P_1 $有限元插值.直接计算可得

$ \begin{equation*} \nabla u_{h}^{(k)} = -\frac{1}{{t}_{k, 2}\cdot{t}_{k, 1}^{\circ} }\left((\bar u_{k}-u_k){t}_{k, 2}^{\circ}-(\bar u_{k+1}-u_k){t}_{k, 1}^{\circ}\right). \end{equation*} $

注意到对于$ Q_0P_k $ ($ Q_0P_{k+1} $), 有$ |Q_0P_k|\mathit{\boldsymbol{n}}_k = \mathit{\boldsymbol{p}}_k^{\circ} $ ($ |Q_0P_{k+1}|\mathit{\boldsymbol{n}}_{k+1} = \mathit{\boldsymbol{p}}_{k+1}^{\circ} $).于是

$ \begin{equation} \left( \begin{array}{l} f_{k, 1}\\[0.2cm] f_{k, 2} \end{array} \right) = \left( \begin{array}{ll} a_{11}^{(k)}(\tau)&a_{12}^{(k)}(\tau)\\[0.2cm] a_{21}^{(k)}(\tau)&a_{22}^{(k)}(\tau) \end{array} \right)\left( \begin{array}{l} \bar{u}_k-u_k\\[0.2cm] \bar{u}_{k+1}-u_k \end{array} \right), \end{equation} $

(3.4)

其中

$ \begin{equation} a_{ij}^{(k)}(\tau) = \frac{(-1)^{i+j}}{{t}_{k, 2}\cdot{t}_{k, 1}^{\circ} } \left(\Lambda_k{t}_{k, 3-j}^{\circ}\right)\cdot{p}_{k+i-1}^{\circ}, \ \ i, j = 1, 2, \end{equation} $

(3.5)

这里记号$ (\tau) $表示相关的量是$ \tau $的函数, 在不引起歧义的情况下将其省略.利用(3.3)式, 可将(3.4)式改写为

$ \begin{equation} \left( \begin{array}{l} f_{k, 1}\\[0.2cm] f_{k, 2} \end{array} \right) = \left( \begin{array}{ll} a_{11}^{(k)}&a_{12}^{(k)}\\[0.2cm] a_{21}^{(k)}&a_{22}^{(k)} \end{array} \right){\mathbb T}_k^T\;\overline{\bf U}-\left( \begin{array}{ll} a_{11}^{(k)}+a_{12}^{(k)}&0\\[0.2cm] a_{21}^{(k)}+a_{22}^{(k)}&0 \end{array} \right){\mathbb T}_k^T {\bf U}. \end{equation} $

(3.6)

再利用边$ Q_0P_{k+1} $上的流连续条件, 即

$ \begin{equation} f_{k+1, 1}+f_{k, 2} = 0, \quad k = 1, \cdots, n, \end{equation} $

(3.7)

或者其等价形式

$ \begin{equation*} \sum\limits_{k = 1}^n{\mathbb T}_k \left( \begin{array}{l} f_{k, 1}\\[0.2cm] f_{k, 2} \end{array} \right) = (0, 0, \cdots, 0)^T, \end{equation*} $

就得到

$ \begin{equation} {\mathbb M}\overline{\bf U} = {\mathbb N}{\bf U}, \end{equation} $

(3.8)

其中

$ \begin{equation} \begin{array}{l} \mathbb{M} = \sum\limits_{k = 1}^n{\mathbb T}_k \left( \begin{array}{ll} a_{11}^{(k)}&a_{12}^{(k)}\\[0.2cm] a_{21}^{(k)}&a_{22}^{(k)} \end{array} \right){\mathbb T}_k^T, \quad \mathbb{N} = \sum\limits_{k = 1}^n{\mathbb T}_k \left( \begin{array}{ll} a_{11}^{(k)}+a_{12}^{(k)}&0\\[0.2cm] a_{21}^{(k)}+a_{22}^{(k)}&0 \end{array} \right){\mathbb T}_k^T. \end{array} \end{equation} $

(3.9)

通过求解(3.8)式, 就能够用单元中心未知量来表示边未知量, 即$ \overline {\bf U} = {\mathbb M}^{-1}{\mathbb N}{\bf U} $.容易看出, 只有当$ {\mathbb M} $可逆时, 上述边未知量插值算法才是可行的.自多点通量逼近算法提出来的二十多年里, 人们还没有在数值计算中遇到过$ {\mathbb M} $奇异的情形, 但这一点至今仍缺乏严格的理论证明.

从(3.1)式和边未知量的定义, 容易发现$ \lim\limits_{\tau\rightarrow 0} {T}_k = Q_0, $ $ \lim\limits_{\tau\rightarrow 0} {\bar u}_k = u_0, $ $ k = 1, \cdots, n. $所以一个十分自然的想法就是对(3.8)式应用极限过程来构建$ u_0 $的插值算法.为此, 首先需要搞清楚$ a_{ij}^{(k)} $对$ \tau $的依赖关系.根据(3.1)和(3.2)式, 有$ {t}_{k, 3-j} = -{s}_k+\tau {p}_{k+2-j}, $ $ j = 1, 2, $进而

$ {t}_{k, 2}\cdot{t}_{k, 1}^{\circ} = \tau^2 {p}_{k}^{\circ}\cdot{p}_{k+1}+\tau{s}_{k}^{\circ}\cdot\left( {p}_{k}-{p}_{k+1}\right) = \tau^2 {p}_{k}^{\circ}\cdot{p}_{k+1}+2\tau(S_{k, 1}+S_{k, 2}), $

其中$ S_{k, 1} $ ($ S_{k, 2} $)表示$ \triangle Q_0P_kO_k $ ($ \triangle Q_0O_kP_{k+1} $)的面积.于是(3.5)式可以写为

$ \begin{equation} a_{ij}^{(k)} = \frac{1}{\tau}\left(a_{ij}^{(k, 1)}+\tau a_{ij}^{(k, 2)}\right), \end{equation} $

(3.10)

其中

$ \begin{equation} \begin{array}{l} a_{ij}^{(k, 1)} = \dfrac{(-1)^{i+j}\left(-\Lambda_k {s}_k^{\circ}\right)\cdot{p}_{k+i-1}^{\circ}}{\tau {p}_{k}^{\circ}\cdot{p}_{k+1}+2(S_{k, 1}+S_{k, 2}) }, \quad a_{ij}^{(k, 2)} = \dfrac{(-1)^{i+j}\left(\Lambda_k {p}_{k+2-j}^{\circ}\right)\cdot{p}_{k+i-1}^{\circ}}{\tau {p}_{k}^{\circ}\cdot{p}_{k+1}+2(S_{k, 1}+S_{k, 2})}. \end{array} \end{equation} $

(3.11)

再利用(3.9)式, 可得

$ \begin{equation} \mathbb{M} = \frac{1}{\tau}\left({\mathbb M}_1+\tau{\mathbb M}_2\right), \quad \mathbb{N} = \frac{1}{\tau}\left({\mathbb N}_1+\tau{\mathbb N}_2\right), \end{equation} $

(3.12)

其中

$ \begin{equation*} \begin{array}{l} \mathbb{M}_l = \sum\limits_{k = 1}^n{\mathbb T}_k \left( \begin{array}{ll} a_{11}^{(k, l)}&a_{12}^{(k, l)}\\[0.2cm] a_{21}^{(k, l)}&a_{22}^{(k, l)} \end{array} \right){\mathbb T}_k^T, \quad \mathbb{N}_l = \sum\limits_{k = 1}^n{\mathbb T}_k \left( \begin{array}{ll} a_{11}^{(k, l)}+a_{12}^{(k, l)}&0\\[0.2cm] a_{21}^{(k, l)}+a_{22}^{(k, l)}&0 \end{array} \right){\mathbb T}_k^T, \quad l = 1, 2. \end{array} \end{equation*} $

需要注意, 由于(3.11)式的分母中含有$ \tau $, $ {\mathbb M}_l $和$ {\mathbb N}_l $依然是$ \tau $的函数.根据(3.11)式, 发现

$ \begin{equation} a_{i1}^{(k, 1)}+a_{i2}^{(k, 1)} = 0, \ \ \ i = 1, 2. \end{equation} $

(3.13)

故有

$ \begin{equation} {\mathbb M}_1{\bf 1} = {\bf 0}, \ \ {\mathbb N}_1 = {\mathbb O}, \end{equation} $

(3.14)

其中$ {\bf 1} $和$ \bf{0} $分别表示所有分量为1和0的向量, $ {\mathbb O} $表示零矩阵.利用(3.12)和(3.14)式, 可证(3.8)式等价于

$ \begin{equation} \left({\mathbb M}_1+\tau{\mathbb M}_2\right)\overline{\bf U} = \tau {\mathbb N}_2{\bf U}. \end{equation} $

(3.15)

注意到解沿着$ Q_0P_k\; (1\leq k\leq n) $是切向连续可微的, 通过泰勒展开, 可以得到

$ \begin{equation} {\overline {\bf U}} = u_0{\bf 1}+\tau{\bf \Gamma} +{\bf O}(\tau^2), \end{equation} $

(3.16)

其中$ {\bf \Gamma} = (\gamma_1, \cdots, \gamma_{n})^T $是与$ \tau $无关的常向量, $ {\bf O}(\tau^2) $表示所有分量为$ O(\tau^2) $的向量.将(3.16)式代入(3.15)式并且利用(3.14)式中的第一个等式, 得到

$ \begin{equation*} u_0\tau{\mathbb M}_2{\bf 1}+ \tau {\mathbb M}_1\left({\bf \Gamma}-\gamma_1{\bf 1}\right) +\tau^2{\mathbb M}_2{\bf \Gamma}+ \left({\mathbb M}_1+\tau{\mathbb M}_2\right){\bf O}(\tau^2) = \tau {\mathbb N}_2{\bf U}. \end{equation*} $

最后将以上方程两端同时除以$ \tau $, 并令$ \tau\rightarrow 0 $, 就得到最终的方程

$ \begin{equation} u_0{\mathbb M}_2(0){1}+{\mathbb M}_1(0)\left({\bf \Gamma}-\gamma_1{\bf 1}\right) = {\mathbb N}_2(0) {\bf U}, \end{equation} $

(3.17)

其中$ {\mathbb M}_1(0) $, $ {\mathbb M}_2(0) $和$ {\mathbb N}_2(0) $表示$ {\mathbb M}_1 $, $ {\mathbb M}_2 $和$ {\mathbb N}_2 $取$ \tau = 0 $的值.将$ {\mathbb M}_1(0) $的第一列替换为$ {\mathbb M}_2(0){1} $其余列不变得到一个新的矩阵, 记为$ {\widetilde {\mathbb M}} $, 则(30)式可以写为与$ \tau $无关的线性系统

$ \begin{equation} \widetilde{\mathbb M}\left(u_0, \gamma_2-\gamma_{1}, \cdots, \gamma_{n}-\gamma_{1}\right)^T = {\mathbb N}_2(0) {\bf U}, \end{equation} $

(3.18)

求解上述局部线性系统就得到一个新的节点插值算法, 即

$ \begin{equation} u_0 = (1, 0, \cdots, 0)\widetilde{\mathbb M}^{-1}{\mathbb N}_2(0){\bf U}. \end{equation} $

(3.19)

这里需要特别强调, 虽然上述算法基于多点通量逼近的边辅助未知量插值, 但它与(21)式中矩阵$ {\mathbb M} $的可逆性没有关系.

虽然(3.8)式的可解性目前尚无理论证明, 但局部线性系统(3.18)式的可解性在一定条件下是可以严格证明的.本节的主要结果依赖以下两个假设:

(H1) $ {p}_{k+i-1}^{\circ}\cdot(\Lambda_k{p}_{k+i-1}^{\circ}-\Lambda_k{p}_{k+2-i}^{\circ})>0, \ 1\leq k\leq n, \ i = 1, 2 $;

(H2) $ \left(\Lambda_k {s}_k^{\circ}\right)\cdot{p}_{k+i-1}^{\circ}>0, \ 1\leq k\leq n, \ i = 1, 2 $.

在扩散系数是标量的情况下, 假设(H1)和(H2)的几何意义是明显的, 此时, 两个假设变为

(H1)' $ {p}_{k+i-1}\cdot({p}_{k+i-1}-{p}_{k+2-i})>0, \ 1\leq k\leq n, \ i = 1, 2 $;

(H2)' $ {s}_k\cdot{p}_{k+i-1}>0, \ 1\leq k\leq n, \ i = 1, 2 $.

如图 3所示, (H1)'表示$ \theta_1^k $和$ \theta_2^k $是锐角. (H2)'表示$ \beta_1^k $和$ \beta_2^k $是锐角.

引理 4.1  在假设(H1)下, $ {\mathbb M}_2(0){1} $的所有元素都是负的, 即

$ \begin{equation} {e}_k^T{\mathbb M}_2(0){1}< 0, \quad k = 1, 2, \cdots, n. \end{equation} $

(4.1)

证  由$ {\mathbb M}_2(0) $的定义, 有

$ \begin{equation*} \begin{array}{rl} {e}_k^T{\mathbb M}_2(0){1}& = \sum\limits_{k^{\prime} = 1}^n{e}_k^T{\mathbb T}_{k^{\prime}} \left( \begin{array}{ll} a_{11}^{({k^{\prime}}, 2)}(0)&a_{12}^{({k^{\prime}}, 2)}(0)\\[0.2cm] a_{21}^{({k^{\prime}}, 2)}(0)&a_{22}^{({k^{\prime}}, 2)}(0) \end{array} \right){\mathbb T}_{k^{\prime}}^T{1}\\[0.7cm] & = \sum\limits_{{k^{\prime}} = 1}^n\left(\delta_{{k^{\prime}}, k}, \delta_{{k^{\prime}}+1, k}\right)\left( \begin{array}{ll} a_{11}^{({k^{\prime}}, 2)}(0)&a_{12}^{({k^{\prime}}, 2)}(0)\\[0.2cm] a_{21}^{({k^{\prime}}, 2)}(0)&a_{22}^{({k^{\prime}}, 2)}(0) \end{array} \right)\left( \begin{array}{l} 1\\[0.2cm] 1 \end{array} \right)\\[0.7cm] & = a_{11}^{(k, 2)}(0)+a_{12}^{(k, 2)}(0)+a_{21}^{(k-1, 2)}(0)+a_{22}^{(k-1, 2)}(0), \end{array} \end{equation*} $

其中$ \delta_{ij} $表示Kronecker delta函数.利用(24)式并且注意到$ \tau = 0 $, 进一步得到

$ \begin{equation*} {e}_k^T{\mathbb M}_2(0){1} = -\frac{{p}_{k}^{\circ}\cdot\left( \Lambda_k{p}_{k}^{\circ}-\Lambda_k{p}_{k+1}^{\circ}\right)}{2(S_{k, 1}+S_{k, 2})}-\frac{{p}_{k}^{\circ}\cdot\left( \Lambda_{k-1}{p}_{k}^{\circ}-\Lambda_{k-1}{p}_{k-1}^{\circ}\right)}{2(S_{k-1, 1}+S_{k-1, 2})}. \end{equation*} $

再结合假设(H1), 立即得到(4.1)式.

引理 4.2  记$ \mathbb A = (a_{ij}) $为$ n\times n $的矩阵, 满足

$ \begin{equation} {\mathbb A} = \sum\limits_{k = 1}^n{\mathbb T}_k \left( \begin{array}{rr} a_k&-a_k\\[0.2cm] -b_k&b_k \end{array} \right) {\mathbb T}_k^T, \end{equation} $

(4.2)

其中$ a_k, b_k\; (1\leq k\leq n) $是正数.则

$ \begin{equation} a_{ij}^* = a_{i1}^*\ \ \hbox{且}\ \ a_{ii}^*>0, \quad 1\leq i, j\leq n, \end{equation} $

(4.3)

其中$ a_{ij}^* $表示$ a_{ij} $对应的代数余子式.

证  根据(4.2)式, 有

$ \begin{equation*} {\mathbb A}{1} = \sum\limits_{k = 1}^n{\mathbb T}_k \left( \begin{array}{rr} a_k&-a_k\\[0.2cm] -b_k&b_k \end{array} \right) {\mathbb T}_k^T{1} = \sum\limits_{k = 1}^n{\mathbb T}_k \left( \begin{array}{rr} a_k&-a_k\\[0.2cm] -b_k&b_k \end{array} \right)\left( \begin{array}{r} 1\\[0.2cm] 1 \end{array} \right) = {0}. \end{equation*} $

由此并根据$ a_{ij}^* $和$ a_{i1}^* $的定义, (35)式的第一部分很容易验证.由$ \mathbb A $的定义有

$ \begin{equation} a_{nn}^* = \left|\begin{array}{ccccccc}a_1+b_n& -a_1 &0&0&\cdots&0 &0\\ -b_1& a_2+b_1 &-a_2&0&\cdots & 0&0\\ 0& -b_2 &a_3+b_2&-a_3&\cdots & 0&0\\ \cdots& \cdots &\cdots&\cdots&\cdots &\cdots& \cdots\\ 0&0&0&0&\cdots&a_{n-2}+b_{n-3}& -a_{n-2}\\ 0&0&0&0&\cdots&-b_{n-2}&a_{n-1}+b_{n-2} \end{array} \right|: = A_{n-1}. \end{equation} $

(4.4)

分裂$ A_{n-1} $的最后一列, 可得

$ \begin{array}{rl} A_{n-1} = &\left|\begin{array}{ccccccc}a_1+b_n& -a_1 &0&0&\cdots&0 &0\\ -b_1& a_2+b_1 &-a_2&0&\cdots & 0&0\\ 0& -b_2 &a_3+b_2&-a_3&\cdots & 0&0\\ \cdots& \cdots &\cdots&\cdots&\cdots &\cdots& \cdots\\ 0&0&0&0&\cdots&a_{n-2}+b_{n-3}& -a_{n-2}\\ 0&0&0&0&\cdots&-b_{n-2}&b_{n-2} \end{array} \right|\\ & +a_{n-1}A_{n-2}. \end{array} $

对于$ A_{n-1} $的第一部分, 通过将最后一列加到第$ (n-2) $列, 再将新的第$ (n-2) $列加到第$ (n-3) $列, 依此类推, 最终可得到一个上三角型行列式, 进而有

$ \begin{equation} A_{n-1} = b_n\prod\limits_{i = 1}^{n-2} b_i+a_{n-1}A_{n-2}. \end{equation} $

(4.5)

用同样的方法可以得到

$ A_{k} = b_n\prod\limits_{i = 1}^{k-1} b_i+a_{k}A_{k-1}, \quad 2\leq k\leq n-2. $

注意到$ A_1 = a_1+b_n>0 $以及$ a_i, b_i>0 $, 通过数学归纳法最终得到$ A_{n-1}>0 $, 从而$ a_{nn}^*>0 $.类似地, 可以证明$ a_{11}^*>0 $.当$ i\ne 1, n $时, 记$ a_{ii} $对应的余子式为$ \mathbb{A}_{ii}^{*} $, 并令

$ \begin{equation*} \mathbb{V}_{i} = \begin{bmatrix} \mathbb{O} & \mathbb{I}_{i-1} \\ \mathbb{I}_{n-i} & \mathbb{O} \end{bmatrix}, \end{equation*} $

其中$ \mathbb{I}_{k} $是$ k $阶单位矩阵.容易验证$ \mathbb{V}_i^{T} \mathbb{A}_{ii}^{*}\mathbb{V}_i $是形如(36)式的三对角矩阵.于是

$ \begin{equation*} a_{ii}^{*} = (-1)^{i+i}\det(\mathbb{A}_{ii}^{*}) = \frac{1}{\left(\det(\mathbb{V}_i)\right)^2}\det(\mathbb{V}_i^{T} \mathbb{A}_{ii}^{*}\mathbb{V}_i)>0, \quad 1<i<n. \end{equation*} $

定理 4.3  在假设(H1)和(H2)的条件下, (31)式中的系数矩阵$ \widetilde{\mathbb M} $是非奇异的.

证  由$ {\mathbb M}_1(0) $的定义和(3.13)式, 有

$ \begin{equation} {\mathbb M}_1(0) = \sum\limits_{k = 1}^n{\mathbb T}_k \left( \begin{array}{rr} a_{11}^{(k, 1)}(0)&-a_{11}^{(k, 1)}(0)\\[0.2cm] -a_{22}^{(k, 1)}(0)&a_{22}^{(k, 1)}(0) \end{array} \right) {\mathbb T}_k^T. \end{equation} $

(4.6)

由(3.11)式和假设(H2), 可以得到

$ \begin{equation} a_{ii}^{(k, 1)}(0) = -\frac{\left(\Lambda_k {s}_k^{\circ}\right)\cdot{p}_{k+i-1}^{\circ}}{2(S_{k, 1}+S_{k, 2}) }<0, \quad i = 1, 2. \end{equation} $

(4.7)

令$ \mathbb A = (a_{ij})_{n\times n} = -{\mathbb M}_1(0) $.根据引理4.2, $ a_{k1}^*>0\; (1\leq k\leq n) $.由$ \widetilde{\mathbb M} $的定义, 并利用引理4.1, 可得

$ \det(-\widetilde{\mathbb M}) = \sum\limits_{k = 1}^n\left(-{e}_k^T{\mathbb M}_2(0){1}\right)a_{k1}^*>0, $

也就是说$ \widetilde{\mathbb M} $是非奇异的.

首先, 定义如下的$ L_{2} $误差和流误差

$ \begin{align*} E_{u} = \left(\sum\limits_{K\in\mathcal{M}}|K|e_K^{2}\right)^{\frac{1}{2}}, \quad E_{q} = \left(\sum\limits_{K\in\mathcal{M}}\sum\limits_{\sigma\in {\mathcal{E}_{{ {K}}}} }\frac{|\sigma|}{N_{\sigma}(d_{K, \sigma}+d_{L, \sigma})}(e_{K}-e_{L})^{2}\right)^{\frac{1}{2}}, \end{align*} $

其中$ \cal M $为单元集合, $ e_{K} $表示单元$ K $中心的误差, $ N_{\sigma} $表示边$ \sigma $的相邻单元数, $ L $表示$ \sigma $的除$ K $之外的另一相邻单元（如果存在的话）, 当$ \sigma $在边界上时, 我们约定$ d_{L, \sigma} = e_L = 0 $.两个网格层之间的收敛率$ R_{\alpha} $ ($ \alpha = u, q $)通过以下公式获得

$ \begin{equation*} {\frac{\log[E_{\alpha}(h_2)/E_{\alpha}(h_1)]}{\log(h_2/h_1)}, } \end{equation*} $

其中$ h_1, h_2 $代表两个网格层的网格尺寸, $ E_{\alpha}(h_1), E_{\alpha}(h_2) $为对应的离散误差.

用BICGSTAB方法求解离散线性系统, 并取收敛准则为$ {\varepsilon}_{{\rm lin}} $ = 1.0E-15.离散流采用公式(13), 节点插值算法分别采用文献[14]中的显式加权算法2和本文的极限加权算法, 对应的格式简记为NPS-EW2和NPS-LW.采用的计算网格见图 4, 且所有的计算均在双精度下完成.

考虑$ \Omega = [0, 1]^2 $上的具有全Dirichlet边界的扩散问题(1).扩散张量和精确解如下

$ \begin{equation*} {\Lambda} = \begin{pmatrix} 1.5 & 0.5 \\ 0.5 & 1.5 \end{pmatrix}, \quad u(x, y) = \frac{1}{2} \left[ \frac{\sin\!\left(\left(1- x\right)\, \left(1 - y\right)\right)} {\sin(1)} + {\left(1 - x\right)}^3\, {\left(1 - y\right)}^2 \right]. \end{equation*} $

这个数值实验是文献[5]中的一个测试的一个简单变形.

从实验结果来看, NPS-LW在各网格上与NPS-EW2的实验结果接近, 这里仅给出前者的计算结果. 表 1给出了NPS-LW在各个网格上数值解和流的误差, 图 5以对数折线图的形式展示了收敛速度.可以看出NPS-LW在这些网格上的数值解的误差随着网格加密以趋近2阶的速度收敛.

在$ [0, 1]^2 $上求解泊松方程$ -\Delta u = 4.0 $, 采用Dirichlet边界条件, 解析解取为

$ \begin{equation*} u(x, y) = -(x^2+y^2) + 2. \end{equation*} $

本算例十分简单, 其设计目的是测试新的节点加权算法与已有显式加权算法在数值表现上的差异. NPS-LW的数值结果在Mesh1, Mesh2和Mesh4上与NPS-EW2非常接近, 但在Mesh3上明显优于NPS-EW2, 在表 2中进行了对比.

在$ \Omega = [0, 1]^2 $上求解扩散问题(1.1)–(1.2), 扩散张量取为

$ \begin{equation*} {\Lambda} = \begin{cases} \mathbb{I}, &0< x\leqslant 0.5, \\ 10^{-3}\, \mathbb{I}, & 0.5 < x < 1, \end{cases} \end{equation*} $

$ \mathbb{I} $为二阶单位矩阵, 精确解取为

$ \begin{equation*} u(x, y) = \begin{cases} 1+x+y+(x-0.5)^2\mathrm{e}^{x+y}, & x\leqslant0.5, \\ 10^3\, x + y + {\left(x - 0.5\right)}^2\mathrm{e}^{x + y} - 498.5, & x >0.5. \end{cases} \end{equation*} $

这里采用Mesh1和Mesh4.对于Mesh1, 所有位于$ x = 0.5 $上的节点在$ x $方向上都不扰动.数值计算结果见表 3.从该表中可以看到, 数值解误差收敛速度均趋于2阶, 流的误差接近1阶, 均为最优阶, 表明新的节点加权算法也能很好地适应间断系数问题.