偏微分方程数值解法的3种核心方法：从有限差分到有限元
偏微分方程数值解法的3种核心方法：从有限差分到有限元
偏微分方程（PDE）是描述物理现象中连续变量随时间和空间变化的数学方程。由于PDE通常无法解析求解，因此需要使用数值方法来获得近似解。本文将介绍偏微分方程数值解法的三种核心方法：有限差分法、有限元法和有限体积法。
1. 偏微分方程数值解法的概述
偏微分方程（PDE）是描述物理现象中连续变量随时间和空间变化的数学方程。由于PDE通常无法解析求解，因此需要使用数值方法来获得近似解。
偏微分方程数值解法的基本原理是将PDE离散化为一个代数方程组，然后使用计算机求解该方程组。常用的数值方法包括有限差分法、有限元法和有限体积法。这些方法各有优缺点，适用于不同的PDE类型和应用领域。
在本章中，我们将介绍偏微分方程数值解法的基本概念、常用方法及其在实际应用中的示例。
2. 有限差分法
有限差分法是一种数值方法，用于求解偏微分方程（PDE）。它基于将偏微分方程离散化成一组代数方程，然后求解这些方程来获得近似解。
2.1 有限差分法的基本原理
2.1.1 时空离散化
有限差分法首先将偏微分方程在时空域上离散化。时空离散化是指将连续的时空域划分为离散的网格点，并在网格点上定义未知函数的值。
对于一个一维偏微分方程，时空离散化可以表示为：
$$
\frac{\partial u}{\partial t} \approx \frac{u_i^{n+1} - u_i^n}{\Delta t}
$$
$$
\frac{\partial u}{\partial x} \approx \frac{u_{i+1}^n - u_i^n}{\Delta x}
$$
其中，$u_i^n$ 表示网格点 $(i, n)$ 处的未知函数值，$\Delta t$ 和 $\Delta x$ 分别表示时间步长和空间步长。
2.1.2 差分格式的构造
离散化后，需要构造差分格式来近似偏微分方程的导数项。差分格式有多种类型，常用的包括：
$$
\frac{\partial u}{\partial x} \approx \frac{u_{i+1}^n - u_i^n}{\Delta x}
$$
$$
\frac{\partial u}{\partial x} \approx \frac{u_i^n - u_{i-1}^n}{\Delta x}
$$
$$
\frac{\partial u}{\partial x} \approx \frac{u_{i+1}^n - u_{i-1}^n}{2\Delta x}
$$
不同的差分格式具有不同的精度和稳定性，需要根据具体问题选择合适的差分格式。
2.2 有限差分法的求解方法
求解有限差分方程组有多种方法，包括：
2.2.1 直接求解法
直接求解法将有限差分方程组转化为一个线性方程组，然后使用高斯消元法或其他直接求解方法求解线性方程组。
2.2.2 迭代求解法
迭代求解法通过迭代的方式逐步逼近解。常用的迭代求解法包括：
$$
u_i^{n+1} = \frac{b_i - \sum_{j \neq i} a_{ij} u_j^n}{a_{ii}}
$$
$$
u_i^{n+1} = \frac{b_i - \sum_{j < i} a_{ij} u_j^{n+1} - \sum_{j > i} a_{ij} u_j^n}{a_{ii}}
$$
$$
u_i^{n+1} = (1 - \omega) u_i^n + \omega \frac{b_i - \sum_{j \neq i} a_{ij} u_j^{n+1}}{a_{ii}}
$$
其中，$\omega$ 为松弛因子，取值范围为 $(0, 2)$。
不同的迭代求解法具有不同的收敛速度和稳定性，需要根据具体问题选择合适的迭代求解法。
3. 有限元法
有限元法是一种将连续问题离散化为离散问题的数值方法。其基本原理是将求解域划分为有限个单元，并用有限个基函数来近似未知解。
3.1 有限元法的基本原理
3.1.1 弱形式的建立
弱形式是偏微分方程的一种积分形式，它可以将偏微分方程转化为一个等价的积分方程。对于二阶椭圆型偏微分方程，其弱形式为：
$$
\int_{\Omega} (\nabla u \cdot \nabla v - fu) d\Omega = 0
$$
其中，$u$ 是未知解，$v$ 是任意光滑函数，$\Omega$ 是求解域，$f$ 是已知函数。
3.1.2 有限元空间的构造
有限元空间是一个由基函数张成的函数空间。对于给定的求解域，可以采用不同的方法构造有限元空间，常用的方法有：
拉格朗日有限元：基函数为局部支持的连续函数，在单元内为多项式。
Hermite 有限元：基函数为局部支持的非连续函数，在单元内为多项式，并在单元边界上满足连续性条件。
3.2 有限元法的求解方法
有限元法求解偏微分方程的步骤如下：
离散化：将求解域划分为有限个单元，并构造有限元空间。
基函数选择：选择有限元空间中的基函数，并用它们来近似未知解。
Galerkin 方法：将弱形式代入基函数，并对基函数进行积分，得到离散方程组。
Petrov-Galerkin 方法：在 Galerkin 方法的基础上，引入权函数，得到加权离散方程组。
3.2.1 Galerkin 方法
Galerkin 方法是一种最常用的有限元求解方法。其基本思想是将弱形式代入基函数，并对基函数进行积分，得到离散方程组：
$$
\int_{\Omega} (\nabla u \cdot \nabla v - fu) d\Omega = 0, \quad \forall v \in V
$$
其中，$V$ 是有限元空间。
3.2.2 Petrov-Galerkin 方法
Petrov-Galerkin 方法是一种改进的 Galerkin 方法。其基本思想是在 Galerkin 方法的基础上，引入权函数 $w$，得到加权离散方程组：
$$
\int_{\Omega} (\nabla u \cdot \nabla w - fw) d\Omega = 0, \quad \forall w \in W
$$
其中，$W$ 是权函数空间。
4. 有限体积法
4.1 有限体积法的基本原理
4.1.1 控制体积的划分
有限体积法将计算域划分为一系列不重叠的控制体积，每个控制体积代表了计算域中的一小部分。控制体积的形状可以是任意多边形，但通常使用矩形或三角形。
4.1.2 守恒律的离散化
有限体积法的基本原理是基于守恒定律的离散化。对于一个控制体积，守恒定律表示控制体积内某一物理量的净流入率等于该物理量在控制体积内的变化率。
对于一个控制体积，守恒定律可以表示为：
$$
\frac{\partial}{\partial t} \int_{\Omega} u d\Omega + \int_{\partial \Omega} F(u) \cdot n dS = 0
$$
其中：
- $\Omega$ 是控制体积
- $u$ 是物理量
- $F(u)$ 是物理量的通量
- $n$ 是控制体积边界上的法向量
将守恒定律离散化后得到：
$$
\frac{\partial u_i}{\partial t} V_i + \sum_{j=1}^{N_f} F_{i,j} = 0
$$
其中：
- $u_i$ 是控制体积 $\Omega_i$ 中物理量的平均值
- $V_i$ 是控制体积 $\Omega_i$ 的体积
- $N_f$ 是控制体积 $\Omega_i$ 的面数
- $F_{i,j}$ 是控制体积 $\Omega_i$ 和 $\Omega_j$ 之间面的通量
4.2 有限体积法的求解方法
4.2.1 有限体积法的求解方法
有限体积法的求解方法通常分为以下几个步骤：
- 空间离散化：将计算域划分为一系列控制体积。
- 守恒律离散化：将守恒定律离散化为控制体积上的一组代数方程。
- 时间离散化：将代数方程组离散化为时间上的一个序列。
- 求解代数方程组：求解每个时间步长上的代数方程组。
4.2.2 数值通量的计算
数值通量是控制体积边界上物理量的近似值。数值通量的计算方法有很多种，常见的方法有：
- 中央差分格式：将通量近似为控制体积中心处的物理量。
- 迎风格式：将通量近似为上游控制体积的物理量。
- 混合格式：将通量近似为中央差分格式和迎风格式的加权平均。
5. 应用实例
偏微分方程数值解法在流体力学中有着广泛的应用，其中最具代表性的便是 Navier-Stokes 方程的求解。Navier-Stokes 方程是一组非线性偏微分方程，描述了粘性流体的运动规律。
5.1.1 Navier-Stokes 方程的求解
求解 Navier-Stokes 方程的数值方法有很多，其中有限差分法、有限元法和有限体积法是最常用的方法。
有限差分法
有限差分法将流体域离散成一个个网格单元，并利用差分格式将偏微分方程转化为代数方程组。这些代数方程组可以通过直接求解法或迭代求解法求解。
有限元法
有限元法将流体域离散成一系列有限元，并利用有限元空间构造弱形式方程。弱形式方程可以通过 Galerkin 方法或 Petrov-Galerkin 方法求解。
有限体积法
有限体积法将流体域离散成一个个控制体积，并利用守恒律对控制体积进行离散化。离散化后的守恒律方程可以通过有限体积法的求解方法求解。
5.1.2 流动可视化
偏微分方程数值解法还可以用于流动可视化。通过对流场数据的后处理，可以生成流线、等值线等可视化结果，直观地展示流体的运动规律。
例如，下图展示了使用有限体积法求解 Navier-Stokes 方程得到的二维流场可视化结果。