奇妙族 (serendipity families)

xsaisigma

奇妙族是一类 Lagrange 型有限元, 它在力学, 建筑学和其它工学的模拟软件里是用得很多的. 对于计算数学爱好者, 它的构造也非常值得学习: 它充分利用了节点的位置关系, 使得不仅在较少的节点下能获得同样好的效果, 并在这个过程中, 把节点都安排在每个元的边界上, 从而规避了节点在内部时函数的光滑性难以控制的弊端.

我们从最简单的 \mathbb R^2 中的奇妙有限元开始. 在这之前, 先快速提一下有限元的定义 "3-字母":

一个在 \mathbb R^n 中的有限元是指 (K, \Sigma, P), 其中

K 是 \mathbb R^n 中有者非空内部和 Lipschitz 边界的子集;
\Sigma 是一些定义在 C^\infty(K) 上的线性无关的线性型 T_i, i = 1, \cdots, N 的集合, 这些线性型就是有限元的自由度;
P 是函数 K \rightarrow \mathbb R 的空间, 它里面的函数满足如下性质: 给定实数 \alpha_i, i = 1,\cdots, N, 存在唯一的 p \in P 满足 T_i(p) = \alpha_i. 换句话说, P 里面有着一组对偶于 \{T_i\} 的基. 这个基就是有限元的基函数.

1) 现在请看下图中的正方形有限元:

这里, a_1 到 a_4 是正方形的的顶点, a_5 到 a_8 是四边中点, a_9 是正方形的中心. 为了方便, 我们先把图中这个正方形看成参考单元, 即单位正方形 [0,1]^2.

如果您想在正方形元上使用分段二次多项式的基函数, 您通常考虑空间 Q_2, 即变量 x_1 和 x_2 各次数不超过 2 的空间. 在这时, 由于
Q_2 = \text{span} \{1, x_1, x_2, x_1x_2, x_1^2, x_2^2, x_1^2 x_2, x_1x_2^2, x_1^2 x_2^2\},
\dim Q_2 = 9 (通常来讲, \mathbb R^n 上的 Q_k 的维数是 (k+1)^n, 很好直接证明), 您直接使用图中的九个点就是正常的. 但这时, 由于 a_9 在元的内部, Lagrange 型元所要求的 C^0 连续性在这个点不好控制; 因为如果是一个完整网格的话, 每个 a_9 都分别在每个正方形内部, 和其它正方形不容易联系.

可是, 这些节点和顶点之间是有中点的位置关系的. 具体的写出来的话: 让 \pi_1 和 \pi_2 分别表示取出第一和第二个直角坐标, 则对于中心
\pi_i(a_9) = \frac{\pi_i(a_1) + \pi_i(a_2) + \pi_i(a_3) + \pi_i(a_4)}{4}, \quad i = 1, 2;
而对于各边中点
\begin{aligned} \pi_i(a_5) &= \frac{\pi_i(a_1) + \pi_i(a_2)}{2}, \quad i = 1, 2;\\ \pi_i(a_6) &= \frac{\pi_i(a_2) + \pi_i(a_3)}{2}, \quad i = 1, 2;\\ \pi_i(a_7) &= \frac{\pi_i(a_3) + \pi_i(a_4)}{2}, \quad i = 1, 2;\\ \pi_i(a_8) &= \frac{\pi_i(a_4) + \pi_i(a_1)}{2}, \quad i = 1, 2.\\ \end{aligned}

也就是说, 中心坐标和顶点坐标之间是线性关系, "系数" 是 4, 边中点坐标和顶点坐标之间也是线性关系, "系数" 是 2. 因此, 应该有一个形如
4 \times \text{ 中心 } + 1 \times \text{ 顶点 } - 2 \times \text{ 边中点 } = 0
的恒等式.

实际上确实是这样. 这个恒等式不仅在 \pi_i 作用下成立, 而且在任意多项式 p \in P_2 作用下都成立. 我们把它表述成

命题: 对每个 p \in P_2, 有恒等式 \psi(p) := 4p(a_9) + \sum_{i=1}^4 p(a_i) - 2\sum_{i=5}^8 p(a_i) = 0.

为了证明它, 只需要把 \psi 当成是 P_2^* (对偶空间) 中的元素, 而
P_2 = \text{span }\{1, x_1, x_2, x_1 x_2, x_1^2, x_2^2\},
验证 \psi 在每个基上为零就行了.

在这里, 您最终代入点 a_i = (a_{1i}, a_{2i}) 的具体坐标. 可是, 即使不把它看成单位正方形, 对于任意正方形, 证明还是成立的. 因为这个性质本身不是依赖于单位正方形, 而是正方形中中心和边中点和顶点的位置关系.

用这个恒等式, 就可以 "消去" a_9. 但是还是有一个矛盾: 恒等式对 P_2 型的多项式成立, 而 "消去" a_9 之后还有八个节点, 显多; 何况您本来是打算用 Q_2 型的多项式的.

所以, 转而使用这个多项式空间 \tilde Q_2 := \{ p \in Q_2: \psi(p) = 0\}.

我们来分析 \tilde Q_2 的维数. 具体的做法是用 Q_2 = P_2 + \{x_1^2 x_2\} + \{x_1 x_2^2\} + \{x_1^2 x_2^2\} (这里的 + 表示直和), 然后看 \psi 在多出来的那三个直因子上的作用情况. 直接演算可以发现 (您可以检验这一点): x_1^2 x_2 和 x_1 x_2^2 可以零化 \psi, 而 x_1^2 x_2^2 不能. 所以, \tilde Q_2 的维数正好是 8, 而 "消去" a_9 后正好有 8 个在边上的节点. 如果 \{p(a_i), i = 1,\cdots,8\} 可以作为自由度, 情况就完美了.

注意人们通常说的 p(a_i) 作为自由度的意思 --- 按照定义, 自由度应该是一个线性泛函. 所以这句话在这里的准确涵义是: 定义
T_i(p) = p(a_i), \quad i = 1,\cdots,8,
这些 T_i 是自由度.

为了验证它们是自由度, 需要证明这些 T_i 线性无关. 利用下列事实:

如果对每个 i 有 p(a_i) = p(a_{1i}, a_{2i}) = 0, 则 p = 0.

这一点很容易直接验证. 而 T_i(p) = p(a_i), 如果设 T = [T_1, \cdots, T_8]^T 为 \tilde Q_2 \rightarrow \mathbb R^8, 说明 T(p) = 0 \Rightarrow p = 0, 即 \ker T = \{0\}. T 进而是单射; 而一个有限维的单射就是双射. 所以对每个 y = [y_1, \cdots, y_8]^T \in \mathbb R^8, 存在一个 p 使得 T(p) = y.
现在作出线性组合
c_1L_1 + c_2L_2 + \cdots + c_8 L_8 = 0.
它等价于
c_1 y_1 + c_2 y_2 + \cdots + c_8 y_8 = 0\quad \text{ 对所有 $y_1,\cdots, y_8 \in \mathbb R$ 成立}.
这时候, 只要取 y_i 为 \mathbb R^8 中的第 i 个标准基向量, 就得到所有的系数 c_i = 0, 线性无关性得证. 从这里, 您可以根据对偶基的关系自然求出这种有限元的基函数.

另外, 此时的有限元空间当然是 C^0 的子空间. 因为只有公共边上有公共节点, 而在这些节点处, 两边的两个二次多项式取值相等; 也就是说, 两个二次多项式在三个点处相等, 从而它们恒等.

xsaisigma

2) 假设您又想在正方形上改用 Q_3 型的多项式. 您现在应该知道, 应该改取每边的三等分点, 然后 "消去" 正方形内部的四个节点, 如下图所示,

您消去 a_{13} 到 a_{16}, 也就是说, 这次要找到四个基于节点位置关系的恒等式 \psi_1, \cdots, \psi_4.

作为例子, 我们来作出用于消去节点 a_{13} 的恒等式 \psi_1. 相对地, 我们为此先只考虑顶点, 以及 a_{13} 所在的一条横线和一条竖线上的三等分点, 即下图中的五点 "十" 字.

按照定比分点的关系, 设 \pi_1 和 \pi_2 分别表示取出第 1, 2 个直角坐标, 则对 a_{13} 而言
\pi_1(a_{13}) = \frac 23 \pi_1(a_{12}) + \frac 13 \pi_1(a_{6}) \Rightarrow 2\pi_1(a_{12}) - 3\pi_1(a_{13}) + \pi_1(a_{6}) = 0
及
\pi_2(a_{13}) = \frac23 \pi_2(a_{5}) + \frac 13 \pi_2(a_{11}) \Rightarrow 2 \pi_2(a_{5}) - 3\pi_2(a_{13}) + \pi_2(a_{11}) = 0.

这里就产生了和上一个奇妙元不同的地方. 虽然还是有着

3\times \text{三等分内点} - 2 \times \text{较低的位于边上的点} - 1 \times \text{较高的位于边上的点} = 0

的线性恒等式, 但在所考虑的坐标分量不同时, 涉及不同的节点, 另一方面又有着相同的系数; 也就是说, 目前有着沿 x_1 和 x_2 的各一个恒等式, 它们分别涉及 x_1 或 x_2 方向的不同节点, 但在对应地位的点上, 系数又是相同的. 这些恒等式中的每一个很明显至少对 \mathbb R^1 上的 P_1 多项式都是成立的. 我们把它表述成

命题: 设 p \in P_1 为一元多项式, t \in [0,1] 是不定元 x = x(t) \in [a,b] 的参数. 则 \Lambda_x(p) = 3p(x(t = 1/3)) - 2p(x(t = 0)) - p(x(t = 1)) = 0.

很容易直接证明, 因为它对每个坐标分量都已经成立了, 而 P_1(\mathbb R) 无非就是 \text{span }\{1, x\}. 剩下的就是要找到一种办法, 把沿着 x_1 和 x_2 的两个这样的恒等式合在一起, 得到一个 \mathbb R^2 上的恒等式.

和之前的例子中做的一样, 我们先考虑 P_3 型的三次多项式, 然后再过渡到 Q_3. 这个空间是
P_3 = \text{span } \{x_1^{r_1} x_2^{r_2}\}, \quad r_1 + r_2 \leq 3.

也就是说, 在每一个基向量里, 至少有一个不定元的次数是 0 或 1. 那只要把沿着 x_1 和 x_2 的 \Lambda(p) "乘" 起来, 不就行了吗? 当然, 您想让这两个一元多项式的泛函 "乘" 出一个二元多项式的泛函, 所以就使用张量积
\begin{aligned} \Lambda_{x_1}(p) \otimes \Lambda_{x_2}(p) &= (3p(x_1(t = 1/3)) - 2p(x_1(t = 0)) - p(x_1(t = 1)))\otimes (3p(x_2(t = 1/3)) - 2p(x_2(t = 0)) - p(x_2(t = 1))) = \\ &= 9\tilde p(x_1(t = 1/3), x_2(t = 1/3)) - 6\tilde p(x_1(t = 1/3), x_2(t = 0)) - 3\tilde p(x_1(t = 1/3), x_2(t = 1)) - \\ &-6\tilde p(x_1(t = 0), x_2(t = 1/3)) + 4\tilde p(x_1(t =0), x_2(t = 0)) + 2\tilde p(x_1(t = 0), x_2(t = 1) -\\ &-3\tilde p(x_1(t = 1), x_2(t = 1/3)) +2 \tilde p(x_1(t = 1), x_2(t = 0)) + \tilde p(x_1(t = 1), x_2(t = 1)) :=\\ &:= \psi_1(\tilde p), \end{aligned}

这里的 \tilde p, 按理来说是定义在 \text{span }\{1 \otimes 1, 1\otimes x_1, 1 \otimes x_2, x_1 \otimes x_2\} \simeq \{1, x_1, x_2, x_1x_2\} 上的: \tilde p(x_1, x_2) = p(x_1) p(x_2), 现在我们把它延拓到整个 P_3 上面, 从而 \psi_1(\tilde p) = \Lambda_{x_1}(p) \otimes \Lambda_{x_2}(p) 也被延拓. 最后, 造成如下的奇观: P_3 型多项式中的每一项, x_1 和 x_2 至少有一个次数 \leq 1. 所以当某项中 x_1 的次数 \leq 1 时, \Lambda_{x_1} 负责零化这项; 而当某项中 x_2 的次数 \leq 1 时, \Lambda_{x_2} 负责零化这项. 综合起来, \psi_1 可以零化所有的 \tilde p \in P_3.

如果您把所考虑的 "十字" 中的坐标代进去, 最后的恒等式就是
\begin{aligned} \psi_1(\tilde p) &= 9 \tilde p(a_{13}) + 4\tilde p(a_1) + 2\tilde p(a_2) + \tilde p(a_3) + 2\tilde p(a_4) - 6\tilde p(a_5) -\\ &- 3\tilde p(a_6) - 3\tilde p(a_{11}) - 6\tilde p(a_{12}) = 0. \end{aligned}
您同样的考虑其它三个内点对应的 "十字", 得到其它三个恒等式. 从形式上讲, 它们就是在上式中循环置换变量
\{a_1, a_2, a_3, a_4\}, \{a_5, a_6, a_7, a_8\}, \{a_9, a_{10}, a_{11}, a_{12}\}, \{a_{13}, a_{14}, a_{15}, a_{16}\}
而得到 (我说的是在每一组内独立循环置换, 且四组都要置换).

接下来的事情就没什么困难了. 您转向子空间 \tilde Q_3 = \{\tilde p \in Q_3: \psi_i(p) = 0, \,i = 1,2,3,4\}, 然后用
Q_3 = P_3 + \{x_1^2 x_2^2\} + \{x_1^3 x_2\} + \{x_1 x_2^3\} + \{x_1^3 x_2^2\} + \{x_1^2 x_2^3\} + \{x_1^3 x_2^3\}
来联系 P_3 和 Q_3, 还检查每个恒等式 \psi_i 在这些多出来的直因子上的作用. 您发现每个 \psi_i 可以零化其中的四个, 而不能零化另两个, 所以 \dim \tilde Q_3 = 12. 另一方面, 在 "消去" 了四个内部节点之后, 正好剩下 12 个在边界上的节点, 而可以验证 \{\tilde p(a_i),\,i = 1,\cdots, 12\} 线性无关, 正好可以作为自由度. 这样形成的有限元空间也是 C^0 的子空间, 因为只有公共边上才有公共节点, 两个三次多项式在四个点上取之相等, 从而必恒等.

奇妙正方形只说这些就够了. 现在:

您可以用这次的更一般化的方法, 重新验证 \tilde Q_2 的情形;
您可以作出任意 \tilde Q_k 的奇妙正方形. 具体地说, 您掌握的已经足够建立下面的

定理: 对 k \in \mathbb N, 为了在正方形元上使用 Q_k 型的多项式, 有 (k-1)^2 个恒等式 \psi_i 使您 "消去" 正方形内部的 (k-1)^2 个节点, 这时, 子空间 \tilde Q_k = \{\tilde p \in Q_k: \psi_i = 0,\,i = 1,\cdots, (k-1)^2\} 的维数是 4k. 在边界上的所有节点的 \tilde p 值的集合是自由度的集合, 每一个多项式 \tilde p \in \tilde Q_k 由这集合唯一确定.