抛物型偏微分方程的有限差分法

字数 4599 2025-12-09 21:31:08

抛物型偏微分方程的有限差分法

我们来学习抛物型偏微分方程的一种基础而重要的数值解法：有限差分法。我将从最简单的模型方程开始，逐步构建出完整的数值格式，并分析其核心性质。

第一步：模型方程与问题设定

我们考虑最简单的抛物型方程——一维热传导方程，作为模型：

\[\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}, \quad 0 < x < L, \quad t > 0 \]

其中 \(\alpha > 0\) 是常数（热扩散系数），\(u(x, t)\) 是未知函数（例如温度）。要唯一确定解，我们需要附加条件：

初始条件： \(u(x, 0) = f(x), \quad 0 \le x \le L\)。
边界条件：在 \(x=0\) 和 \(x=L\) 上给出。常见的有狄利克雷条件（给定函数值）或诺伊曼条件（给定导数值）。例如：
\(u(0, t) = g_0(t), \quad u(L, t) = g_1(t), \quad t > 0\)。

我们的目标是找到 \(u(x, t)\) 在定义域内各点的近似值。

第二步：区域的离散化（网格生成）

有限差分法的核心思想是用离散的网格点逼近连续的区域。我们在空间区间 \([0, L]\) 和时间区间 \([0, T]\) 上分别进行划分。

空间离散：将 \([0, L]\) 分为 \(M\) 等份，空间步长 \(\Delta x = h = L/M\)。离散点坐标为 \(x_i = ih, \quad i = 0, 1, ..., M\)。
时间离散：将 \([0, T]\) 分为 \(N\) 等份，时间步长 \(\Delta t = k = T/N\)。离散时刻为 \(t_n = nk, \quad n = 0, 1, ..., N\)。

这样，连续区域 \((x, t)\) 被一张矩形网格覆盖，网格节点为 \((x_i, t_n)\)。我们用 \(U_i^n\) 表示在节点 \((x_i, t_n)\) 处对精确解 \(u(x_i, t_n)\) 的数值近似。

第三步：微分算子的差分近似

有限差分法的关键是用函数在邻近节点值的代数组合（差分商）来近似导数。

对时间偏导 \(\frac{\partial u}{\partial t}\)：在点 \((x_i, t_n)\) 处，最常用的近似是向前差分：

\[ \frac{\partial u}{\partial t}(x_i, t_n) \approx \frac{U_i^{n+1} - U_i^n}{k} \]

其截断误差为 \(O(k)\)（一阶精度）。也可以使用更精确的格式，但这是最基础的。

对空间二阶偏导 \(\frac{\partial^2 u}{\partial x^2}\)：在点 \((x_i, t_n)\) 处，使用中心差分：

\[ \frac{\partial^2 u}{\partial x^2}(x_i, t_n) \approx \frac{U_{i-1}^n - 2U_i^n + U_{i+1}^n}{h^2} \]

其截断误差为 \(O(h^2)\)（二阶精度）。

第四步：构造显式差分格式（FTCS格式）

将上述差分近似代入原方程 \(\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}\)，一种最直接的方式是在时间层 \(t_n\) 上近似空间导数，得到：

\[\frac{U_i^{n+1} - U_i^n}{k} = \alpha \frac{U_{i-1}^n - 2U_i^n + U_{i+1}^n}{h^2} \]

整理后，得到著名的FTCS格式（Forward Time, Centered Space）：

\[U_i^{n+1} = U_i^n + r (U_{i-1}^n - 2U_i^n + U_{i+1}^n) \]

其中 \(r = \alpha k / h^2\) 是一个无量纲常数，称为网格比。

“显式”的含义：在这个格式中，新时间层 \(n+1\) 上某点 \(i\) 的值 \(U_i^{n+1}\)，只依赖于旧时间层 \(n\) 上相邻三点 \(U_{i-1}^n, U_i^n, U_{i+1}^n\) 的已知值。因此，如果我们已知第 \(n\) 时间层所有点的值，可以直接、独立地（无需解方程组）计算出第 \(n+1\) 时间层所有内部点 \((i=1,...,M-1)\) 的值。边界点 \((i=0, M)\) 的值由边界条件直接给出。

第五步：数值格式的稳定性分析（冯·诺依曼分析）

显式格式虽然简单，但有一个致命弱点：条件稳定。即时间步长 \(k\) 和空间步长 \(h\) 必须满足一定关系，否则计算中微小的舍入误差会指数级放大，使结果完全失真。

最常用的稳定性分析方法是冯·诺依曼（von Neumann）稳定性分析。其核心思想是考察单频波（傅里叶模式）的放大情况。

假设数值解具有形式 \(U_i^n = G^n e^{j\beta (ih)}\)，其中 \(j\) 是虚数单位，\(\beta\) 是波数，\(G\) 是放大因子。
将其代入 FTCS 格式的齐次部分（忽略边界条件）：

\[ G^{n+1} e^{j\beta ih} = G^n e^{j\beta ih} + r [G^n e^{j\beta (i-1)h} - 2G^n e^{j\beta ih} + G^n e^{j\beta (i+1)h}] \]

两边除以 \(G^n e^{j\beta ih}\)，得到放大因子：

\[ G = 1 + r(e^{-j\beta h} - 2 + e^{j\beta h}) = 1 - 4r \sin^2(\frac{\beta h}{2}) \]

稳定性要求：对于所有波数 \(\beta\)，误差模不增长，即 \(|G| \le 1\)。由此推出：

\[ |1 - 4r \sin^2(\frac{\beta h}{2})| \le 1 \quad \Rightarrow \quad -1 \le 1 - 4r \sin^2(\frac{\beta h}{2}) \le 1 \]

右边不等式自动成立。左边不等式要求 \(1 - 4r \sin^2(\frac{\beta h}{2}) \ge -1\)，即 \(4r \sin^2(\frac{\beta h}{2}) \le 2\)。对最坏情况 \(\sin^2(...)=1\)，得到稳定性条件：

\[ r = \frac{\alpha k}{h^2} \le \frac{1}{2} \]

这就是 FTCS 格式的稳定性限制。它意味着时间步长 \(k\) 必须小于等于 \(h^2/(2\alpha)\)。当需要高空间精度（\(h\) 很小）时，时间步长 \(k\) 必须取得非常小，计算代价激增。

第六步：构造隐式差分格式（BTCS格式）

为了克服显式格式的苛刻稳定性限制，我们发展出隐式格式。一个经典的方法是BTCS格式（Backward Time, Centered Space）。它在时间层 \(t_{n+1}\) 上近似空间导数：

\[\frac{U_i^{n+1} - U_i^n}{k} = \alpha \frac{U_{i-1}^{n+1} - 2U_i^{n+1} + U_{i+1}^{n+1}}{h^2} \]

整理后得到：

\[- r U_{i-1}^{n+1} + (1+2r)U_i^{n+1} - r U_{i+1}^{n+1} = U_i^n \]

“隐式”的含义：新时间层 \(n+1\) 上某点 \(i\) 的值 \(U_i^{n+1}\)，不仅与旧时间层 \(n\) 的值有关，还与同时间层 \(n+1\) 的邻近点 \(U_{i-1}^{n+1}, U_{i+1}^{n+1}\) 有关。因此，要得到第 \(n+1\) 时间层上所有内部点 \((i=1,...,M-1)\) 的值，必须联立求解一个线性方程组。这个方程组的系数矩阵是三对角矩阵，可以用高效、稳定的托马斯算法求解。
稳定性：对 BTCS 格式进行冯·诺依曼分析，可得其放大因子为 \(G = 1 / [1 + 4r \sin^2(\beta h/2)]\)。显然，对任意的网格比 \(r>0\)，都有 \(|G| < 1\)。因此 BTCS 格式是无条件稳定的。这意味着我们可以基于精度（而非稳定性）的要求相对独立地选择 \(k\) 和 \(h\)，通常可以取比显式格式大得多的时间步长。

第七步：更精确的隐式格式（Crank-Nicolson格式）

FTCS（一阶时间精度）和 BTCS（一阶时间精度）在时间上都是一阶精度。为了获得二阶时间精度，常用的格式是Crank-Nicolson格式。它在时间层 \(t_{n+1/2}\) 上对方程进行中心差分近似，等效于对空间导数在 \(t_n\) 和 \(t_{n+1}\) 两个时间层上取算术平均：

\[\frac{U_i^{n+1} - U_i^n}{k} = \frac{\alpha}{2} \left( \frac{U_{i-1}^n - 2U_i^n + U_{i+1}^n}{h^2} + \frac{U_{i-1}^{n+1} - 2U_i^{n+1} + U_{i+1}^{n+1}}{h^2} \right) \]

整理后：

\[- \frac{r}{2} U_{i-1}^{n+1} + (1+r)U_i^{n+1} - \frac{r}{2} U_{i+1}^{n+1} = \frac{r}{2} U_{i-1}^n + (1-r)U_i^n + \frac{r}{2} U_{i+1}^n \]

特点： Crank-Nicolson格式在时间和空间上都具有二阶精度 \(O(k^2 + h^2)\)，并且也是无条件稳定的。它结合了高精度和良好稳定性的优点，是求解抛物型方程最受欢迎的标准方法之一。同样，它也需要在每一个时间步求解一个三对角线性方程组。

总结：
有限差分法为求解抛物型方程提供了系统性的数值框架。我们从最简单的显式FTCS格式入手，理解了离散化、差分近似和稳定性分析的基本概念。为了克服显式格式的稳定性限制，我们引入了无条件稳定的隐式BTCS格式。最后，为了在时间和空间上都获得二阶精度，我们学习了经典的Crank-Nicolson格式。这些方法是数值求解抛物型偏微分方程的基石，可以进一步推广到更复杂的方程（如对流-扩散方程）和高维问题。

抛物型偏微分方程的有限差分法我们来学习抛物型偏微分方程的一种基础而重要的数值解法：有限差分法。我将从最简单的模型方程开始，逐步构建出完整的数值格式，并分析其核心性质。第一步：模型方程与问题设定我们考虑最简单的抛物型方程——一维热传导方程，作为模型： \[ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}, \quad 0 < x < L, \quad t > 0 \] 其中 \(\alpha > 0\) 是常数（热扩散系数），\(u(x, t)\) 是未知函数（例如温度）。要唯一确定解，我们需要附加条件：初始条件： \(u(x, 0) = f(x), \quad 0 \le x \le L\)。边界条件：在 \(x=0\) 和 \(x=L\) 上给出。常见的有狄利克雷条件（给定函数值）或诺伊曼条件（给定导数值）。例如： \(u(0, t) = g_ 0(t), \quad u(L, t) = g_ 1(t), \quad t > 0\)。我们的目标是找到 \(u(x, t)\) 在定义域内各点的近似值。第二步：区域的离散化（网格生成）有限差分法的核心思想是用离散的网格点逼近连续的区域。我们在空间区间 \([ 0, L]\) 和时间区间 \([ 0, T ]\) 上分别进行划分。空间离散：将 \([ 0, L]\) 分为 \(M\) 等份，空间步长 \(\Delta x = h = L/M\)。离散点坐标为 \(x_ i = ih, \quad i = 0, 1, ..., M\)。时间离散：将 \([ 0, T]\) 分为 \(N\) 等份，时间步长 \(\Delta t = k = T/N\)。离散时刻为 \(t_ n = nk, \quad n = 0, 1, ..., N\)。这样，连续区域 \((x, t)\) 被一张矩形网格覆盖，网格节点为 \((x_ i, t_ n)\)。我们用 \(U_ i^n\) 表示在节点 \((x_ i, t_ n)\) 处对精确解 \(u(x_ i, t_ n)\) 的数值近似。第三步：微分算子的差分近似有限差分法的关键是用函数在邻近节点值的代数组合（差分商）来近似导数。对时间偏导 \(\frac{\partial u}{\partial t}\) ：在点 \((x_ i, t_ n)\) 处，最常用的近似是向前差分： \[ \frac{\partial u}{\partial t}(x_ i, t_ n) \approx \frac{U_ i^{n+1} - U_ i^n}{k} \] 其截断误差为 \(O(k)\)（一阶精度）。也可以使用更精确的格式，但这是最基础的。对空间二阶偏导 \(\frac{\partial^2 u}{\partial x^2}\) ：在点 \((x_ i, t_ n)\) 处，使用中心差分： \[ \frac{\partial^2 u}{\partial x^2}(x_ i, t_ n) \approx \frac{U_ {i-1}^n - 2U_ i^n + U_ {i+1}^n}{h^2} \] 其截断误差为 \(O(h^2)\)（二阶精度）。第四步：构造显式差分格式（FTCS格式）将上述差分近似代入原方程 \(\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}\)，一种最直接的方式是在时间层 \(t_ n\) 上近似空间导数，得到： \[ \frac{U_ i^{n+1} - U_ i^n}{k} = \alpha \frac{U_ {i-1}^n - 2U_ i^n + U_ {i+1}^n}{h^2} \] 整理后，得到著名的 FTCS格式（Forward Time, Centered Space）： \[ U_ i^{n+1} = U_ i^n + r (U_ {i-1}^n - 2U_ i^n + U_ {i+1}^n) \] 其中 \(r = \alpha k / h^2\) 是一个无量纲常数，称为网格比。 “显式”的含义：在这个格式中，新时间层 \(n+1\) 上某点 \(i\) 的值 \(U_ i^{n+1}\)，只依赖于旧时间层 \(n\) 上相邻三点 \(U_ {i-1}^n, U_ i^n, U_ {i+1}^n\) 的已知值。因此，如果我们已知第 \(n\) 时间层所有点的值，可以直接、独立地（无需解方程组）计算出第 \(n+1\) 时间层所有内部点 \((i=1,...,M-1)\) 的值。边界点 \((i=0, M)\) 的值由边界条件直接给出。第五步：数值格式的稳定性分析（冯·诺依曼分析）显式格式虽然简单，但有一个致命弱点：条件稳定。即时间步长 \(k\) 和空间步长 \(h\) 必须满足一定关系，否则计算中微小的舍入误差会指数级放大，使结果完全失真。最常用的稳定性分析方法是冯·诺依曼（von Neumann）稳定性分析。其核心思想是考察单频波（傅里叶模式）的放大情况。假设数值解具有形式 \(U_ i^n = G^n e^{j\beta (ih)}\)，其中 \(j\) 是虚数单位，\(\beta\) 是波数，\(G\) 是放大因子。将其代入 FTCS 格式的齐次部分（忽略边界条件）： \[ G^{n+1} e^{j\beta ih} = G^n e^{j\beta ih} + r [ G^n e^{j\beta (i-1)h} - 2G^n e^{j\beta ih} + G^n e^{j\beta (i+1)h} ] \] 两边除以 \(G^n e^{j\beta ih}\)，得到放大因子： \[ G = 1 + r(e^{-j\beta h} - 2 + e^{j\beta h}) = 1 - 4r \sin^2(\frac{\beta h}{2}) \] 稳定性要求：对于所有波数 \(\beta\)，误差模不增长，即 \(|G| \le 1\)。由此推出： \[ |1 - 4r \sin^2(\frac{\beta h}{2})| \le 1 \quad \Rightarrow \quad -1 \le 1 - 4r \sin^2(\frac{\beta h}{2}) \le 1 \] 右边不等式自动成立。左边不等式要求 \(1 - 4r \sin^2(\frac{\beta h}{2}) \ge -1\)，即 \(4r \sin^2(\frac{\beta h}{2}) \le 2\)。对最坏情况 \(\sin^2(...)=1\)，得到稳定性条件： \[ r = \frac{\alpha k}{h^2} \le \frac{1}{2} \] 这就是 FTCS 格式的稳定性限制。它意味着时间步长 \(k\) 必须小于等于 \(h^2/(2\alpha)\)。当需要高空间精度（\(h\) 很小）时，时间步长 \(k\) 必须取得非常小，计算代价激增。第六步：构造隐式差分格式（BTCS格式）为了克服显式格式的苛刻稳定性限制，我们发展出隐式格式。一个经典的方法是 BTCS格式（Backward Time, Centered Space）。它在时间层 \(t_ {n+1}\) 上近似空间导数： \[ \frac{U_ i^{n+1} - U_ i^n}{k} = \alpha \frac{U_ {i-1}^{n+1} - 2U_ i^{n+1} + U_ {i+1}^{n+1}}{h^2} \] 整理后得到： \[ r U_ {i-1}^{n+1} + (1+2r)U_ i^{n+1} - r U_ {i+1}^{n+1} = U_ i^n \] “隐式”的含义：新时间层 \(n+1\) 上某点 \(i\) 的值 \(U_ i^{n+1}\)，不仅与旧时间层 \(n\) 的值有关，还与同时间层 \(n+1\) 的邻近点 \(U_ {i-1}^{n+1}, U_ {i+1}^{n+1}\) 有关。因此，要得到第 \(n+1\) 时间层上所有内部点 \((i=1,...,M-1)\) 的值，必须联立求解一个线性方程组。这个方程组的系数矩阵是三对角矩阵，可以用高效、稳定的托马斯算法求解。稳定性：对 BTCS 格式进行冯·诺依曼分析，可得其放大因子为 \(G = 1 / [ 1 + 4r \sin^2(\beta h/2)]\)。显然，对任意的网格比 \(r>0\)，都有 \(|G| < 1\)。因此 BTCS 格式是无条件稳定的。这意味着我们可以基于精度（而非稳定性）的要求相对独立地选择 \(k\) 和 \(h\)，通常可以取比显式格式大得多的时间步长。第七步：更精确的隐式格式（Crank-Nicolson格式） FTCS（一阶时间精度）和 BTCS（一阶时间精度）在时间上都是一阶精度。为了获得二阶时间精度，常用的格式是 Crank-Nicolson格式。它在时间层 \(t_ {n+1/2}\) 上对方程进行中心差分近似，等效于对空间导数在 \(t_ n\) 和 \(t_ {n+1}\) 两个时间层上取算术平均： \[ \frac{U_ i^{n+1} - U_ i^n}{k} = \frac{\alpha}{2} \left( \frac{U_ {i-1}^n - 2U_ i^n + U_ {i+1}^n}{h^2} + \frac{U_ {i-1}^{n+1} - 2U_ i^{n+1} + U_ {i+1}^{n+1}}{h^2} \right) \] 整理后： \[ \frac{r}{2} U_ {i-1}^{n+1} + (1+r)U_ i^{n+1} - \frac{r}{2} U_ {i+1}^{n+1} = \frac{r}{2} U_ {i-1}^n + (1-r)U_ i^n + \frac{r}{2} U_ {i+1}^n \] 特点： Crank-Nicolson格式在时间和空间上都具有二阶精度 \(O(k^2 + h^2)\)，并且也是无条件稳定的。它结合了高精度和良好稳定性的优点，是求解抛物型方程最受欢迎的标准方法之一。同样，它也需要在每一个时间步求解一个三对角线性方程组。总结：有限差分法为求解抛物型方程提供了系统性的数值框架。我们从最简单的显式FTCS格式入手，理解了离散化、差分近似和稳定性分析的基本概念。为了克服显式格式的稳定性限制，我们引入了无条件稳定的隐式BTCS格式。最后，为了在时间和空间上都获得二阶精度，我们学习了经典的Crank-Nicolson格式。这些方法是数值求解抛物型偏微分方程的基石，可以进一步推广到更复杂的方程（如对流-扩散方程）和高维问题。