数值抛物型方程的变分形式与有限元方法
字数 4540 2025-11-06 12:40:40

数值抛物型方程的变分形式与有限元方法

好的，我们开始学习“数值抛物型方程的变分形式与有限元方法”。这是一个将有限元法应用于时间相关问题的核心主题。我们将从最基础的概念开始，逐步深入。

第一步：回顾抛物型方程及其经典形式

首先，我们明确研究对象。一个典型的线性抛物型方程是热传导方程：

\[\frac{\partial u}{\partial t} - \nabla \cdot (a(\mathbf{x}) \nabla u) = f(\mathbf{x}, t) \]

其中：

  • \(u = u(\mathbf{x}, t)\) 是未知函数（例如温度）。
  • \(a(\mathbf{x}) > 0\) 是扩散系数（例如热导率）。
  • \(f(\mathbf{x}, t)\) 是源项（例如热源）。
  • \(\frac{\partial u}{\partial t}\) 是时间导数项，体现了方程的“抛物”特性，即解随时间演化。
    这个方程需要初始条件 \(u(\mathbf{x}, 0) = u_0(\mathbf{x})\) 和边界条件（例如狄利克雷边界条件 \(u = g\) 在边界上）。

我们之前讨论的“有限差分法”是直接在偏微分方程上对导数进行离散。而“有限元法”则需要一个不同的出发点——变分形式（或弱形式）。

第二步：从经典形式到变分形式（空间离散）

变分形式是有限元法的基石。其核心思想是：我们不要求偏微分方程在空间每一点都严格成立，而是要求它在某种“平均”意义下成立。

  1. 乘以试验函数：我们将方程的两边同时乘以一个光滑的“试验函数” \(v(\mathbf{x})\)。这个函数在满足齐次边界条件的边界上为零（即如果我们的边界条件是 \(u=g\)，那么 \(v\) 在边界上满足 \(v=0\)）。

\[ \int_{\Omega} \left( \frac{\partial u}{\partial t} \right) v \, d\mathbf{x} - \int_{\Omega} \nabla \cdot (a \nabla u) v \, d\mathbf{x} = \int_{\Omega} f v \, d\mathbf{x} \]

  1. 应用格林公式（散度定理）：第二项包含二阶导数，为了降低对函数光滑性的要求（使得解可以存在于更广泛的函数空间中），我们对其应用格林公式。公式告诉我们：

\[ \int_{\Omega} \nabla \cdot (a \nabla u) v \, d\mathbf{x} = - \int_{\Omega} a \nabla u \cdot \nabla v \, d\mathbf{x} + \int_{\partial \Omega} (a \nabla u \cdot \mathbf{n}) v \, ds \]

其中 \(\partial \Omega\) 是区域边界，\(\mathbf{n}\) 是外法向量。由于我们的试验函数 \(v\) 在边界上为0，边界积分项 \(\int_{\partial \Omega} ... ds\) 消失了。

  1. 得到变分形式：将上述结果代回原方程，我们得到抛物型方程的变分形式（或弱形式）：

\[ \int_{\Omega} \frac{\partial u}{\partial t} v \, d\mathbf{x} + \int_{\Omega} a \nabla u \cdot \nabla v \, d\mathbf{x} = \int_{\Omega} f v \, d\mathbf{x} \]

这个形式的美妙之处在于，它只要求 \(u\) 和 \(v\) 具有一阶导数（而不是原来的二阶导数），放宽了解的存在性要求。

第三步：有限元离散——从连续到有限维

变分形式仍然是关于连续函数 \(u\) 和 \(v\) 的。有限元法的下一步是将其投影到一个有限的维度空间。

  1. 构造有限元空间：我们将求解区域 \(\Omega\) 划分成许多小的单元（如三角形、四边形），构成网格。然后在每个单元上定义简单的基础函数（最常见的是分段线性多项式，即“帐篷函数”）。所有这些基础函数 \(\{\phi_1(\mathbf{x}), \phi_2(\mathbf{x}), ..., \phi_N(\mathbf{x})\}\) 张成一个有限维的函数空间 \(V_h \subset V\)（\(V\) 是满足边界条件的连续函数空间）。

  2. 近似解：我们在有限元空间 \(V_h\) 中寻找近似解 \(u_h(\mathbf{x}, t)\)。它可以表示为基函数的线性组合：

\[ u_h(\mathbf{x}, t) = \sum_{j=1}^{N} c_j(t) \phi_j(\mathbf{x}) \]

注意，系数 \(c_j\) 现在是时间的函数，而基函数 \(\phi_j\) 只与空间有关。这是我们求解随时间演化问题的关键。

  1. 伽辽金方法：我们将近似解 \(u_h\) 代入变分形式，并要求变分形式对有限元空间中的所有试验函数 \(v_h \in V_h\) 都成立。由于基函数 \(\{\phi_i\}\) 张成了整个空间 \(V_h\)，我们只需要求它对每一个基函数 \(\phi_i\) 都成立即可。

将 \(u_h\) 和 \(v = \phi_i\) 代入变分形式：

\[\int_{\Omega} \frac{\partial}{\partial t} \left( \sum_{j=1}^{N} c_j(t) \phi_j(\mathbf{x}) \right) \phi_i \, d\mathbf{x} + \int_{\Omega} a \nabla \left( \sum_{j=1}^{N} c_j(t) \phi_j(\mathbf{x}) \right) \cdot \nabla \phi_i \, d\mathbf{x} = \int_{\Omega} f \phi_i \, d\mathbf{x} \]

交换求和与积分顺序，并利用 \(\frac{\partial}{\partial t}\) 只作用于 \(c_j(t)\)：

\[\sum_{j=1}^{N} \frac{dc_j(t)}{dt} \int_{\Omega} \phi_i \phi_j \, d\mathbf{x} + \sum_{j=1}^{N} c_j(t) \int_{\Omega} a \nabla \phi_i \cdot \nabla \phi_j \, d\mathbf{x} = \int_{\Omega} f \phi_i \, d\mathbf{x} \]

对于每一个 \(i = 1, 2, ..., N\)，这都构成一个方程。

第四步：得到常微分方程组系统

上一步的方程系统可以写成简洁的矩阵形式。我们定义：

  • 质量矩阵 \(M\)： 其元素为 \(M_{ij} = \int_{\Omega} \phi_i \phi_j \, d\mathbf{x}\)。
  • 刚度矩阵 \(K\)： 其元素为 \(K_{ij} = \int_{\Omega} a \nabla \phi_i \cdot \nabla \phi_j \, d\mathbf{x}\)。
  • 载荷向量 \(F(t)\)： 其元素为 \(F_i(t) = \int_{\Omega} f(\mathbf{x}, t) \phi_i \, d\mathbf{x}\)。
  • 解向量 \(C(t)\)： \(C(t) = [c_1(t), c_2(t), ..., c_N(t)]^T\)。

于是，我们得到了一个关于时间 \(t\) 的常微分方程组（ODEs）：

\[M \frac{dC(t)}{dt} + K C(t) = F(t) \]

这个系统被称为半离散系统，因为空间已经被有限元法离散了（体现在矩阵 \(M\) 和 \(K\) 上），但时间仍然是连续的。

第五步：时间积分方法

现在我们需要一个数值方法来求解这个常微分方程组。这被称为时间离散。

  1. θ-方法： 最常用的一类方法是θ-方法。我们将时间区间划分为小步长 \(\Delta t\)，并在时间点 \(t_n\) 和 \(t_{n+1}\) 之间建立近似。θ-方法将方程表示为：

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

其中 \(C^n\) 是 \(t_n\) 时刻的解向量近似。

  1. 常见格式：
  • 显式欧拉法（θ=0）： \((M / \Delta t) C^{n+1} = (M / \Delta t - K) C^n + F^n\)。需要求解 \(M\) 的逆，通常 \(M\) 被“集中”（对角化）以使计算简单，但时间步长 \(\Delta t\) 受稳定性限制很严（CFL条件）。
  • 隐式欧拉法（θ=1）： \((M / \Delta t + K) C^{n+1} = (M / \Delta t) C^n + F^{n+1}\)。这是无条件稳定的，但每一步都需要求解一个线性方程组。
    • 克兰克-尼科尔森法（θ=1/2）： 这是二阶精度且无条件稳定的方法，在精度和稳定性之间取得了很好的平衡，是最常用的选择之一。

\[ (M / \Delta t + \frac{K}{2}) C^{n+1} = (M / \Delta t - \frac{K}{2}) C^n + \frac{1}{2}(F^{n+1} + F^n) \]

第六步：总结与优势

总结一下，用有限元法求解抛物型方程的流程是：
建立PDE → 推导变分形式 → 空间有限元离散 → 得到ODEs系统 → 时间离散求解。

这种方法的优势在于：

  • 复杂几何： 能自然处理复杂的区域形状。
  • 自适应性： 可以方便地进行空间网格自适应（\(h\)-自适应）和多项式次数自适应（\(p\)-自适应）。
  • 边界条件： 自然边界条件（如诺伊曼条件）能自动融入变分形式中。
  • 理论基础扎实： 有成熟的数学理论进行误差分析和稳定性证明。

通过这六个步骤，我们完成了从抛物型方程的经典形式出发，最终得到一套完整的、可在计算机上实现的数值算法的全过程。

数值抛物型方程的变分形式与有限元方法 好的，我们开始学习“数值抛物型方程的变分形式与有限元方法”。这是一个将有限元法应用于时间相关问题的核心主题。我们将从最基础的概念开始，逐步深入。 第一步：回顾抛物型方程及其经典形式 首先，我们明确研究对象。一个典型的线性抛物型方程是热传导方程： \[ \frac{\partial u}{\partial t} - \nabla \cdot (a(\mathbf{x}) \nabla u) = f(\mathbf{x}, t) \] 其中： \( u = u(\mathbf{x}, t) \) 是未知函数（例如温度）。 \( a(\mathbf{x}) > 0 \) 是扩散系数（例如热导率）。 \( f(\mathbf{x}, t) \) 是源项（例如热源）。 \( \frac{\partial u}{\partial t} \) 是时间导数项，体现了方程的“抛物”特性，即解随时间演化。 这个方程需要初始条件 \( u(\mathbf{x}, 0) = u_ 0(\mathbf{x}) \) 和边界条件（例如狄利克雷边界条件 \( u = g \) 在边界上）。 我们之前讨论的“有限差分法”是直接在偏微分方程上对导数进行离散。而“有限元法”则需要一个不同的出发点——变分形式（或弱形式）。 第二步：从经典形式到变分形式（空间离散） 变分形式是有限元法的基石。其核心思想是：我们不要求偏微分方程在空间每一点都严格成立，而是要求它在某种“平均”意义下成立。 乘以试验函数 ：我们将方程的两边同时乘以一个光滑的“试验函数” \( v(\mathbf{x}) \)。这个函数在满足齐次边界条件的边界上为零（即如果我们的边界条件是 \( u=g \)，那么 \( v \) 在边界上满足 \( v=0 \)）。 \[ \int_ {\Omega} \left( \frac{\partial u}{\partial t} \right) v \, d\mathbf{x} - \int_ {\Omega} \nabla \cdot (a \nabla u) v \, d\mathbf{x} = \int_ {\Omega} f v \, d\mathbf{x} \] 应用格林公式（散度定理） ：第二项包含二阶导数，为了降低对函数光滑性的要求（使得解可以存在于更广泛的函数空间中），我们对其应用格林公式。公式告诉我们： \[ \int_ {\Omega} \nabla \cdot (a \nabla u) v \, d\mathbf{x} = - \int_ {\Omega} a \nabla u \cdot \nabla v \, d\mathbf{x} + \int_ {\partial \Omega} (a \nabla u \cdot \mathbf{n}) v \, ds \] 其中 \( \partial \Omega \) 是区域边界，\( \mathbf{n} \) 是外法向量。由于我们的试验函数 \( v \) 在边界上为0，边界积分项 \( \int_ {\partial \Omega} ... ds \) 消失了。 得到变分形式 ：将上述结果代回原方程，我们得到抛物型方程的 变分形式（或弱形式） ： \[ \int_ {\Omega} \frac{\partial u}{\partial t} v \, d\mathbf{x} + \int_ {\Omega} a \nabla u \cdot \nabla v \, d\mathbf{x} = \int_ {\Omega} f v \, d\mathbf{x} \] 这个形式的美妙之处在于，它只要求 \( u \) 和 \( v \) 具有一阶导数（而不是原来的二阶导数），放宽了解的存在性要求。 第三步：有限元离散——从连续到有限维 变分形式仍然是关于连续函数 \( u \) 和 \( v \) 的。有限元法的下一步是将其投影到一个有限的维度空间。 构造有限元空间 ：我们将求解区域 \( \Omega \) 划分成许多小的单元（如三角形、四边形），构成网格。然后在每个单元上定义简单的基础函数（最常见的是分段线性多项式，即“帐篷函数”）。所有这些基础函数 \( \{\phi_ 1(\mathbf{x}), \phi_ 2(\mathbf{x}), ..., \phi_ N(\mathbf{x})\} \) 张成一个有限维的函数空间 \( V_ h \subset V \)（\( V \) 是满足边界条件的连续函数空间）。 近似解 ：我们在有限元空间 \( V_ h \) 中寻找近似解 \( u_ h(\mathbf{x}, t) \)。它可以表示为基函数的线性组合： \[ u_ h(\mathbf{x}, t) = \sum_ {j=1}^{N} c_ j(t) \phi_ j(\mathbf{x}) \] 注意，系数 \( c_ j \) 现在是时间的函数，而基函数 \( \phi_ j \) 只与空间有关。这是我们求解随时间演化问题的关键。 伽辽金方法 ：我们将近似解 \( u_ h \) 代入变分形式，并要求变分形式对有限元空间中的所有试验函数 \( v_ h \in V_ h \) 都成立。由于基函数 \( \{\phi_ i\} \) 张成了整个空间 \( V_ h \)，我们只需要求它对每一个基函数 \( \phi_ i \) 都成立即可。 将 \( u_ h \) 和 \( v = \phi_ i \) 代入变分形式： \[ \int_ {\Omega} \frac{\partial}{\partial t} \left( \sum_ {j=1}^{N} c_ j(t) \phi_ j(\mathbf{x}) \right) \phi_ i \, d\mathbf{x} + \int_ {\Omega} a \nabla \left( \sum_ {j=1}^{N} c_ j(t) \phi_ j(\mathbf{x}) \right) \cdot \nabla \phi_ i \, d\mathbf{x} = \int_ {\Omega} f \phi_ i \, d\mathbf{x} \] 交换求和与积分顺序，并利用 \( \frac{\partial}{\partial t} \) 只作用于 \( c_ j(t) \)： \[ \sum_ {j=1}^{N} \frac{dc_ j(t)}{dt} \int_ {\Omega} \phi_ i \phi_ j \, d\mathbf{x} + \sum_ {j=1}^{N} c_ j(t) \int_ {\Omega} a \nabla \phi_ i \cdot \nabla \phi_ j \, d\mathbf{x} = \int_ {\Omega} f \phi_ i \, d\mathbf{x} \] 对于每一个 \( i = 1, 2, ..., N \)，这都构成一个方程。 第四步：得到常微分方程组系统 上一步的方程系统可以写成简洁的矩阵形式。我们定义： 质量矩阵 \( M \) ： 其元素为 \( M_ {ij} = \int_ {\Omega} \phi_ i \phi_ j \, d\mathbf{x} \)。 刚度矩阵 \( K \) ： 其元素为 \( K_ {ij} = \int_ {\Omega} a \nabla \phi_ i \cdot \nabla \phi_ j \, d\mathbf{x} \)。 载荷向量 \( F(t) \) ： 其元素为 \( F_ i(t) = \int_ {\Omega} f(\mathbf{x}, t) \phi_ i \, d\mathbf{x} \)。 解向量 \( C(t) \) ： \( C(t) = [ c_ 1(t), c_ 2(t), ..., c_ N(t) ]^T \)。 于是，我们得到了一个关于时间 \( t \) 的 常微分方程组（ODEs） ： \[ M \frac{dC(t)}{dt} + K C(t) = F(t) \] 这个系统被称为 半离散系统 ，因为空间已经被有限元法离散了（体现在矩阵 \( M \) 和 \( K \) 上），但时间仍然是连续的。 第五步：时间积分方法 现在我们需要一个数值方法来求解这个常微分方程组。这被称为 时间离散 。 θ-方法 ： 最常用的一类方法是θ-方法。我们将时间区间划分为小步长 \( \Delta t \)，并在时间点 \( t_ n \) 和 \( t_ {n+1} \) 之间建立近似。θ-方法将方程表示为： \[ M \frac{C^{n+1} - C^n}{\Delta t} + K [ \theta C^{n+1} + (1-\theta) C^n ] = \theta F^{n+1} + (1-\theta)F^n \] 其中 \( C^n \) 是 \( t_ n \) 时刻的解向量近似。 常见格式 ： 显式欧拉法（θ=0） ： \( (M / \Delta t) C^{n+1} = (M / \Delta t - K) C^n + F^n \)。需要求解 \( M \) 的逆，通常 \( M \) 被“集中”（对角化）以使计算简单，但时间步长 \( \Delta t \) 受稳定性限制很严（CFL条件）。 隐式欧拉法（θ=1） ： \( (M / \Delta t + K) C^{n+1} = (M / \Delta t) C^n + F^{n+1} \)。这是无条件稳定的，但每一步都需要求解一个线性方程组。 克兰克-尼科尔森法（θ=1/2） ： 这是二阶精度且无条件稳定的方法，在精度和稳定性之间取得了很好的平衡，是最常用的选择之一。 \[ (M / \Delta t + \frac{K}{2}) C^{n+1} = (M / \Delta t - \frac{K}{2}) C^n + \frac{1}{2}(F^{n+1} + F^n) \] 第六步：总结与优势 总结一下，用有限元法求解抛物型方程的流程是： 建立PDE → 推导变分形式 → 空间有限元离散 → 得到ODEs系统 → 时间离散求解。 这种方法的优势在于： 复杂几何 ： 能自然处理复杂的区域形状。 自适应性 ： 可以方便地进行空间网格自适应（\( h \)-自适应）和多项式次数自适应（\( p \)-自适应）。 边界条件 ： 自然边界条件（如诺伊曼条件）能自动融入变分形式中。 理论基础扎实 ： 有成熟的数学理论进行误差分析和稳定性证明。 通过这六个步骤，我们完成了从抛物型方程的经典形式出发，最终得到一套完整的、可在计算机上实现的数值算法的全过程。