抛物型偏微分方程的有限元方法基础（Finite Element Method for Parabolic Partial Differential Equations）
字数 4398 2025-12-21 09:36:59

抛物型偏微分方程的有限元方法基础（Finite Element Method for Parabolic Partial Differential Equations）

我们来系统学习抛物型偏微分方程（如热传导方程）的有限元方法。这是一个从连续方程到离散数值求解的关键过程。

第一步：从一个经典模型问题出发

我们考虑一个最简单的模型——带有齐次狄利克雷边界条件的一维热传导方程：

  1. 控制方程：在空间-时间区域 \(\Omega \times (0, T] = (0, L) \times (0, T]\) 上

\[ \frac{\partial u}{\partial t} - \kappa \frac{\partial^2 u}{\partial x^2} = f(x, t)， \quad \kappa > 0 \]

这里 \(u(x, t)\) 是未知函数（如温度），\(\kappa\) 是扩散系数，\(f\) 是源项。

  1. 初始条件：\(u(x, 0) = u_0(x)\)，对所有 \(x \in (0, L)\)。

  2. 边界条件：\(u(0, t) = u(L, t) = 0\)，对所有 \(t \in (0, T]\)。

这是典型的抛物型方程，其解随时间演化，具有“平滑效应”。

第二步：从强形式到弱形式（空间离散准备）

有限元方法不从原始的偏微分方程（强形式）直接离散，而是先将其转化为弱形式（或变分形式）。这降低了导数阶数的要求，是处理复杂边界和区域的基础。

  1. 引入试探函数空间：我们寻找的解 \(u\) 在任意时刻 \(t\) 都应属于一个函数空间 \(V\)。由于边界条件是 \(u=0\)，我们定义空间：

\[ V = H^1_0(0, L) = \{ v \in H^1(0, L) \mid v(0)=v(L)=0 \} \]

其中 \(H^1\) 是索伯列夫空间，包含了函数及其一阶导数都平方可积的函数。

  1. 与一个检验函数作内积：在强形式方程两边同时乘以一个任意的检验函数 \(v(x) \in V\)，并在空间域 \([0, L]\) 上积分：

\[ \int_0^L \left( \frac{\partial u}{\partial t} v - \kappa \frac{\partial^2 u}{\partial x^2} v \right) dx = \int_0^L f v dx \]

  1. 应用分部积分：对二阶导数项进行分部积分（这是关键一步），利用 \(v\) 在边界上也为零的条件，消除边界项：

\[ - \int_0^L \kappa \frac{\partial^2 u}{\partial x^2} v dx = \int_0^L \kappa \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} dx \]

  1. 得到半离散弱形式：将上式代回，得到关于空间的弱形式：对于每个固定的 \(t > 0\)，寻找 \(u(\cdot, t) \in V\)，使得对所有的 \(v \in V\)，有

\[ \int_0^L \frac{\partial u}{\partial t} v dx + \int_0^L \kappa \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} dx = \int_0^L f v dx \]

这个形式只包含一阶导数，对解的正则性要求更低，是有限元离散的出发点。

第三步：空间离散——有限元逼近

现在，我们将无限维的函数空间 \(V\) 近似为一个有限维的子空间 \(V_h \subset V\)。

  1. 构造有限元空间 \(V_h\)：
  • 将区间 \([0, L]\) 划分为 \(N\) 个子区间（单元），节点为 \(x_0=0, x_1, \dots, x_N=L\)。
  • 我们选择最简单的分段线性有限元空间。\(V_h\) 由所有在网格节点上连续、在每个单元上是线性多项式、且在边界 \(x=0, L\) 处取值为零的函数构成。
  • 这个空间有一组基函数 \(\{\phi_j(x)\}_{j=1}^{N-1}\)，其中 \(\phi_j(x)\) 是“帽函数”：在节点 \(x_j\) 处值为1，在其他节点处值为0，且在相邻单元上线性变化。
  1. 在半离散弱形式中进行伽辽金逼近：
  • 寻找近似解 \(u_h(x, t) \in V_h\)，使其满足弱形式：对所有 \(v_h \in V_h\)，有

\[ \int_0^L \frac{\partial u_h}{\partial t} v_h dx + \int_0^L \kappa \frac{\partial u_h}{\partial x} \frac{\partial v_h}{\partial x} dx = \int_0^L f v_h dx \]

  • 将 \(u_h\) 用基函数展开：\(u_h(x, t) = \sum_{j=1}^{N-1} U_j(t) \phi_j(x)\)。这里 \(U_j(t)\) 是随时间变化的节点系数（即我们要求的未知量）。
  • 依次取检验函数 \(v_h = \phi_i(x) \ (i=1, \dots, N-1)\)，代入上式。由于方程对所有基函数都成立，我们就得到了一个关于系数 \(U_j(t)\) 的常微分方程组。

第四步：得到半离散系统（矩阵形式）

将展开式代入弱形式，整理后得到：

  1. 质量矩阵（Mass Matrix） \(M\)：元素为 \(M_{ij} = \int_0^L \phi_i(x) \phi_j(x) dx\)。它来自于时间导数项：

\[ \int_0^L \frac{\partial u_h}{\partial t} \phi_i dx = \sum_j \frac{dU_j}{dt} \int_0^L \phi_j \phi_i dx = \sum_j M_{ij} \frac{dU_j}{dt} \]

  1. 刚度矩阵（Stiffness Matrix） \(K\)：元素为 \(K_{ij} = \int_0^L \kappa \frac{d\phi_i}{dx} \frac{d\phi_j}{dx} dx\)。它来自于扩散项：

\[ \int_0^L \kappa \frac{\partial u_h}{\partial x} \frac{d\phi_i}{dx} dx = \sum_j U_j \int_0^L \kappa \frac{d\phi_j}{dx} \frac{d\phi_i}{dx} dx = \sum_j K_{ij} U_j \]

  1. 载荷向量（Load Vector） \(F(t)\)：元素为 \(F_i(t) = \int_0^L f(x, t) \phi_i(x) dx\)。

  2. 半离散常微分方程组：

\[ M \frac{dU(t)}{dt} + K U(t) = F(t) \]

其中 \(U(t) = [U_1(t), \dots, U_{N-1}(t)]^T\)。这是一个关于时间 \(t\) 的线性常微分方程组，初始条件 \(U(0)\) 由初始函数 \(u_0(x)\) 在节点上的值（或其投影）给出。

第五步：时间离散——差分格式

现在，我们需要用数值方法求解上面的常微分方程组。常用的方法是有限差分法。

  1. 将时间区间 \([0, T]\) 离散：设时间步长为 \(\Delta t\)，时间层 \(t_n = n\Delta t\)。记 \(U^n \approx U(t_n)\)， \(F^n = F(t_n)\)。

  2. θ-方法：这是一族单步格式，在 \(t_n\) 和 \(t_{n+1}\) 之间对方程进行近似：

\[ M \frac{U^{n+1} - U^n}{\Delta t} + K [\theta U^{n+1} + (1-\theta) U^n] = \theta F^{n+1} + (1-\theta) F^n \]

其中 \(\theta \in [0, 1]\) 是一个参数。

  • \(\theta = 0\)：显式欧拉法。格式简单，但通常要求时间步长非常小（\(\Delta t = O(h^2)\)， \(h\) 为空间步长）才能稳定。
  • \(\theta = 1/2\)：Crank-Nicolson 格式。具有二阶时间精度，且是无条件稳定的（对抛物问题），是常用的选择。
  • \(\theta = 1\)：隐式欧拉法。仅有一阶时间精度，但无条件稳定，且具有强烈的数值耗散性，能抑制振荡。
  1. 得到全离散代数系统：将上述方程整理，得到在每一时间步需求解的线性方程组：

\[ (M + \theta \Delta t K) U^{n+1} = (M - (1-\theta) \Delta t K) U^n + \Delta t [\theta F^{n+1} + (1-\theta) F^n] \]

当 \(\theta > 0\) 时，系统矩阵 \((M + \theta \Delta t K)\) 是对称正定、稀疏的（因为 \(M\) 和 \(K\) 都是），可以用高效的稀疏矩阵求解器（如共轭梯度法）求解。

第六步：方法特性总结

  1. 优点：

    • 几何灵活性：有限元天然适应复杂区域和不规则网格。
    • 边界条件处理：自然边界条件（诺伊曼型）在弱形式中自动包含，狄利克雷条件易于施加。
    • 理论基础扎实：有成熟的误差分析和收敛性理论。
    • 高精度潜力：通过提高多项式阶数（p-型）或加密网格（h-型）可系统提高精度。

  2. 核心概念：

    • 弱形式：降低导数要求，是有限元的起点。
    • 有限维子空间：用简单多项式（在单元上）的组合逼近复杂解。
    • 基函数：通常具有局部紧支集，导致稀疏矩阵。
    • 质量矩阵与刚度矩阵：分别体现解的“惯性”和“刚度”，是离散系统的核心。
    • 空间半离散 + 时间全离散：这是求解发展方程的典型策略。

以上就是抛物型方程有限元方法的核心骨架。从连续的物理方程出发，通过空间弱形式化、有限维逼近得到半离散系统，再通过时间差分得到最终可计算的代数系统，每一步都体现了将连续问题转化为离散计算问题的数学思想。

抛物型偏微分方程的有限元方法基础（Finite Element Method for Parabolic Partial Differential Equations） 我们来系统学习抛物型偏微分方程（如热传导方程）的有限元方法。这是一个从连续方程到离散数值求解的关键过程。 第一步：从一个经典模型问题出发 我们考虑一个最简单的模型——带有齐次狄利克雷边界条件的一维热传导方程： 控制方程 ：在空间-时间区域 \( \Omega \times (0, T] = (0, L) \times (0, T ] \) 上 \[ \frac{\partial u}{\partial t} - \kappa \frac{\partial^2 u}{\partial x^2} = f(x, t)， \quad \kappa > 0 \] 这里 \(u(x, t)\) 是未知函数（如温度），\(\kappa\) 是扩散系数，\(f\) 是源项。 初始条件 ：\(u(x, 0) = u_ 0(x)\)，对所有 \(x \in (0, L)\)。 边界条件 ：\(u(0, t) = u(L, t) = 0\)，对所有 \(t \in (0, T ]\)。 这是典型的抛物型方程，其解随时间演化，具有“平滑效应”。 第二步：从强形式到弱形式（空间离散准备） 有限元方法不从原始的偏微分方程（ 强形式 ）直接离散，而是先将其转化为 弱形式 （或变分形式）。这降低了导数阶数的要求，是处理复杂边界和区域的基础。 引入试探函数空间 ：我们寻找的解 \(u\) 在任意时刻 \(t\) 都应属于一个函数空间 \(V\)。由于边界条件是 \(u=0\)，我们定义空间： \[ V = H^1_ 0(0, L) = \{ v \in H^1(0, L) \mid v(0)=v(L)=0 \} \] 其中 \(H^1\) 是索伯列夫空间，包含了函数及其一阶导数都平方可积的函数。 与一个检验函数作内积 ：在强形式方程两边同时乘以一个任意的 检验函数 \(v(x) \in V\)，并在空间域 \([ 0, L ]\) 上积分： \[ \int_ 0^L \left( \frac{\partial u}{\partial t} v - \kappa \frac{\partial^2 u}{\partial x^2} v \right) dx = \int_ 0^L f v dx \] 应用分部积分 ：对二阶导数项进行分部积分（这是关键一步），利用 \(v\) 在边界上也为零的条件，消除边界项： \[ \int_ 0^L \kappa \frac{\partial^2 u}{\partial x^2} v dx = \int_ 0^L \kappa \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} dx \] 得到半离散弱形式 ：将上式代回，得到 关于空间的弱形式 ：对于每个固定的 \(t > 0\)，寻找 \(u(\cdot, t) \in V\)，使得对所有的 \(v \in V\)，有 \[ \int_ 0^L \frac{\partial u}{\partial t} v dx + \int_ 0^L \kappa \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} dx = \int_ 0^L f v dx \] 这个形式只包含一阶导数，对解的正则性要求更低，是有限元离散的出发点。 第三步：空间离散——有限元逼近 现在，我们将无限维的函数空间 \(V\) 近似为一个有限维的子空间 \(V_ h \subset V\)。 构造有限元空间 \(V_ h\)： 将区间 \([ 0, L]\) 划分为 \(N\) 个子区间（单元），节点为 \(x_ 0=0, x_ 1, \dots, x_ N=L\)。 我们选择最简单的 分段线性 有限元空间。\(V_ h\) 由所有在网格节点上连续、在每个单元上是线性多项式、且在边界 \(x=0, L\) 处取值为零的函数构成。 这个空间有一组 基函数 \(\{\phi_ j(x)\}_ {j=1}^{N-1}\)，其中 \(\phi_ j(x)\) 是“帽函数”：在节点 \(x_ j\) 处值为1，在其他节点处值为0，且在相邻单元上线性变化。 在半离散弱形式中进行伽辽金逼近 ： 寻找近似解 \(u_ h(x, t) \in V_ h\)，使其满足弱形式：对所有 \(v_ h \in V_ h\)，有 \[ \int_ 0^L \frac{\partial u_ h}{\partial t} v_ h dx + \int_ 0^L \kappa \frac{\partial u_ h}{\partial x} \frac{\partial v_ h}{\partial x} dx = \int_ 0^L f v_ h dx \] 将 \(u_ h\) 用基函数展开：\(u_ h(x, t) = \sum_ {j=1}^{N-1} U_ j(t) \phi_ j(x)\)。这里 \(U_ j(t)\) 是 随时间变化的节点系数 （即我们要求的未知量）。 依次取检验函数 \(v_ h = \phi_ i(x) \ (i=1, \dots, N-1)\)，代入上式。由于方程对所有基函数都成立，我们就得到了一个关于系数 \(U_ j(t)\) 的 常微分方程组 。 第四步：得到半离散系统（矩阵形式） 将展开式代入弱形式，整理后得到： 质量矩阵（Mass Matrix） \(M\)：元素为 \(M_ {ij} = \int_ 0^L \phi_ i(x) \phi_ j(x) dx\)。它来自于时间导数项： \[ \int_ 0^L \frac{\partial u_ h}{\partial t} \phi_ i dx = \sum_ j \frac{dU_ j}{dt} \int_ 0^L \phi_ j \phi_ i dx = \sum_ j M_ {ij} \frac{dU_ j}{dt} \] 刚度矩阵（Stiffness Matrix） \(K\)：元素为 \(K_ {ij} = \int_ 0^L \kappa \frac{d\phi_ i}{dx} \frac{d\phi_ j}{dx} dx\)。它来自于扩散项： \[ \int_ 0^L \kappa \frac{\partial u_ h}{\partial x} \frac{d\phi_ i}{dx} dx = \sum_ j U_ j \int_ 0^L \kappa \frac{d\phi_ j}{dx} \frac{d\phi_ i}{dx} dx = \sum_ j K_ {ij} U_ j \] 载荷向量（Load Vector） \(F(t)\)：元素为 \(F_ i(t) = \int_ 0^L f(x, t) \phi_ i(x) dx\)。 半离散常微分方程组 ： \[ M \frac{dU(t)}{dt} + K U(t) = F(t) \] 其中 \(U(t) = [ U_ 1(t), \dots, U_ {N-1}(t)]^T\)。这是一个关于时间 \(t\) 的线性常微分方程组，初始条件 \(U(0)\) 由初始函数 \(u_ 0(x)\) 在节点上的值（或其投影）给出。 第五步：时间离散——差分格式 现在，我们需要用数值方法求解上面的常微分方程组。常用的方法是有限差分法。 将时间区间 \([ 0, T]\) 离散 ：设时间步长为 \(\Delta t\)，时间层 \(t_ n = n\Delta t\)。记 \(U^n \approx U(t_ n)\)， \(F^n = F(t_ n)\)。 θ-方法 ：这是一族单步格式，在 \(t_ n\) 和 \(t_ {n+1}\) 之间对方程进行近似： \[ M \frac{U^{n+1} - U^n}{\Delta t} + K [ \theta U^{n+1} + (1-\theta) U^n ] = \theta F^{n+1} + (1-\theta) F^n \] 其中 \(\theta \in [ 0, 1 ]\) 是一个参数。 \(\theta = 0\)： 显式欧拉法 。格式简单，但通常要求时间步长非常小（\(\Delta t = O(h^2)\)， \(h\) 为空间步长）才能稳定。 \(\theta = 1/2\)： Crank-Nicolson 格式 。具有二阶时间精度，且是无条件稳定的（对抛物问题），是常用的选择。 \(\theta = 1\)： 隐式欧拉法 。仅有一阶时间精度，但无条件稳定，且具有强烈的数值耗散性，能抑制振荡。 得到全离散代数系统 ：将上述方程整理，得到在每一时间步需求解的线性方程组： \[ (M + \theta \Delta t K) U^{n+1} = (M - (1-\theta) \Delta t K) U^n + \Delta t [ \theta F^{n+1} + (1-\theta) F^n ] \] 当 \(\theta > 0\) 时，系统矩阵 \((M + \theta \Delta t K)\) 是 对称正定、稀疏 的（因为 \(M\) 和 \(K\) 都是），可以用高效的稀疏矩阵求解器（如共轭梯度法）求解。 第六步：方法特性总结 优点 ： 几何灵活性 ：有限元天然适应复杂区域和不规则网格。 边界条件处理 ：自然边界条件（诺伊曼型）在弱形式中自动包含，狄利克雷条件易于施加。 理论基础扎实 ：有成熟的误差分析和收敛性理论。 高精度潜力 ：通过提高多项式阶数（p-型）或加密网格（h-型）可系统提高精度。 核心概念 ： 弱形式 ：降低导数要求，是有限元的起点。 有限维子空间 ：用简单多项式（在单元上）的组合逼近复杂解。 基函数 ：通常具有局部紧支集，导致稀疏矩阵。 质量矩阵与刚度矩阵 ：分别体现解的“惯性”和“刚度”，是离散系统的核心。 空间半离散 + 时间全离散 ：这是求解发展方程的典型策略。 以上就是抛物型方程有限元方法的核心骨架。从连续的物理方程出发，通过 空间弱形式化 、 有限维逼近 得到半离散系统，再通过 时间差分 得到最终可计算的代数系统，每一步都体现了将连续问题转化为离散计算问题的数学思想。