Numerical Solutions of PDE by FEM

# **Numerical Solutions of** 
# **Partial Differential Equations by** 
# **Finite Element Method**

</div>
</div>

# **Ch. 1, Sec. 2**
# FEM for the model problem 
# with piecewise linear functions

</div>
</div>

## 回顾

在 1.1 节中，我们已经讨论了如下的模型问题：

$$\begin{aligned}
-u''(x) = f(x)&, \quad x \in (0,1), \\\\
u(0) = 0&, \quad u(1) = 0.
\end{aligned}$$

在接下来的部分中，我们将构造有限的子空间来近似该问题的解（Ch. 1, Sec. 2, 本节），并分析所得近似解的误差（Ch. 1, Sec. 3）。

## 分区与节点

我们将区间 $[0,1]$ 分成 $M+1$ 个子区间：

$$0 = x_0 < x_1 < x_2 < \cdots < x_{M+1} = 1,$$

其中 $x_i$ 称为节点，区间为$I_{j}=(x_{j-1}, x_j)$，$h_j = x_j - x_{j-1}$ 称为第 $j$ 个子区间的长度。

定义最大区间长度为
$$h = \max_{1 \leq i \leq M} h_i.$$

我们使用这个参数 $h$ 来衡量分区的精细程度的。

## 分段线性函数空间

在1.1节中，我们定义了如下的函数空间 $V$：

> $V=${$v: v$在$[0,1]$上是连续的函数，$v'$在$[0,1]$上是分段连续且有界的，并且$v(0)=v(1)=0$}

</div>
<div class="col">

我们现在构造一个有限维子空间 $V_h$，它由分段线性函数组成，并且是上述空间 $V$ 的子空间。$V_h$ 为如下的函数集合：在每个子区间 $I_j$ 上，$v$ 是线性的；在区间 $[0,1]$ 上，$v$ 是连续的；并且满足边界条件 $v(0) = v(1) = 0$。

这就是满足给定边界条件的分段线性函数空间。如果去除边界条件，则可以得到在 $[0,1]$ 上的给定分段的分段线性函数空间。

</div>
</div>

## 分段线性函数的基函数

在 $V_h$ 中，我们可以定义如下的基函数 $\\{\phi_i\\}_{i=1}^{M}$ ：

$$\phi_i(x) = \begin{cases}
\frac{x - x_{i-1}}{h_i}, & x \in [x_{i-1}, x_i], \\\\
\frac{x_{i+1} - x}{h_{i+1}}, & x \in [x_i, x_{i+1}], \\\\
0, & \text{otherwise}.
\end{cases}$$

这些基函数 $\phi_i$ 在节点 $x_i$ 处取值为 $1$，在其他节点处取值为 $0$。因此，任意的 $v \in V_h$ 都可以表示为
$$v(x) = \sum_{i=1}^{M} v(x_i) \phi_i(x).$$

## 基函数的简单示例

![](CODE/FEM_01_Raw_Basis.png)

</div>

如图所示，基函数在每个节点处呈现出帽子形状。不相邻的节点处，基函数的值为零。

## FEM 近似解的表示

所以，对于 $v \in V_h$，我们有

$$v(x) = \sum_{i=1}^{M} \eta_{i} \phi_i(x).$$

其中，$$\eta_i = v(x_i).$$

每一个 $v \in V_h$ 都可以唯一地表示为基函数 $\phi_i$ 的线性组合。因此，$V_h$ 是一个维数为 $M$ 的线性空间，其基为 $\{\phi_i\}_{i=1}^{M}$。

## FEM 近似解的表示

![加权基函数](CODE/FEM_02_Weighted_Basis.png)

</div>

## FEM 近似解的表示

![最终结果](CODE/FEM_03_Final_Result.png)

</div>

## FEM 近似解的表示

于是，针对给定边界问题的有限元法可以改写成如下的形式：

> 对于函数空间 $V$ 和 $V_h$，我们有如下的问题：
> $$
> \begin{aligned}
> (M\ )&: Find \ u \in V \ such \ that \ F(u)\le F(v)&, \ &\forall \ v \in V. \\\\
> (M_h)&: Find \ u_h \in V_h \ such \ that \ F(u_h)\le F(v)&, \ &\forall \ v \in V_h.
> \end{aligned}
> $$
>

> 或者等价地，可以写成如下的弱形式：
> $$
> \begin{aligned}
> (V\ )&: Find \ u \in V \ such \ that \ (u',v')=(f,v)&, \ &\forall \ v \in V. \\\\
> (V_h)&: Find \ u_h \in V_h \ such \ that \ (u_h',v')=(f,v)&, \ &\forall \ v \in V_h.
> \end{aligned}
> $$
>

## FEM 近似解的表示

我们接下来针对 $(V_h)$ 问题进行讨论。如之前所示，由于 $V_h$ 的维数为 $M$，我们可以将 $u_h$ 表示为基函数的线性组合：
$$u_h(x) = \sum_{i=1}^{M} \eta_i \phi_i(x).$$
将 $u_h$ 代入 $(V_h)$ 问题中，并选择测试函数 $v$ 为每一个基函数 $\phi_j$，我们可以得到如下的线性方程组：
$$\sum_{i=1}^{M} \eta_i (\phi_i', \phi_j') = (f, \phi_j), \quad j = 1, 2, \ldots, M.$$

## FEM 近似解的表示

于是这个线性方程组可以写成矩阵形式为
$$A \eta = b,$$
其中矩阵 $A$ 的元素为 $A_{ji} = (\phi_i',   \phi_j')$，向量 $b$ 的元素为 $b_j = (f, \phi_j)$，向量 $\eta$ 包含了未知系数 $\eta_i$。

具体地说，线性方程组 $A\eta=b$ 形如：
$$
A\eta=\begin{bmatrix}
(\phi_1', \phi_1') & (\phi_2', \phi_1') & \cdots & (\phi_M', \phi_1') \\\\
(\phi_1', \phi_2') & (\phi_2', \phi_2') & \cdots & (\phi_M', \phi_2') \\\\
\vdots & \vdots & \ddots & \vdots \\\\
(\phi_1', \phi_M') & (\phi_2', \phi_M') & \cdots & (\phi_M', \phi_M')
\end{bmatrix}
\begin{bmatrix}
\eta_1 \\\\
\eta_2 \\\\
\vdots \\\\
\eta_M
\end{bmatrix}=
\begin{bmatrix}
(f, \phi_1) \\\\
(f, \phi_2) \\\\
\vdots \\\\
(f, \phi_M)
\end{bmatrix}
=b.
$$

在多维情形下，FEM 近似解的表示也是类似的。
$A_{ij}=(\nabla\phi_i, \nabla\phi_j)$，$b_j=(\phi_j,f)$，并且 $u_h$ 仍然可以表示为基函数的线性组合。

## 矩阵 $A$ 的结构

对于这个给定的问题，矩阵 $A$ 是三对角矩阵的，也就是说，只有主对角线和相邻的两条对角线上的元素是非零的。这是因为基函数 $\phi_i$ 和 $\phi_j$ 只有在它们对应的节点 $x_i$ 和 $x_j$ 相邻时，它们的导数才会有重叠的支持。而且，矩阵 $A$ 是对称的。

对于这个给定的问题，

$$
A_{ij} = \begin{cases}
\frac{1}{h_i} + \frac{1}{h_{i+1}}, & i = j, \\\\
-\frac{1}{h_{i+1}}, & j = i + 1, \\\\
-\frac{1}{h_i}, & j = i - 1, \\\\
0, & \text{otherwise}.
\end{cases}
$$

</div>
<div class="col">

均匀划分的情况下，矩阵 $A$ 的形式为：
$$
A = \frac{1}{h}\begin{bmatrix}
2 & -1 & 0 & \cdots & 0 \\\\
-1 & 2 & -1 & \cdots & 0 \\\\
0 & -1 & 2 & \cdots & 0 \\\\
\vdots & \vdots & \vdots & \ddots & \vdots \\\\
0 & 0 & 0 & \cdots & 2
\end{bmatrix}.
$$

再利用 $b_j = (f, \phi_j)$，我们就可以组装出完整的线性方程组 $A\eta = b$，并最终求出数值解。

</div>
</div>

<span style="color: blue;">在下一节中，我们将讨论有限元法近似解的误差问题。</span>

> - 在给定的 $V_h$ 中，我们可不可以找到更好的近似解？
> - 有限元法的近似解与真实解之间的误差有多大？
> - 随着分区的细化，误差会如何变化？
> - $\cdots$

</div>
</div>

# **Ch. 1, Sec. 3**
# An error estimate for FEM
# for the model problem

</div>
</div>

## 回顾

在上一节中，我们讨论了如下的模型问题：

$$\begin{aligned}
-u''(x) = f(x)&, \quad x \in (0,1), \\\\
u(0) = 0&, \quad u(1) = 0.
\end{aligned}$$

并给出了该问题的有限元法近似解的表示方法，即
$$u_h(x) = \sum_{i=1}^{M} \eta_i \phi_i(x),$$
其中 $\phi_i$ 是分段线性基函数，系数 $\eta_i$ 通过解线性方程组 $A\eta = b$ 得到。

> 在本节中，我们将分析有限元法近似解 $u_h$ 在其空间 $V_h$ 中的最佳性，并给出误差估计。

## 有限元法近似解的最佳性

我们已经知道，有限元法近似解 $u_h$ 满足如下的弱形式：

$$(u_h', v') = (f, v), \quad \forall v \in V_h.$$

同时，真实解 $u$ 满足：

$$(u', v') = (f, v), \quad \forall v \in V.$$

因此，对于任意的 $v \in V_h$，我们有：

$$(u' - u_h', v') = 0.$$

这表明误差 $u - u_h$ 在 $V_h$ 中是正交于所有测试函数 $v$ 的。

## 在给定函数空间中的最佳性

> 对于任意的 $v \in V_h$，我们有：$ \|\|(u-u_{h})'\|\|\leq \|\|(u-v)'\|\|. $
>

**证明：**

$$
\begin{aligned}
\|\|(u - v)'\|\|^2 &= ( (u - v)', (u - v)' ) \\\\
&= ( (u - u_h + u_h - v)', (u - u_h + u_h - v)' ) \\\\
&= \|\|(u - u_h)'\|\|^2 + 2( (u - u_h)', (u_h - v)' ) + \|\|(u_h - v)'\|\|^2 \\\\
&= \|\|(u - u_h)'\|\|^2 + \|\|(u_h - v)'\|\|^2 \\\\
&\geq \|\|(u - u_h)'\|\|^2.
\end{aligned}
$$

于是，我们得到了有限元法近似解在其空间 $V_h$ 中的最佳性。也就是说在 $V_h$ 中，没有其他函数能比 $u_h$ 更接近真实解 $u$。

## 与真解之间的误差分析

通过估计 $\|\|(u - v)'\|\|$，我们可以得到有限元法近似解 $u_h$ 与真实解 $u$ 之间的误差估计。我们选择$u_{h} \in V_h$为 $u$ 的插值函数，即在每个节点 $x_i$ 处，$u_h(x_i) = u(x_i)$。

如果我们采用这样的线性差值函数 $u_h$，则可以证明如下的误差估计：

> 对于给定的 $[0,1]$ 区间划分，我们有这样的误差估计：
> 
> - $\displaystyle\|u'(x)-u_{h}'(x)\| \leq h\max_{0\leq x \leq 1}\|u'(x)\|,$
> - $\displaystyle\|u(x)-u_{h}(x)\| \leq \frac{h^{2}}{8}\max_{0\leq x \leq 1}\|u''(x)\|,$
>
> 证明在数学分析课程中已经给出了。
>

于是可以推出 $\displaystyle\|\|u'(x)-u_{h}'(x)\|\| \leq h\max_{0\leq x \leq 1}\|u''(x)\|$ .

## 与真解之间的误差分析

结合有限元法近似解的最佳性，我们可以得到如下的误差估计：

> 对于有限元法近似解 $u_h$，我们有如下的误差估计：
> - $\displaystyle\|u'(x)-u_{h}'(x)\| \leq h\max_{0\leq x \leq 1}\|u''(x)\|,$
> - $\displaystyle\|u(x)-u_{h}(x)\| \leq \frac{h^{2}}{8}\max_{0\leq x \leq 1}\|u''(x)\|,$
> - $\displaystyle\|\|u'(x)-u_{h}'(x)\|\| \leq h\max_{0\leq x \leq 1}\|u''(x)\|.$
> 
> 这表明随着分区的细化（即 $h$ 的减小），有限元法近似解 $u_h$ 会更接近真实解 $u$。
>

这给我们提供了有限元法在数值求解偏微分方程中的有效性和收敛性的理论基础。我们并不需要解出真实解 $u$，就可以通过误差估计来评估有限元法近似解的质量。

# **Ch. 1, Sec. 4**
# FEM for the Possion equation

</div>
</div>

## Possion 方程简介

物理学和力学中的许多问题都可以用 Possion 方程来建模；

例如，可以表示温度、电磁势，或者在给定横向载荷强度的作用下、边界固定的弹性膜的位移。

> 在二维区域 $\Omega$ 上，一类特殊边界条件的 Possion 方程的一般形式为：
> $$\begin{aligned}
> -\Delta u(x,y) = f(x,y), &\quad (x,y) \in \Omega, \\\\
> u(x,y) = 0, &\quad (x,y) \in \partial \Omega.
> \end{aligned}$$
> 其中，$\Delta$ 是拉普拉斯算子，定义为 $\displaystyle\Delta = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$。
>

我们可以看到， Possion 方程的形式比较强（$-\Delta u=f$）为了离散化地解决问题，我们首先需要将其转变为更弱的形式，换句话说，就是削弱二阶导数的要求。

## 弱形式、Green 公式与分部积分

对于给定的一个测试函数 $v$，我们将 Possion 方程与 $v$ 相乘并在区域 $\Omega$ 上积分：
$$\int_{\Omega} -\Delta u \cdot v \, d\Omega = \int_{\Omega} f \cdot v \, d\Omega.$$
通过应用 Green 公式（二维情形下的分部积分），我们可以将左侧的二阶导数项转换为一阶导数项：

$$\int_{\Omega} \nabla u \cdot \nabla v \, d\Omega - \int_{\partial \Omega} \frac{\partial u}{\partial n} v \, dS = \int_{\Omega} f \cdot v \, d\Omega.$$

由于边界条件 $u = 0$ 在 $\partial \Omega$ 上成立，测试函数 $v$ 也满足 $v = 0$ 在 $\partial \Omega$ 上，因此边界项消失了。最终，我们得到 Possion 方程的弱形式：

$$\int_{\Omega} \nabla u \cdot \nabla v \, d\Omega = \int_{\Omega} f \cdot v \, d\Omega, \quad \forall v \in V.$$

## 有限元法的离散化

在 `Ch.1, Sec.2` 中，我们已经介绍了如何为一维问题构造有限元空间 $V_h$。对于二维问题，我们可以采用类似的方法构造一个函数空间 $V_h$，该空间由在区域 $\Omega$ 上的单元上的分段线性函数组成.

我们利用这些基函数来表达在这个选定空间中的最优近似解 $u_h$：

$$u_h(x,y) = \sum_{i=1}^{N} \eta_i \phi_i(x,y),$$

代入弱形式，并选择测试函数 $v$ 为每一个基函数 $\phi_j$，我们可以得到如下的线性方程组：

$$\sum_{i=1}^{N} \eta_i (\nabla \phi_i, \nabla \phi_j) = (f, \phi_j), \quad j = 1, 2, \ldots, N.$$

## 组装线性方程组

于是这个线性方程组可以写成矩阵形式为
$$A \eta = b,$$
其中矩阵 $A$ 的元素为 $A_{ji} = (\nabla \phi_i, \nabla \phi_j)$，向量 $b$ 的元素为 $b_j = (f, \phi_j)$，向量 $\eta$ 包含了未知系数 $\eta_i$。

也就是说，线性方程组 $A\eta=b$ 形如：
$$
\underbrace{\begin{bmatrix}
(\nabla \phi_1, \nabla \phi_1) & (\nabla \phi_2, \nabla \phi_1) & \cdots & (\nabla \phi_N, \nabla \phi_1) \\\\
(\nabla \phi_1, \nabla \phi_2) & (\nabla \phi_2, \nabla \phi_2) & \cdots & (\nabla \phi_N, \nabla \phi_2) \\\\
\vdots & \vdots & \ddots & \vdots \\\\
(\nabla \phi_1, \nabla \phi_N) & (\nabla \phi_2, \nabla \phi_N) & \cdots & (\nabla \phi_N, \nabla \phi_N)
\end{bmatrix}}_{\text{刚度矩阵 (Stiffness Matrix)A}}
\begin{bmatrix}
\eta_1 \\\\
\eta_2 \\\\
\vdots \\\\
\eta_N
\end{bmatrix}=
\begin{bmatrix}
(f, \phi_1) \\\\
(f, \phi_2) \\\\
\vdots \\\\
(f, \phi_N)
\end{bmatrix}
$$

## 其余参考

其余部分内容请参考：

[点我](/SLIDES/FECK/index.html)

</div>

# **Ch. 1, Sec. 5**
# The Hilbert spaces 
# $L_2(\Omega)$, $H^1(\Omega)$ and $H_0^1(\Omega)$

</div>
</div>

## 回顾

在前面的章节中，我们已经通过一些具体的一维和二维问题，看到了有限元方法是如何通过分片多项式函数来逼近微分方程的解的。当时我们隐含地做了一个假设：解是连续的，甚至在每个小区间上有连续的导数。从数值构造的角度看，这样的假设并不会立刻出问题，但如果我们希望从数学上真正理解“解是什么”“为什么解存在”“有限元解到底在逼近什么”，那么继续停留在“连续函数、连续导数”的层面就显得非常不够了。第 1.5 小节的全部目的，其实就是为有限元方法寻找一个正确而自然的函数空间背景，以便我们能够严谨地讨论有限元方法的收敛性和误差估计。

## 内积空间简介

我们称 $L$ 为 $V$ 上的线性范数，如果对于所有的 $u, v \in V$ 和所有的标量 $a, b \in \mathbb{R}$，都有
$$L(au + bv) = aL(u) + bL(v).$$

给定一个二元函数 $a(\cdot, \cdot): V \times V \to \mathbb{R}$，如果对于所有的 $u, v, w \in V$ 和所有的标量 $a, b \in \mathbb{R}$，都有
1. **对称性：**$a(u, v) = a(v, u)$,
2. **线性：**$a(au + bv, w) = a a(u, w) + b a(v, w)$, 
3. **正定性：**$a(u, u) \geq 0$，且当且仅当 $u = 0$ 时取等号，

那么我们称 $a(\cdot, \cdot)$ 是 $V$ 上的一个**内积**。

由内积，我们可以定义对应的范数 $\|\|u\|\| = \sqrt{a(u, u)}$。

## 内积空间简介

我们进一步回顾，如果 $V$ 是一个带有标量积及对应范数 $\|\|\cdot\|\|$ 的线性空间，且 $V$ 是完备的，即如果关于 $\|\|\cdot\|\|$ 的每一个柯西列都是收敛的，那么 $V$ 被称为一个希尔伯特空间（**H**ilbert space）。

对于空间 $V$ 中具有范数 $\|\|\cdot\|\|$ 的元素 $v_i$ 组成的序列 $v_1, v_2, v_3, \dots$，如果对于任意 $\varepsilon > 0$，都存在一个自然数 $N$，使得当 $i, j > N$ 时 $\|\|v_i - v_j\|\| < \varepsilon$，则称该序列为柯西序列。此外，如果当 $i \to \infty$ 时 $\|\|v - v_i\|\| \to 0$，则称 $v_i$ **收敛**于 $v$。

> - 我们可以将希尔伯特空间视为一个带有内积的完备空间。
> - 希尔伯特空间的完备性确保了在该空间中，所有“收敛趋势合理的函数序列”都能找到极限。
>
> 对于分段函数空间，这一点在数值分析中极其关键，因为有限元方法本质上是构造一列近似函数并研究其极限行为。

## $L^2$ 空间

我们先从最简单、也是最基础的空间 $L^2$ 开始。

对于一个开区间 $\Omega \subset \mathbb{R}^n$，我们定义 $L^2(\Omega)$ 为如下的函数集合：

$$L^2(\Omega) = \{v: \Omega \to \mathbb{R} \ | \ v \text{ 是可测函数且 } \int_{\Omega} v^2 \, d\Omega < \infty \}.$$

我们在 $L^2(\Omega)$ 上定义内积为：
$$(v, w) = \int_{\Omega} v(x) w(x) \, d\Omega.$$

由此，我们可以定义对应的范数为：
$$\|\|v\|\| = \sqrt{(v, v)} = \left(\int_{\Omega} v^2 \, d\Omega\right)^{\frac{1}{2}}.$$

## $L^2$ 空间

> 直观地说，$L^2$空间描述的是这样一类函数：它们不一定连续，不一定有导数，甚至可以在某些点上变得很大以至于到达 $\infty$，但只要“平方以后在整体上是（勒贝格）可积的”，我们就认为它们是可以接受的。

然而，$L^2$仍然不足以描述偏微分方程。原因很简单：微分方程中出现了导数，而$L^2$并不要求导数存在。在此基础上，我们可以进一步定义几个 $Sobolev$ 空间:

>  $H^1(\Omega)$ 空间：$H^1(\Omega) = \{v \in L^2(\Omega) \ | \ \nabla v \in (L^2(\Omega))^n \}$，即函数 $v$ 及其一阶弱导数都属于 $L^2(\Omega)$。
> 
>  $H_0^1(\Omega)$ 空间：$H_0^1(\Omega) = \{v \in H^1(\Omega) \ | \ v = 0 \text{ 在 } \partial \Omega\}$，即函数 $v$ 属于 $H^1(\Omega)$ 且在边界上取零值。

## Sobolev 空间

待补充

## 本节完

# 敬请老师和同学们批评指正！

<hr/>

参考资料：

[Numerical Solution of Partial Differential Equations<br>by the Finite Element Method, Claes Johnson](https://books.google.com/books?id=0IFCAwAAQBAJ&lpg=PP1&ots=q3tIHf9A8I&dq=info%3AzjpNBiIw5coJ%3Ascholar.google.com&lr&pg=PP1#v=onepage&q&f=false)

我的部分联系方式：

[hsingyu.yen@icloud.com](mailto:hsingyu.yen@icloud.com)

</div>
</div>

## 参考书目

Numerical Solution of Partial Differential Equations by the Finite Element Method, Claes Johnson: