FE

# 很简单的有限元方法介绍
# A Very Basic Introduction to
# Finite Element Method

</div>
</div>

# Part.1
# 什么是有限元方法（FEM）？

</div>
</div>

## 用于求解复杂的工程和物理问题的计算方法

无论是计算桥梁的应力、飞机的气动外形，还是热传导分布，有限元无处不在。

[![图一](image_1.png)](https://www.comsol.com/multiphysics/finite-element-method)

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

比如，要计算一个形状极其不规则的齿轮在受力后的变形。我们无法给出一个简单的数学公式（**解析解**）来直接算出答案。但是，如果你：

> 1. 把这个齿轮切分成成千上万个极小的、形状规则的小块；
> 2. 算出每个小块的受力情况；
> 3. 组装整合每一个小块的解，就能得到整个齿轮的近似的**数值解**；
>
> 以此可以看出我所以为的有限元方法的基本思想。

</div>
</div>

## 为什么要用有限元方法？

**传统方法的局限性：**

- 现实世界的复杂性： 现实中的汽车底盘、飞机机翼、人体骨骼，其形状是不规则的，材料可能是非线性的，受力情况是动态变化的。
- 为什么不用解析解？比如给定一个经典的 Poisson 方程 $-\Delta u = f$ 。如果可以求出解的数学公式，那非常完美。但通常只能解决几何形状极其简单（如圆柱、立方体）且边界条件简单的问题。
- 在一部分问题中，我们能做的最好结果，就是找到一个足够好的**近似解**。

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

**几个问题**

- 什么样的方法能用有限的数据极好地近似一个函数？
- 如何在不知道真解 $u(x,y)$ 的情况下找到这个近似解？
- 我们已经有了很多手段，比如傅里叶级数、泰勒展开等，为什么还要用有限元方法？

</div>

## 有限元方法的核心思想

有限元方法的一个比较重要的思想是使用**分段线性函数**来近似复杂的解。如果你观察典型 PDE 的解（例如 Stokes 流中的压力分布、冲击波中的密度梯度、......），你会发现它们通常具有以下特征：

- 在大部分区域是平滑的。
- 但在某些微小区域（如角落、界面）变化剧烈。
- 可能存在折点、[奇点（Singularities）](https://en.wikipedia.org/wiki/Singular_point_of_a_curve)。

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

![alt text](image_2.png)

</div>

## 分段线性函数的优势

**为什么不用傅里叶级数？**

傅里叶级数形式为 $\sum U_k \sin(kx)$。它只有在函数全局光滑时效果才好。如果解中有一个尖点，傅里叶级数需要无穷多项才能逼近，且效率极低。

![alt text](image.png)

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

**为什么不用泰勒级数？**

形式为 $\sum U_k (x-x_0)^k$ 。泰勒展开通常只在展开点 $x_0$ 附近准确，远离该点误差巨大。

![alt text](image-3.png)

</div>
</div>

# Part.2
# 从一个例子说起：
# 有限元方法求解 Poisson 方程

</div>
</div>

## The Game is Afoot

我们先从一类特殊的 Poisson 方程入手：

> 在区域 $\Omega$ 上，求解：
> 
> $$
> -\Delta u = f \quad \text{in } \Omega
> $$
> 
> 边界条件：
> 
> $$
>  u = 0 \quad \text{on } \partial \Omega
> $$
>
> 这里，$f$ 是已知函数，$u$ 是我们要求解的函数。
>

我们的目的是找到 $u$ 在 $\Omega$ 上满足边界条件和上述方程的近似解。

## 弱解

**目标：**寻找函数 $u$（来自合适的函数空间，且边界值为零），使得上式对所有合适的测试函数 $\varphi$ 都成立。

这样就得到一个弱解了，因为我们不能直接比较函数，但可以比较对所有测试函数的投影（积分）是否相等。

例如：

$$
\int g = \int h , \quad \int x\cdot g = \int x\cdot h , \quad \int x^2 g = \int x^2 h , \quad \ldots \Rightarrow g = h
$$

如果对于所有的**测试函数 (Test Functions)** $\varphi(x)$，其积分都相等，则认为这是一个弱解。

## 弱形式

如果我们需要求解一个弱解，就将原方程 $-\Delta u = f$ 两边乘以一个测试函数 $\varphi$，并在求解域 $\Omega$ 上积分：

$$
-\int_{\Omega} \varphi \Delta u \, d\mathbf{x} = \int_{\Omega} \varphi f \, d\mathbf{x}
$$

**分部积分：**

对左侧项 $-\int_{\Omega} \varphi \Delta u \, d\mathbf{x}$ 使用分部积分法：

$$
\int_{\Omega} \nabla \varphi \cdot \nabla u \, d\mathbf{x} - \int_{\partial \Omega} \varphi \mathbf{n} \cdot \nabla u \, dS = \int_{\Omega} \varphi f \, d\mathbf{x}
$$

**边界项处理：**

由于原问题 $u$ 满足齐次 Dirichlet 边界条件 ($u=0$)，测试函数 $\varphi$ 也必须满足同样的条件 （$\varphi=0$ on $\partial \Omega$ ）。因此，边界积分项 $\int_{\partial \Omega} \varphi \mathbf{n} \cdot \nabla u \, dS$ 消失。

## 离散化

在计算机上无法找到无限维空间中的精确解 $u$。引入 [Galerkin 方法](#gotojlj)，要求测试函数空间与解函数空间完全相同，把测试函数缩小到有限维子空间中，用给定的基函数表示出解 $u$ 的近似。

**我们将其近似为有限维解 $u_h$：**

将解 $u$ 近似为一组 $N$ 个形函数（Shape Functions）$\{\varphi_{j=0\cdots N-1}\}$ 的线性组合：

$$
u_h(\mathbf{x}) = \sum_{j=0}^{N-1} U_j \varphi_j(\mathbf{x})
$$

其中 $U_j$ 是待定的展开系数，也就是我们要解的未知量。

（ 或在 Deal.II 中这个东西被称为自由度，Degrees of Freedom, DoFs ）

## Shape Functions

有限元形函数（Shape Functions）是定义在每个单元（Element）上的局部函数，通常具有以下特点：
- 每个形函数在其对应的单元内非零，在其他单元上为零。
- 在节点处，形函数的值为1，而在其他节点处为0。
- 形函数可以是线性的或高阶多项式等。

**一个直观理解：** $\varphi_j(x)$ 是像帐篷一样的形函数，在某个点为1，周围线性下降到0。整个近似解 $u_h(x)$ 就是对这些帐篷函数赋予不同高度的叠加。

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

![alt text](image-1.png)

要确定 $u_h(x)$，我们只需要算出有限个系数（Coefficients）$U_j$ 。

</div>

## 建立线性系统

**记号：**$(f, g)= \int_{\Omega} f g \, d\mathbf{x}$

于是，离散解 $u_h$ 必须满足弱形式，但我们只需要它对有限个基函数 $\varphi_i$（作为测试函数）成立：
$$
(\nabla u_h, \nabla \varphi_i) = (\varphi_i, f), \quad i = 0, \dots, N-1
$$

将 $u_h$ 的表达式代入上式：
$$
\left(\nabla \left(\sum_{j=0}^{N-1} U_j \varphi_j\right), \nabla \varphi_i\right) = (\varphi_i, f)
$$
$$
\sum_{j=0}^{N-1} U_j (\nabla \varphi_j, \nabla \varphi_i) = (\varphi_i, f)
$$

## 建立线性系统

这可以写成一个线性代数方程组

$$
\mathbf{A} \mathbf{U} = \mathbf{F}
$$

刚度矩阵 (Stiffness Matrix) $A$:

$$
A_{ij} = (\nabla \varphi_i, \nabla \varphi_j) = \int_{\Omega} \nabla \varphi_i \cdot \nabla \varphi_j \, d\mathbf{x}
$$

右端向量 (Right Hand Side Vector) $F$:
$$
F_i = (\varphi_i, f) = \int_{\Omega} \varphi_i f \, d\mathbf{x}
$$

解向量 (Solution Vector) $U$:
$$
U = [U_0, U_1, \dots, U_{N-1}]^T
$$

## 求解线性系统

使用数值线性代数方法（如共轭梯度法、直接法等）求解线性系统 $\mathbf{A} \mathbf{U} = \mathbf{F}$，得到系数向量 $\mathbf{U}$。

一旦得到 $\mathbf{U}$，我们就可以构造出近似解 $u_h(\mathbf{x})$，并在整个域 $\Omega$ 上进行评估和可视化。

> 未尽事宜：
> - 给定一个不那么特殊的微分方程，如何给出其弱形式等以使用有限元方法求解？
> - 近似解 $u_h$ 是否接近精确解 $u$？
> - 达到一定精度所需的计算量是多少？ 是否最优？
>
> <p> </p>
>

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

![SORRY](nihao.png)

</div>

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

<hr/>

参考资料有：

[Wolfgang Bangerth's Video Lectures](https://www.math.colostate.edu/~bangerth/videos.html)

[The deal.II Library](https://dealii.org/current/doxygen/deal.II/)

如果您有任何意见或者建议，欢迎通过电子邮件联系我：

[yxyou@mail.bnu.edu.cn](mailto:yxyou@mail.bnu.edu.cn)

</div>
</div>

<p style="color: #ccccccff;">此页被有意地留作空白。</p>

## COPYRIGHT

**图片来源：**

除了在前页中标明的参考资料外，其他图片均来自于：

- [いらすとや](https://www.irasutoya.com/) ；
- [Wikipedia](https://www.wikipedia.org/) ；

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

**其他说明：**

- 本站幻灯片及相关样式/模板部分改编自 TonyCrane 的 [slide-template](https://github.com/TonyCrane/slide-template) ；

</div>

## Galerkin方法

**基本思想：**
Galerkin方法是一种用于求解偏微分方程和边界值问题的变分方法。它通过将待求解函数表示为一组试验函数的线性组合，然后将这些试验函数与待求解函数的乘积在整个求解域上积分，并令该积分为零，从而将原方程转化为一个求解线性方程组的问题。

摘自[北太天元科普: 从微分方程导出变分问题以及微分方程的弱解](https://www.bilibili.com/opus/927310798199455784)

[返航](#jlj)

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

> **应用步骤：**
> 
> - 根据问题的边界条件和物理特性，选择一组适当的试验函数。
> - 在整个求解域上建立积分方程，将原方程转化为积分形式的方程。
> - 通过求解线性代数方程组得到待求解函数的近似解。
> 
> <p> </p>
</div>

</div>

## 三角形的运算符们与分部积分公式

$$
\nabla = \left( \frac{\partial}{\partial x}, \frac{\partial}{\partial y} \right)
$$
$$
\Delta = \nabla \cdot \nabla = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}
$$

</div>

<div class="col">
散度定理的公式：

$$
\int_{\Omega} \nabla \cdot \mathbf{F} \, d\mathbf{x} = \int_{\partial \Omega} \mathbf{F} \cdot \mathbf{n} \, dS
$$
</div>

</div>

在高维空间中，分部积分公式的推广形式如下：

$$
\int_\Omega \nabla u \cdot \mathbf{v} , d\mathbf{x} = \int_{\partial \Omega} u , \mathbf{v} \cdot \mathbf{n} , dS - \int_\Omega u , \nabla \cdot \mathbf{v} , d\mathbf{x}
$$

针对我们的问题，则有：
$$
-\int_{\Omega} \varphi \Delta u \, d\mathbf{x} = \int_{\Omega} \nabla \varphi \cdot \nabla u \, d\mathbf{x} - \int_{\partial \Omega} \varphi \mathbf{n} \cdot \nabla u \, dS
$$