2021-5-18 有限元方法从0开始学习笔记 第二天 (一维线性元简介, 与差分格式的联系)

又是新的一天, 继续开始有限元从0开始学习, 参考资料
Brenner S, Scott R. The mathematical theory of finite element methods[M]. Springer Science & Business Media, 2007.

昨天内容回顾

2021-5-17 有限元方法从0开始学习笔记 第一天

• 利用模型问题, 构造了其变分形式, 并且验证了当 $f$ 和 $u$ 的光滑性足够时, 变分问题能够得到边值问题对古典解.
• 阐述了 Ritz-Galerkin 框架的逼近过程, 并给出了简单的离散系统的解的存在唯一性的证明.(到后面可以直接使用 Lax-Milgram 定理).
• 给出了一维模型问题 Ritz-Galerkin 框架下的能量误差以及 $L^2$ 范数误差估计的示范.

下面开始今天内容的学习.

0.4 分片多项式空间-有限元方法

前面介绍 Ritz-Galerkin 框架时, 没有给出子空间 $S$ 的具体形式, 实际上对于选取不同的子空间, 可以给出不同的方法(谱方法, 有限元方法) 等, 这里主要考虑有限元方法, 在本章会给出第一个有限元基函数空间, 并且定义简单的插值算子, 以及给出插值误差估计.
继续考虑模型问题
$$\{\begin{array}{l}-\frac{d^{2}u}{d{x}^2}=f\in (0,1),\\ u(0)=0,u^{{}^{\prime }}(1)=0.\end{array}.$$
对应的弱形式
$$\text{find}u\in Vs.t.,a(u,v)=(f,v)\forall v\in V.(\text{WP})$$
以及 Ritz-Galerkin 方法的离散问题
$$\text{find}{u}_s\in S\subset Vs.t.,a(u_s,v)=(f,v)\forall v\in S\subset V.(\text{R-G})$$
其中, 空间 $V$ 定义如下:
$$V=H_{\diamond }^1(0,1)=\left\{u\in H^1(0,1)|u(0)=0\right\}.$$
前面只是抽像地给出了问题在子空间 $S$ 下的计算, 下面定义第一个有限元空间:
记 $0=x_0<x_1<\cdots <x_n=1$ 为区间 $[0,1]$ 上的剖分, 定义空间 $S$ 为满足如下条件的函数 $v$ 构成的线性空间:

• $v\in C^0[0,1]$,
• ${v|}_{[x_{i-1},x_i]}$ 为线性多项式,
• $v(0)=0$

下面自然的会有两个问题

• 上面定义的 $S$ 是不是 $V$ 的子空间,
• 如何找出 $S$ 的基函数, 应该取什么样的基函数,

这里先回答第二个问题, 对于第一个问题留在后面回答.
定义 $S$ 的基函列 { ϕ i } i = 1 n \left\lbrace \phi_{i} \right\rbrace_{i = 1}^{n} {ϕi​}i=1n​, 满足 ${\varphi }_i(x_j)={\delta }_{ij}$ , 这个性质称为插值性质.
引理: 满足上面条件的函数序列存在唯一, 并且构成子空间 $S$ 的基.
注1 称上面的 { ϕ i } \left\lbrace \phi_{i} \right\rbrace {ϕi​} 为空间 $S$ 的节点基函数, 称 { v ( x i ) } \left\lbrace v(x_{i}) \right\rbrace {v(xi​)} 为函数 $v$ 在节点处的值.
定义 对 $v\in C^0([0,1])$, 定义插值算子 $I:C^0([0,1])\to S$
$$v_I(x):=Iv(x)=\sum _{i=1}^{n}v(x_i)\varphi (x).$$
由上面插值算子的定义, 可以知道, 显然 ${\varphi }_i$ 可以张成整个空间 $S$, 根据基函数的插值性质, 容易验证 { ϕ i } i = 1 n \left\lbrace \phi_{i} \right\rbrace_{i=1}^{n} {ϕi​}i=1n​ 是线性无关的.
引理: 设 $v\in S$, 则 $v_I=v$.
证明, 只要验证 $\forall x\in [0,1],v_I(x)=v(x)$ 即可.显然, 有
$$\begin{gathered} v_I(x_{i-1})-v(x_{i-1})=0 \\ {v}_I(x_i)-v(x_i)=0 \end{gathered}$$
注意到 $v_I(x),v(x)$ 限制在每个区间上都是分片多项式. 因此 $v_I(x)=v(x)\text{in}[x_{i-1},x_i]\forall i$, 从而 $v_I(x)=v(x)$,
下面给出插值误差在能量范数下的误差估计.
定理: 设 $h={\max }_{1\le i\le n}\left(x_i-x_{i-1}\right)$ ,则有如下能量范数估计:
$$\Vert u-u_I{\Vert }_E\le Ch\Vert u^{{}^{\prime \prime }}\Vert$$
$\forall u\in V$, 并且常数 $C$ 不依赖于网格大小 $h$ 以及 $u$.

证明, 利用 scaling 技巧, 只要想办法证明上面的估计分片成立即可, 即想办法做出如下估计:
$${\int }_{x_{j-1}}^{x_j}{\left((u-u_I{)}^{\prime }(x)\right)}^{2}dx\le c{\left(x_j-x_{j-1}\right)}^2{\int }_{x_{j-1}}^{x_j}{u}^{{}^{\prime \prime }}(x{)}^{2}dx$$
记误差函数 $e(x)=(u-u_I)(x)$, 注意 $u_I(x)$ 是分片线性多项式, 从而上述估计等价于证明:
$${\int }_{x_{j-1}}^{x_j}{e}^{{}^{\prime }}(x{)}^{2}dx\le c{\left(x_j-x_{j-1}\right)}^2{\int }_{x_{j-1}}^{x_j}{e}^{{}^{\prime \prime }}(x{)}^{2}dx$$
为了使得估计不依赖于网格, 利用 scaling 技术, 将上述估计用仿射变换变换到参考单元上, 记
$$\hat{x}=\frac{x-x_{j-1}}{x_j-x_{j-1}}$$
再有
$$\begin{array}{rl}& x=x_{j-1}+\hat{x}(x_j-x_{j-1})\\ & dx=\left(x_j-x_{j-1}\right)d\hat{x}\\ & \hat{e}(\hat{x})=e(x_{j-1}+\hat{x}(x_j-x_{j-1}))=e(x)\\ & \frac{d\hat{e}(\hat{x})}{d\hat{x}}=\frac{de(x)}{dx}\frac{dx}{d\hat{x}}=(x_j-x_{j-1})\frac{de(x)}{dx}\\ & \frac{d{\hat{e}}^{{}^{\prime }}(\hat{x})}{d\hat{x}}=\frac{d{\hat{e}}^{{}^{\prime }}(\hat{x})}{dx}\frac{dx}{d\hat{x}}=(x_j-x_{j-1}{)}^2\frac{d{e}^{{}^{\prime }}(x)}{dx}\end{array}$$
代入上式, 从而只要在参考单元上证明:
$${\int }_0^1{\hat{e}}^{{}^{\prime }}(\hat{x}{)}^{2}d\hat{x}\le c{\int }_0^1{\hat{e}}^{{}^{\prime \prime }}(\hat{x}{)}^{2}d\hat{x}$$
下面只要证明这一点即可, 为了简化记号, 记 $w(x)=\hat{e}(\hat{x})$, 下面证明
$${\int }_0^1{w}^{{}^{\prime }}(x{)}^{2}dx\le c{\int }_0^1{w}^{{}^{\prime \prime }}(x{)}^{2}dx$$
左边是一阶导数, 右边是二阶导数, 显然应该利用 Newton-Leibniz 公式建立两边的联系, 注意到 $w(x)$ 在端点上为 0, 利用 Rolle 定理可知, 存在 $\xi \in (0,1),s.t.,w^{{}^{\prime }}(\xi )=0$, 从而有
$$\begin{array}{rl}|w^{{}^{\prime }}(y)|& =\left|w^{{}^{\prime }}(y)-w^{{}^{\prime }}(\xi )\right|=\left|{\int }_{\xi }^y{w}^{{}^{\prime \prime }}(x)dx\right|\\ & \le {\left|{\int }_{\xi }^y{1}^{2}dx\right|}^{1/2}{\left|{\int }_{\xi }^y{w}^{{}^{\prime \prime }}(x{)}^{2}dx\right|}^{1/2}\\ & =|y-\xi {|}^{1/2}{\left({\int }_0^1{w}^{{}^{\prime \prime }}(x{)}^{2}dx\right)}^{1/2}\end{array}$$
从而有:
$$\begin{array}{rl}{\int }_0^1{w}^{{}^{\prime }}(y{)}^{2}dy& \le {\int }_0^1\left|y-\xi \right|dy{\int }_0^1{w}^{{}^{\prime \prime }}(x{)}^{2}dx\\ & \le c{\int }_0^1{w}^{{}^{\prime \prime }}(x{)}^{2}dx\end{array}$$
其中
$$c=\underset{0<\xi <1}{\sup }{\int }_0^1|y-\xi |dy=\frac{1}{2}$$
从而证明了
$${\int }_{x_{j-1}}^{x_j}{e}^{{}^{\prime }}(x{)}^{2}dx\le c{\left(x_j-x_{j-1}\right)}^2{\int }_{x_{j-1}}^{x_j}{e}^{{}^{\prime \prime }}(x{)}^{2}dx$$
上式两边对 $j$ 从 0 到 n 进行求和:
$${\int }_0^1{e}^{{}^{\prime }}(x{)}^{2}dx\le c{h}^2{\int }_0^1{e}^{{}^{\prime \prime }}(x{)}^{2}dx$$
从而有
$$\Vert u-u_I{\Vert }_E\le Ch\Vert u^{{}^{\prime \prime }}\Vert$$
推论: $\Vert u-u_S\Vert +Ch\Vert u-u_S{\Vert }_E\le 2(Ch{)}^2\Vert u^{{}^{\prime \prime }}\Vert$.
根据前面对 Ritz-Galerkin 框架的误差估计, 马上可以得到.
在后面的内容, 仍然会继续考虑类似的误差估计, 到时候就是高维情形, 高维情形会和一维情形在细节上会有比较大的差距,但是就整体思路而言, 我认为基本上是一样的.

0.5 和差分方法的关系

下面计算出由上面基函数所定义的线性系统的刚度矩阵, 由 ${\mathbf{K}}_{ij}=a({\varphi }_j,{\varphi }_i)$, 注意这里的双线性型是对称的, 这将直接导致刚度矩阵也是对称的, 对 $1\le i\le n-1$:
$$\begin{array}{rl}& K_{i,i}={\int }_{x_{i-1}}^{x_i}{\varphi }_i^{{}^{\prime }}(x){\varphi }_i^{{}^{\prime }}(x)+{\int }_{x_i}^{x_{i+1}}{\varphi }_{i+1}^{{}^{\prime }}(x){\varphi }_{i+1}^{{}^{\prime }}(x)=\frac{1}{h_i}+\frac{1}{h_{i+1}}\\ & K_{i,i+1}=K_{i+1,i}=a({\varphi }_i,{\varphi }_{i+1})={\int }_{x_i}^{x_{i+1}}{\varphi }_i^{{}^{\prime }}(x){\varphi }_{i+1}^{{}^{\prime }}(x)=-\frac{1}{h_{i+1}}\end{array}$$

此外, $K_{nn}=\frac{1}{h_n}$不难知道, 对于其它位置, 刚度矩阵都是 0.
对右端项用如下逼近(将 $f(x)$ 在 $x_i$ 处进行 Taylor 展开, 之后对基函数直接积分):
$$\begin{array}{rl}{\int }_0^{1}f(x){\varphi }_i(x)dx& ={\int }_{x_{i-1}}^{x_i}f(x){\varphi }_i(x)dx+{\int }_{x_i}^{x_{i+1}}f(x){\varphi }_i(x)dx\\ & =(f(x_i)+O(h))\left({\int }_{x_{i-1}}^{x_i}{\varphi }_i(x)dx+{\int }_{x_i}^{x_{i+1}}{\varphi }_i(x)dx\right)\\ & =(f(x_i)+O(h))\left({\int }_{x_{i-1}}^{x_i}\frac{1}{h_i}(x-x_{i-1})dx+{\int }_{x_i}^{x_{i+1}}\frac{-1}{h_{i+1}}(x-x_{i+1})dx\right)\\ & =\frac{1}{2}(h_i+h_{i+1})(f(x_i)+O(h)).\end{array}$$
其中 h = max ⁡ { h i , h i + 1 } h = \max\{h_{i}, h_{i+1}\} h=max{hi​,hi+1​}, 利用 Taylor 展开马上就可以验证, 如果想要让逼近精度达到二阶, 则网格步长需要满足条件, 即 $1-(h_i/h_{i+1})=O(h)$. 考虑方程组 $\mathbf{K}\mathbf{U}=\mathbf{F}$ 的第 $i$ 个分量方程
$$K_{i,i-1}{U}_{i-1}+K_{i,i}{U}_i+K_{i,i+1}{U}_{i+1}=\left(f,{\varphi }_i\right)$$
代入前面的计算:
$$-\frac{1}{h_i}{U}_{i-1}+\left(\frac{1}{h_i}+\frac{1}{h_{i+1}}\right)U_i-\frac{1}{h_{i+1}}{U}_{i+1}=\frac{1}{2}(h_i+h_{i+1})(f(x_i)+O(h))$$
整理为如下形式得:
$$\frac{-2}{h_i+h_{i+1}}\left[\frac{U_{i+1}-U_i}{h_{i+1}}-\frac{U_i-U_{i-1}}{h_i}\right]=f(x_i)+O(h)$$
注意, 当网格取一直剖分的时候, 上述格式退化为中心差分格式, 并且由上面条件可知, 截断误差变为 $O(h^2)$, 即:
$$-\frac{U_{i+1}-2U_i+U_{i-1}}{h^2}=f(x_i)+O(h^2).$$
关于其它注解, 可以看参考资料 第十页.
今天就到这里啦, 明天有课, 可能会停一天,或者少打一点 后面希望继续坚持.