有限元简述(5)：有限元基函数与网格划分

上期回顾

程序

例题 1
利用拉格朗日多项式做基函数编程实现: 设函数 $f(x)=10(x-1)^2-1$ , 在线性函数空间中找到最佳逼近函数 $u(x)$ , 求解域 $\Omega = [1, 2]$ , 两个真解点为: $x_0 = 1+1/3, x_1 = 1+2/3$ .

主程序

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
f = @(x)(10*(x-1).^2 - 1);
syms x
psi = [1, x];
points = [1+1/3, 2-1/3];

[ A, b ] = interpolation( f, psi, points );
c = A\b;
u = psi*c;

% 可视化
omega = [1, 2];
X = omega(1) : 0.01 : omega(2);
U = inline(vectorize(u), 'x');
Plot(X, f, U)

% 结果展示
result(A, b, c, u)

例题 2
假设 $f(x)=\sin(2\pi x), \Omega=[0, 1]$ , 利用拉格朗日多项式做基函数, 分别用最小二乘法和插值法寻找最佳逼近函数 $u(x)$ .
注: 可假设基函数个数为 4 个, 拉格朗日点为等距节点

最小二乘法

插值法

程序过长, 源代码见原文链接. 肉眼可见, 最小二乘法逼近效果更好, 但其系数矩阵更复杂一些.

小结

上次推送了重要的拉格朗日多项式作基函数, 结合上述例题, 知用于插值法中, 可使得系数矩阵 $A=I$ 为单位矩阵, 但是, 应用到最小二乘法中, 系数矩阵却不是稀疏矩阵的, 因为

A_{i, j} =\left(\psi_{i}, \psi_{j}\right) \neq 0.

因此, 本次继续讨论逼近函数问题, 考虑 ..有限元基函数.. .

有限元基函数

引子

在之前 3 期推送(逼近函数)中, 基函数 $\psi_i(x), i \in {\mathcal{I}_s}$ 在整个区域 $\Omega$ 上一般都是非零的, 如下图所示.

而有限元法要求基函数有紧支撑, 即它们只有在一小部分区域内是非零的, 在大部分区域都是零. 另外, 一般要求基函数是分段多项式函数, 如下图所示

单元和节点

$\color{gray}{\textit{Elements and nodes}}$

将区域 $\Omega$ 剖分为有限个互补重叠的子区域 $\Omega^{(e)}, e=0,\dots, N_e$ :

\Omega=\Omega^{(0)} \cup \cdots \cup \Omega^{\left(N_{e}\right)}

称这样的子区域为单元, 每个单元有一个区分其他单元的上标 $^{(e)}$ . 然后, 在每个单元 $\Omega^{(e)}$ 上定义点的集合, 命名为节点. 为区分不同节点, 我们还需要一套节点的全局编号, 并且在每个小单元上定义局部编号.

单元和节点完全定义了有限元网格, 即完成了对空间的离散. 上图为常见的均匀剖分网格, 即单元的规格相同, 节点间的距离也相同.

一个例子
$\Omega=[0, 1]$ 内剖分为两个单元:
2 个单元 5 个节点的单位有限元网格
\begin{aligned} \Omega^{(0)}=[0, 0.4]\ \Omega^{(1)}=[0.4, 1] \end{aligned}
每个单元内有 3 个节点. 则有全局编号
全局编号 0 1 2 3 4
坐标 0 0.2 0.4 0.7 1
因此, 节点信息用向量表示
1
nodes = [0, 0.2, 0.4, 0.7, 1];
每个单元由三个节点组成
局部编号 0 1 2
$\Omega^{(0)}$ 0 1 2
$\Omega^{(1)}$ 2 3 4
因此, 单元信息用矩阵信息
1
2
elements= [0, 1, 2;
           2, 3, 4];

全局编号	0	1	2	3	4
坐标	0	0.2	0.4	0.7	1

局部编号	0	1	2
$\Omega^{(0)}$	0	1	2
$\Omega^{(1)}$	2	3	4

基函数

有限元基函数用 $\varphi_i(x)$ 表示, 索引 $i$ 为节点的全局编码, 即基函数的个数等于节点的个数. 在逼近函数问题时, 我们另 $\psi_i(x) = \varphi_i(x)$ .

构造原理

为方便, 在单元 $e$ 中, 令 $i$ 为全局编号, $r$ 为局部编号, 有限元基函数构造遵循以下规则:

$r$ 为单元 $e$ 内部编号:
选取 $\varphi_i$ 在单元 $e$ 内 $r$ 处为 1, 其他节点为 0; 在其他单元上 $\varphi_i(x) \equiv 0$ .
$r$ 为单元 $e$ 边界编号:
选取 $\varphi_i$ 在单元 $e$ 和单元 $e+1$ 上不恒为 0 ( $\varphi_i \not\equiv0$ ), 具体的, 在 $r$ 处等于 1, 其他节点处等于 0. 而在其他单元上, $\varphi_i(x) \equiv 0$ .

下图展示对这两条规则直观的理解

$\varphi_i$ 的性质

帽子

\varphi_{i}\left(x_{j}\right)=\delta_{i j} \quad \delta_{i j}= \begin{cases} {1,} & {i=j} \\\\[3pt] {0,} & {i \neq j} \end{cases}

从而

u\left(x_{i}\right)=\sum_{j \in \mathcal{I}_{s}} c_{j} \psi_{j}\left(x_{i}\right)=\sum_{j \in \mathcal{I}_{s}} c_{j} \varphi_{j}\left(x_{i}\right)=c_{i} \varphi_{i}\left(x_{i}\right)=c_{i}

即..系数.. $c_i$ 的值等于 $u$ 在 $i$ 处的函数值(数值解).

$\varphi_i(x)$ 在区域 $\Omega$ 的大部分区域等于 0.
- $\varphi_i(x)\neq 0$ 仅在有全局节点 $i$ 的单元中成立.
- $\varphi_i(x)\varphi_j(x)\neq 0$ 仅在同时有全局节点 $i$ 和全局节点 $j$ 的单元中成立.

因此, $A*{ij} = \int*{\Omega^{(e)}} \varphi \_i(x) \varphi_j(x) {\rm d} x$ 在系数矩阵中大多等于 0.

假设每个单元有 $d+1$ 个节点, 那么局部拉格朗日多项式的次数是 $d$ .

分段线性有限元基函数

\varphi_{i}(x)= \begin{cases} 0, \quad & x < x_{i-1} \\\\[3pt] (x-x_{i-1})/h, & x_{i-1} \leq x < x_i \\\\[3pt] 1-(x-x_{i-1})/h, & x_{i} \leq x < x_{i+1} \\\\[3pt] 0, & x \geq x_{i+1} \end{cases}

分段二次有限元基函数

在每一个单元上, 都是拉格朗日多项式作基函数

建立线性方程组

组装系数矩阵

[A_{i,j}] = \big(\varphi_j(x), \varphi_i(x)\big).

由基函数的性质, 可知只有 $A_{i, i}$ 或者 $A_{i,i-1}=A_{i.i+1}$ 不为 0, 即 $A$ 为三对角阵.

当 $0<i<N$

\begin{aligned} A_{i,i} & = \int_{\Omega} \varphi_i^2 {\rm d} x \\\\[3pt] & = \int_{i-1}^{i} \varphi_i^2 {\rm d} x + \int_{i}^{i-1} \varphi_i^2 {\rm d} x \\\\[3pt] &=\int_{i-1}^{i} \left(\frac{x-x_{i-1}}{h}\right)^2 {\rm d} x + \int_{i}^{i-1} \left(1-\frac{x-x_{i}}{h}\right)^2 {\rm d} x \\\\[3pt] & = \frac{2h}{3}. \end{aligned}

以及 $0 < i \leq N$

$A_{i,i+1}=A_{i,i-1}$ .最后显然,

A_{0,0}=\int_{\Omega} \varphi_{0}^{2} \mathrm{d} x=\int_{x_{0}}^{x_{1}}\left(1-\frac{x-x_{0}}{h}\right)^{2} \mathrm{d} x=\frac{h}{3}.

$A_{N, N}=A_{0, 0}.$

载入右端项

\begin{aligned} b_{i}&=\int_{\Omega} \varphi_{i}(x) f(x) \mathrm{d} x \\\\[3pt] & = \int_{x_{i-1}}^{x_{i}} \frac{x-x_{i-1}}{h} f(x) \mathrm{d} x+\int_{x_{i}}^{x_{i+1}}\left(1-\frac{x-x_{i}}{h}\right) f(x) \mathrm{d} x \end{aligned}

例题

将 $\Omega = [0, 1]$ 划分为两个相同的的子区域, 设 $f(x)=x(1-x)$ , 利用一阶有限元基函数(P1)计算最佳逼近 $u(x)$ .

利用上述公式, 得到

A=\frac{h}{6}\left[\begin{array}{lll} {2} & {1} & {0} \\\\[3pt] {1} & {4} & {1} \\\\[3pt] {0} & {1} & {2} \end{array}\right], \quad b=\frac{h^{2}}{12}\left[\begin{array}{c} {2-3 h} \\\\[3pt] {12-14 h} \\\\[3pt] {10-17 h} \end{array}\right].

进而得到

c_{0}=\frac{h^{2}}{6}, \quad c_{1}=h-\frac{5}{6} h^{2}, \quad c_{2}=2 h-\frac{23}{6} h^{2}.

所以,

u(x)=c_{0} \varphi_{0}(x)+c_{1} \varphi_{1}(x)+c_{2} \varphi_{2}(x).

利用计算机绘图得到

小结

有限元基函数可以保证最小二乘法的系数矩阵是稀疏的, 而且基函数是简单的多项式.

下节预告

有限元计算在实际中不可能是手算的, 所以下次介绍一下编程技巧.