数值抛物型方程的有限元法

字数 4225 2025-11-08 10:03:13

数值抛物型方程的有限元法

好的，我们开始学习“数值抛物型方程的有限元法”。这个方法结合了处理空间导数的有限元离散和处理时间导数的差分离散，是求解与时间相关的扩散、热传导等问题的强大工具。

第一步：理解问题原型

我们考虑一个典型的线性抛物型方程——一维热传导方程，并配以初边值条件：

控制方程: \(\frac{\partial u}{\partial t} - \alpha \frac{\partial^2 u}{\partial x^2} = f(x, t), \quad x \in (0, L), \quad t \in (0, T]\)
这里，\(u(x, t)\) 是未知函数（例如温度），\(\alpha > 0\) 是扩散系数，\(f(x, t)\) 是源项。
初始条件: \(u(x, 0) = u_0(x), \quad x \in [0, L]\)
这给出了整个区域在初始时刻 \(t=0\) 的状态。
边界条件: 例如狄利克雷边界条件 \(u(0, t) = g_0(t), \quad u(L, t) = g_L(t)\)
- 这规定了区域边界在所有时刻的状态。

我们的目标是在整个时空区域 \([0, L] \times [0, T]\) 上找到满足上述条件的函数 \(u(x, t)\) 的近似解。

第二步：空间离散的弱形式与有限元逼近

有限元法的核心第一步是将偏微分方程转化为积分形式的“弱形式”，从而降低对解的光滑性要求。

构造弱形式：

将控制方程两边同时乘以一个“试探函数” \(v(x)\)，该函数在满足齐次边界条件的函数空间中选取（例如，若边界条件是 \(u=0\)，则要求 \(v(0)=v(L)=0\)）。
在空间域 \([0, L]\) 上积分： \(\int_0^L \frac{\partial u}{\partial t} v \, dx - \alpha \int_0^L \frac{\partial^2 u}{\partial x^2} v \, dx = \int_0^L f v \, dx\)。
对含二阶导数的项应用分部积分（这是关键一步）： \(-\alpha \int_0^L \frac{\partial^2 u}{\partial x^2} v \, dx = \alpha \int_0^L \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} \, dx - \alpha \left[ \frac{\partial u}{\partial x} v \right]_0^L\)。
如果边界条件使得边界项消失（例如 \(v=0\) 在边界上），我们得到弱形式： \(\int_0^L \frac{\partial u}{\partial t} v \, dx + \alpha \int_0^L \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} \, dx = \int_0^L f v \, dx\)。
这个形式只要求解 \(u\) 具有一阶导数，比原始方程要求二阶导数更“弱”。

有限元离散空间：

我们将区间 \([0, L]\) 分割成许多小单元（子区间）。
- 在每个单元上，我们定义简单的多项式函数（最常见的是线性函数），这些函数在单元内部光滑，在单元连接处连续。这些多项式构成了分片多项式函数空间。
我们选择一组基函数 \(\{\phi_1(x), \phi_2(x), \dots, \phi_N(x)\}\) 来张成这个空间。每个基函数 \(\phi_i(x)\) 通常只在包含第 \(i\) 个节点的少数几个单元上非零。
我们用这组基函数的线性组合来近似真实解： \(u(x, t) \approx u_h(x, t) = \sum_{j=1}^N u_j(t) \phi_j(x)\)。
这里的关键是，我们将解的时空变量分离了：空间变化由基函数 \(\phi_j(x)\) 描述，而时间变化由待求的系数 \(u_j(t)\) 描述。

伽辽金方法：

将近似解 \(u_h(x, t)\) 代入弱形式方程，并要求该方程对所有试探函数 \(v(x)\) 都成立（在实际计算中，我们让 \(v(x)\) 取遍所有的基函数 \(\phi_i(x)\)）。
这样，我们为每个 \(i = 1, \dots, N\) 得到一个方程： \(\int_0^L \frac{\partial}{\partial t} \left( \sum_{j=1}^N u_j(t) \phi_j(x) \right) \phi_i(x) \, dx + \alpha \int_0^L \frac{\partial}{\partial x} \left( \sum_{j=1}^N u_j(t) \phi_j(x) \right) \frac{\partial \phi_i(x)}{\partial x} \, dx = \int_0^L f \phi_i \, dx\)。
由于求和与积分可交换，且 \(u_j(t)\) 与 \(x\) 无关，我们得到： \(\sum_{j=1}^N \frac{du_j}{dt} \int_0^L \phi_i \phi_j \, dx + \alpha \sum_{j=1}^N u_j \int_0^L \frac{\partial \phi_i}{\partial x} \frac{\partial \phi_j}{\partial x} \, dx = \int_0^L f \phi_i \, dx\)。

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

上一步的推导将一个偏微分方程转化为了一个关于时间 \(t\) 的常微分方程组。

定义：
质量矩阵 \(M\): \(M_{ij} = \int_0^L \phi_i \phi_j \, dx\)
刚度矩阵 \(K\): \(K_{ij} = \alpha \int_0^L \frac{\partial \phi_i}{\partial x} \frac{\partial \phi_j}{\partial x} \, dx\)
载荷向量 \(F(t)\): \(F_i(t) = \int_0^L f(x, t) \phi_i(x) \, dx\)
解向量 \(U(t) = [u_1(t), u_2(t), \dots, u_N(t)]^T\)
上述方程组可以写成简洁的矩阵形式：

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

这是一个线性常微分方程组（如果原问题是线性的）。质量矩阵 \(M\) 和刚度矩阵 \(K\) 都是对称正定的稀疏矩阵（大多数元素为0）。初始条件 \(U(0)\) 可以由初始函数 \(u_0(x)\) 在基函数上的投影确定。

第四步：时间离散——从ODE系统到代数方程组

现在我们需要求解这个关于时间 \(t\) 的常微分方程组。这通常通过时间步进法完成，将连续时间离散为 \(t^0, t^1, \dots, t^n\)。

θ-方法：
- 这是一种通用的单步格式，许多常见格式是它的特例。

思想是在时间区间 \([t^n, t^{n+1}]\) 上，用 \(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\)。
这里，\(U^n\) 是 \(U(t^n)\) 的近似，\(\Delta t = t^{n+1} - t^n\)，参数 \(\theta \in [0, 1]\)。

常见特例：

显式欧拉法 (\(\theta = 0\)): \(M U^{n+1} = (M - \Delta t K) U^n + \Delta t F^n\)。
* 优点：求解简单，只需矩阵乘法。
缺点：通常需要非常小的时间步长 \(\Delta t\) 来保证稳定性（条件稳定）。
隐式欧拉法 (\(\theta = 1\)): \((M + \Delta t K) U^{n+1} = M U^n + \Delta t F^{n+1}\)。
优点：无条件稳定，意味着时间步长 \(\Delta t\) 可以取得较大（受精度限制）。
* 缺点：每一步都需要求解一个线性方程组。
克兰克-尼科尔森方法 (\(\theta = 1/2\)): \((M + \frac{\Delta t}{2} K) U^{n+1} = (M - \frac{\Delta t}{2} K) U^n + \frac{\Delta t}{2} (F^n + F^{n+1})\)。
* 优点：无条件稳定，且时间精度是二阶的（比一阶的隐式/显式欧拉法更精确）。

第五步：求解与总结

在每一步时间推进中，我们需要求解一个形如 \(A U^{n+1} = b\) 的线性代数系统，其中矩阵 \(A = M + \theta \Delta t K\)，右端项 \(b\) 由已知量构成。
由于质量矩阵 \(M\) 和刚度矩阵 \(K\) 是稀疏的，矩阵 \(A\) 通常也是大型稀疏矩阵。因此，高效的迭代法（如共轭梯度法）或直接法（如针对稀疏矩阵的LU分解）被用于求解。
总结流程：抛物型方程的有限元法将一个时空耦合的偏微分方程求解问题，通过“空间离散（有限元）+ 时间离散（差分）”的算子分裂策略，分解为两个相对独立的步骤：1）在空间上用有限元法进行离散，得到关于时间的常微分方程组；2）在时间上用差分法进行离散，最终转化为需要逐步求解的代数方程组。这种方法结合了有限元法处理复杂几何区域的灵活性和差分法推进时间的简洁性。

数值抛物型方程的有限元法好的，我们开始学习“数值抛物型方程的有限元法”。这个方法结合了处理空间导数的有限元离散和处理时间导数的差分离散，是求解与时间相关的扩散、热传导等问题的强大工具。第一步：理解问题原型我们考虑一个典型的线性抛物型方程——一维热传导方程，并配以初边值条件：控制方程 : \(\frac{\partial u}{\partial t} - \alpha \frac{\partial^2 u}{\partial x^2} = f(x, t), \quad x \in (0, L), \quad t \in (0, T ]\) 这里，\(u(x, t)\) 是未知函数（例如温度），\(\alpha > 0\) 是扩散系数，\(f(x, t)\) 是源项。初始条件 : \(u(x, 0) = u_ 0(x), \quad x \in [ 0, L ]\) 这给出了整个区域在初始时刻 \(t=0\) 的状态。边界条件 : 例如狄利克雷边界条件 \(u(0, t) = g_ 0(t), \quad u(L, t) = g_ L(t)\) 这规定了区域边界在所有时刻的状态。我们的目标是在整个时空区域 \([ 0, L] \times [ 0, T ]\) 上找到满足上述条件的函数 \(u(x, t)\) 的近似解。第二步：空间离散的弱形式与有限元逼近有限元法的核心第一步是将偏微分方程转化为积分形式的“弱形式”，从而降低对解的光滑性要求。构造弱形式：将控制方程两边同时乘以一个“试探函数” \(v(x)\)，该函数在满足齐次边界条件的函数空间中选取（例如，若边界条件是 \(u=0\)，则要求 \(v(0)=v(L)=0\)）。在空间域 \([ 0, L]\) 上积分： \(\int_ 0^L \frac{\partial u}{\partial t} v \, dx - \alpha \int_ 0^L \frac{\partial^2 u}{\partial x^2} v \, dx = \int_ 0^L f v \, dx\)。对含二阶导数的项应用分部积分（这是关键一步）： \(-\alpha \int_ 0^L \frac{\partial^2 u}{\partial x^2} v \, dx = \alpha \int_ 0^L \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} \, dx - \alpha \left[ \frac{\partial u}{\partial x} v \right]_ 0^L\)。如果边界条件使得边界项消失（例如 \(v=0\) 在边界上），我们得到弱形式： \(\int_ 0^L \frac{\partial u}{\partial t} v \, dx + \alpha \int_ 0^L \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} \, dx = \int_ 0^L f v \, dx\)。这个形式只要求解 \(u\) 具有一阶导数，比原始方程要求二阶导数更“弱”。有限元离散空间：我们将区间 \([ 0, L ]\) 分割成许多小单元（子区间）。在每个单元上，我们定义简单的多项式函数（最常见的是线性函数），这些函数在单元内部光滑，在单元连接处连续。这些多项式构成了分片多项式函数空间。我们选择一组基函数 \(\{\phi_ 1(x), \phi_ 2(x), \dots, \phi_ N(x)\}\) 来张成这个空间。每个基函数 \(\phi_ i(x)\) 通常只在包含第 \(i\) 个节点的少数几个单元上非零。我们用这组基函数的线性组合来近似真实解： \(u(x, t) \approx u_ h(x, t) = \sum_ {j=1}^N u_ j(t) \phi_ j(x)\)。这里的关键是，我们将解的时空变量分离了：空间变化由基函数 \(\phi_ j(x)\) 描述，而时间变化由待求的系数 \(u_ j(t)\) 描述。伽辽金方法：将近似解 \(u_ h(x, t)\) 代入弱形式方程，并要求该方程对所有试探函数 \(v(x)\) 都成立（在实际计算中，我们让 \(v(x)\) 取遍所有的基函数 \(\phi_ i(x)\)）。这样，我们为每个 \(i = 1, \dots, N\) 得到一个方程： \(\int_ 0^L \frac{\partial}{\partial t} \left( \sum_ {j=1}^N u_ j(t) \phi_ j(x) \right) \phi_ i(x) \, dx + \alpha \int_ 0^L \frac{\partial}{\partial x} \left( \sum_ {j=1}^N u_ j(t) \phi_ j(x) \right) \frac{\partial \phi_ i(x)}{\partial x} \, dx = \int_ 0^L f \phi_ i \, dx\)。由于求和与积分可交换，且 \(u_ j(t)\) 与 \(x\) 无关，我们得到： \(\sum_ {j=1}^N \frac{du_ j}{dt} \int_ 0^L \phi_ i \phi_ j \, dx + \alpha \sum_ {j=1}^N u_ j \int_ 0^L \frac{\partial \phi_ i}{\partial x} \frac{\partial \phi_ j}{\partial x} \, dx = \int_ 0^L f \phi_ i \, dx\)。第三步：得到常微分方程组系统上一步的推导将一个偏微分方程转化为了一个关于时间 \(t\) 的常微分方程组。定义：质量矩阵 \(M\): \(M_ {ij} = \int_ 0^L \phi_ i \phi_ j \, dx\) 刚度矩阵 \(K\): \(K_ {ij} = \alpha \int_ 0^L \frac{\partial \phi_ i}{\partial x} \frac{\partial \phi_ j}{\partial x} \, dx\) 载荷向量 \(F(t)\): \(F_ i(t) = \int_ 0^L f(x, t) \phi_ i(x) \, dx\) 解向量 \(U(t) = [ u_ 1(t), u_ 2(t), \dots, u_ N(t) ]^T\) 上述方程组可以写成简洁的矩阵形式： \[ M \frac{dU}{dt} + K U = F(t) \] 这是一个线性常微分方程组（如果原问题是线性的）。质量矩阵 \(M\) 和刚度矩阵 \(K\) 都是对称正定的稀疏矩阵（大多数元素为0）。初始条件 \(U(0)\) 可以由初始函数 \(u_ 0(x)\) 在基函数上的投影确定。第四步：时间离散——从ODE系统到代数方程组现在我们需要求解这个关于时间 \(t\) 的常微分方程组。这通常通过时间步进法完成，将连续时间离散为 \(t^0, t^1, \dots, t^n\)。 θ-方法：这是一种通用的单步格式，许多常见格式是它的特例。思想是在时间区间 \([ t^n, t^{n+1} ]\) 上，用 \(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\)。这里，\(U^n\) 是 \(U(t^n)\) 的近似，\(\Delta t = t^{n+1} - t^n\)，参数 \(\theta \in [ 0, 1 ]\)。常见特例：显式欧拉法 (\(\theta = 0\)): \(M U^{n+1} = (M - \Delta t K) U^n + \Delta t F^n\)。优点：求解简单，只需矩阵乘法。缺点：通常需要非常小的时间步长 \(\Delta t\) 来保证稳定性（条件稳定）。隐式欧拉法 (\(\theta = 1\)): \((M + \Delta t K) U^{n+1} = M U^n + \Delta t F^{n+1}\)。优点：无条件稳定，意味着时间步长 \(\Delta t\) 可以取得较大（受精度限制）。缺点：每一步都需要求解一个线性方程组。克兰克-尼科尔森方法 (\(\theta = 1/2\)): \((M + \frac{\Delta t}{2} K) U^{n+1} = (M - \frac{\Delta t}{2} K) U^n + \frac{\Delta t}{2} (F^n + F^{n+1})\)。优点：无条件稳定，且时间精度是二阶的（比一阶的隐式/显式欧拉法更精确）。第五步：求解与总结在每一步时间推进中，我们需要求解一个形如 \(A U^{n+1} = b\) 的线性代数系统，其中矩阵 \(A = M + \theta \Delta t K\)，右端项 \(b\) 由已知量构成。由于质量矩阵 \(M\) 和刚度矩阵 \(K\) 是稀疏的，矩阵 \(A\) 通常也是大型稀疏矩阵。因此，高效的迭代法（如共轭梯度法）或直接法（如针对稀疏矩阵的LU分解）被用于求解。总结流程：抛物型方程的有限元法将一个时空耦合的偏微分方程求解问题，通过“空间离散（有限元）+ 时间离散（差分）”的算子分裂策略，分解为两个相对独立的步骤：1）在空间上用有限元法进行离散，得到关于时间的常微分方程组；2）在时间上用差分法进行离散，最终转化为需要逐步求解的代数方程组。这种方法结合了有限元法处理复杂几何区域的灵活性和差分法推进时间的简洁性。