上期回顾
程序
例题 1
利用拉格朗日多项式做基函数编程实现: 设函数 $f(x)=10(x-1)^2-1$, 在线性函数空间中找到最佳逼近函数 $u(x)$, 求解域 $\Omega = [1, 2]$, 两个真解点为: $x_0 = 1+1/3, x_1 = 1+2/3$.
主程序
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]$ 内剖分为两个单元:
\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 因此,
节点信息用向量表示nodes = [0, 0.2, 0.4, 0.7, 1];每个单元由三个节点组成
局部编号 0 1 2 $\Omega^{(0)}$ 0 1 2 $\Omega^{(1)}$ 2 3 4 因此,
单元信息用矩阵信息elements= [0, 1, 2; 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). $$
利用计算机绘图得到
小结
有限元基函数可以保证最小二乘法的系数矩阵是稀疏的, 而且基函数是简单的多项式.
下节预告
有限元计算在实际中不可能是手算的, 所以下次介绍一下编程技巧.