块三对角矩阵的追赶法（Thomas算法）
题目描述
考虑一个块三对角线性方程组 \(A\mathbf{x} = \mathbf{b}\)，其中矩阵 \(A\) 的结构如下：
\[A = \begin{bmatrix} B_1 & C_1 & & & \\ A_2 & B_2 & C_2 & & \\ & A_3 & B_3 & C_3 & \\ & & \ddots & \ddots & \ddots \\ & & & A_m & B_m \end{bmatrix} \]
这里，\(A_i, B_i, C_i\) 都是 \(n \times n\) 的子块矩阵（假设非奇异），\(\mathbf{x}\) 和 \(\mathbf{b}\) 是分块向量，形式为：
\[\mathbf{x} = \begin{bmatrix} \mathbf{x}_1 \\ \mathbf{x}_2 \\ \vdots \\ \mathbf{x}_m \end{bmatrix}, \quad \mathbf{b} = \begin{bmatrix} \mathbf{b}_1 \\ \mathbf{b}_2 \\ \vdots \\ \mathbf{b}_m \end{bmatrix} \]
其中每个 \(\mathbf{x}_i, \mathbf{b}_i\) 是 \(n\) 维向量。块追赶法（Thomas算法）是标准追赶法（用于标量三对角矩阵）在块矩阵情形的推广，旨在高效求解此类特殊结构的线性方程组。本题目将详细讲解块Thomas算法的推导和求解过程。
解题过程
步骤1：算法思想
追赶法的核心思想是将系数矩阵 \(A\) 做块LU分解，即分解为一个块下三角矩阵 \(L\) 和一个块上三角矩阵 \(U\) 的乘积。由于 \(A\) 的特殊块三对角结构，可以假设 \(L\) 和 \(U\) 具有以下形式：
\[L = \begin{bmatrix} I & & & & \\ \hat{A}_2 & I & & & \\ & \hat{A}_3 & I & & \\ & & \ddots & \ddots & \\ & & & \hat{A}_m & I \end{bmatrix}, \quad U = \begin{bmatrix} \hat{B}_1 & C_1 & & & \\ & \hat{B}_2 & C_2 & & \\ & & \hat{B}_3 & C_3 & \\ & & & \ddots & \ddots \\ & & & & \hat{B}_m \end{bmatrix} \]
其中 \(I\) 是 \(n \times n\) 单位矩阵，\(\hat{A}_i\) 和 \(\hat{B}_i\) 是待求的 \(n \times n\) 矩阵。目标是使得 \(A = LU\)。
步骤2：推导块LU分解的递推关系
将 \(L\) 和 \(U\) 相乘，并与 \(A\) 的对应块进行比较：
- 第一行块：\(\hat{B}_1 = B_1\)。
- 第一列块（第二行开始）：\(\hat{A}_2 \hat{B}_1 = A_2\)。由于已得 \(\hat{B}_1 = B_1\)，所以 \(\hat{A}_2 = A_2 \hat{B}_1^{-1}\)。
- 对角线块（i ≥ 2）：由 \(LU\) 乘积的第 i 行、第 i 列块可得：
\(\hat{A}_i C_{i-1} + \hat{B}_i = B_i\)。
因此，\(\hat{B}_i = B_i - \hat{A}_i C_{i-1}\)。 - 上对角线块：\(LU\) 中上对角线块自然就是 \(C_i\)，与 \(A\) 一致，无需额外计算。
- 下对角线块（i ≥ 2）：\(LU\) 中下对角线块为 \(\hat{A}_i \hat{B}_{i-1}\)，应等于 \(A_i\)。但注意，我们在计算 \(\hat{A}_i\) 时实际上是通过前一步的 \(\hat{B}_{i-1}\) 来定义的，即：
\(\hat{A}_i = A_i \hat{B}_{i-1}^{-1}\)。
总结递推公式（追过程，即前向消去）：
\[\begin{aligned} \hat{B}_1 &= B_1, \\ \text{对 } i &= 2, 3, \dots, m: \\ \hat{A}_i &= A_i \hat{B}_{i-1}^{-1}, \\ \hat{B}_i &= B_i - \hat{A}_i C_{i-1}. \end{aligned} \]
这里需要求 \(\hat{B}_{i-1}\) 的逆矩阵，这是算法的核心计算量所在。
步骤3：求解 \(L\mathbf{y} = \mathbf{b}\)（前向代入）
令 \(U\mathbf{x} = \mathbf{y}\)，先解 \(L\mathbf{y} = \mathbf{b}\)。写出方程：
\[\begin{bmatrix} I & & & & \\ \hat{A}_2 & I & & & \\ & \hat{A}_3 & I & & \\ & & \ddots & \ddots & \\ & & & \hat{A}_m & I \end{bmatrix} \begin{bmatrix} \mathbf{y}_1 \\ \mathbf{y}_2 \\ \vdots \\ \mathbf{y}_m \end{bmatrix} = \begin{bmatrix} \mathbf{b}_1 \\ \mathbf{b}_2 \\ \vdots \\ \mathbf{b}_m \end{bmatrix} \]
逐块求解：
\[\begin{aligned} \mathbf{y}_1 &= \mathbf{b}_1, \\ \mathbf{y}_i &= \mathbf{b}_i - \hat{A}_i \mathbf{y}_{i-1}, \quad i = 2, 3, \dots, m. \end{aligned} \]
步骤4：求解 \(U\mathbf{x} = \mathbf{y}\)（后向代入）
现在解：
\[\begin{bmatrix} \hat{B}_1 & C_1 & & & \\ & \hat{B}_2 & C_2 & & \\ & & \hat{B}_3 & C_3 & \\ & & & \ddots & \ddots \\ & & & & \hat{B}_m \end{bmatrix} \begin{bmatrix} \mathbf{x}_1 \\ \mathbf{x}_2 \\ \vdots \\ \mathbf{x}_m \end{bmatrix} = \begin{bmatrix} \mathbf{y}_1 \\ \mathbf{y}_2 \\ \vdots \\ \mathbf{y}_m \end{bmatrix} \]
从最后一个块开始回代：
\[\begin{aligned} \mathbf{x}_m &= \hat{B}_m^{-1} \mathbf{y}_m, \\ \mathbf{x}_i &= \hat{B}_i^{-1} (\mathbf{y}_i - C_i \mathbf{x}_{i+1}), \quad i = m-1, m-2, \dots, 1. \end{aligned} \]
步骤5：算法总结与计算量分析
块Thomas算法流程如下：
-
追过程（分解与前向消去）：
- 计算 \(\hat{B}_1 = B_1\)。
- 对 \(i = 2\) 到 \(m\)：
a. 求解 \(\hat{A}_i\) 满足 \(\hat{A}_i \hat{B}_{i-1} = A_i\)（即 \(\hat{A}_i = A_i \hat{B}_{i-1}^{-1}\)）。
b. 计算 \(\hat{B}_i = B_i - \hat{A}_i C_{i-1}\)。
-
前向代入：
- \(\mathbf{y}_1 = \mathbf{b}_1\)。
- 对 \(i = 2\) 到 \(m\)：\(\mathbf{y}_i = \mathbf{b}_i - \hat{A}_i \mathbf{y}_{i-1}\)。
-
后向代入：
- 求解 \(\hat{B}_m \mathbf{x}_m = \mathbf{y}_m\) 得 \(\mathbf{x}_m\)。
- 对 \(i = m-1\) 到 1：求解 \(\hat{B}_i \mathbf{x}_i = \mathbf{y}_i - C_i \mathbf{x}_{i+1}\) 得 \(\mathbf{x}_i\)。
计算量：每一步都需要求解一个以 \(\hat{B}_i\) 为系数的 \(n \times n\) 线性方程组（即求逆或分解 \(\hat{B}_i\)）。如果 \(n\) 不大，可以直接对每个 \(\hat{B}_i\) 做LU分解（\(O(n^3)\) 运算）；如果 \(n\) 很大，但 \(\hat{B}_i\) 有特殊结构（如稀疏、对角占优），可使用迭代法求解。整体复杂度约为 \(O(m \cdot \text{解一个} n \times n \text{系统的代价})\)。
步骤6：数值稳定性说明
块Thomas算法在 \(\hat{B}_i\) 非奇异时是可行的。数值稳定性依赖于 \(\hat{B}_i\) 的条件数。若 \(A\) 是块对角占优或对称正定，通常 \(\hat{B}_i\) 良态，算法稳定。否则，可能需要部分主元或完全主元策略，但会破坏块三对角结构，增加计算量。
总结
块追赶法将大规模块三对角系统分解为一系列小块矩阵的求解，显著降低计算复杂度，是求解偏微分方程离散化后线性系统的经典方法。核心在于块LU分解的递推计算，以及前后两次代入求解。