5.3矩阵三角分解法

矩阵直接三角分解法（杜利特尔Doolittle分解）详细讲解

各位同学，今天我们系统讲解数值线性代数中核心的直接解法——矩阵直接三角分解法，它是高斯消去法的矩阵形式化与工程优化实现，是求解线性方程组的核心方法之一。我会从核心思想、前提条件、公式推导、唯一性证明、求解流程、算例验证、方法特性全链条展开，最后用表格做系统归纳。

---

一、方法核心思想与基础前提

1. 核心思想

直接三角分解法的本质，是将n阶非奇异系数矩阵A，分解为两个三角矩阵的乘积：

\[\boldsymbol{A} = \boldsymbol{L}\boldsymbol{U} \]

其中：

• \(\boldsymbol{L}\) 是单位下三角矩阵：对角线元素全为1，对角线以上元素全为0；
• \(\boldsymbol{U}\) 是上三角矩阵：对角线以下元素全为0。

一旦完成上述分解，线性方程组 \(\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}\) 就可以等价拆解为两个极易求解的三角方程组：

1. 单位下三角方程组 \(\boldsymbol{L}\boldsymbol{y}=\boldsymbol{b}\)，通过前向替换（前代）求解中间变量 \(\boldsymbol{y}\)；
2. 上三角方程组 \(\boldsymbol{U}\boldsymbol{x}=\boldsymbol{y}\)，通过后向替换（回代）求解最终解 \(\boldsymbol{x}\)。

这个方法的核心优势，是将高斯消去法的“消元+回代”两步，转化为“一次分解+多次前代回代”，尤其适合求解同系数矩阵、多右端项的批量方程组，工程计算中效率提升极为显著。

2. 分解存在且唯一的充要条件

在讲解推导前，必须先明确方法的适用边界：

定理：n阶方阵\(\boldsymbol{A}\)存在唯一的杜利特尔分解（\(\boldsymbol{L}\)单位下三角，\(\boldsymbol{U}\)上三角）的充要条件是：\(\boldsymbol{A}\)的前\(n-1\)阶顺序主子式全不为零，即

\[D_k = \begin{vmatrix} a_{11} & \dots & a_{1k} \\ \vdots & & \vdots \\ a_{k1} & \dots & a_{kk} \end{vmatrix} \neq 0, \quad k=1,2,\dots,n-1 \]

这个定理是所有推导的前提，后续公式中的除法操作，都依赖于顺序主子式非零带来的主元非零特性。

---

二、杜利特尔分解公式的详细推导

我们从矩阵乘法的基本定义出发，一步步推导\(\boldsymbol{L}\)和\(\boldsymbol{U}\)元素的递推计算公式。

设n阶矩阵：

\[\boldsymbol{A}=(a_{ij})_{n\times n}, \quad \boldsymbol{L}=(l_{ij})_{n\times n}, \quad \boldsymbol{U}=(u_{ij})_{n\times n} \]

其中：

• 单位下三角矩阵\(\boldsymbol{L}\)满足：\(l_{ii}=1\)，\(l_{ij}=0 \ (i<j)\)；
• 上三角矩阵\(\boldsymbol{U}\)满足：\(u_{ij}=0 \ (i>j)\)。

由矩阵乘法定义，\(\boldsymbol{A}=\boldsymbol{L}\boldsymbol{U}\)等价于对任意\(i,j \in [1,n]\)，有：

\[a_{ij} = \sum_{k=1}^n l_{ik}u_{kj} \]

结合\(\boldsymbol{L}\)和\(\boldsymbol{U}\)的三角结构，求和项中大量元素为0，我们可以分步骤简化这个等式，得到递推公式。

步骤1：计算\(\boldsymbol{U}\)的第1行、\(\boldsymbol{L}\)的第1列

（1）\(\boldsymbol{U}\)的第1行推导

当\(i=1\)时，\(\boldsymbol{L}\)的第1行只有\(l_{11}=1\)，其余\(l_{1k}=0 \ (k>1)\)，因此：

\[a_{1j} = \sum_{k=1}^n l_{1k}u_{kj} = l_{11}u_{1j} = u_{1j} \]

直接得到\(\boldsymbol{U}\)第1行的计算公式：

\[\boldsymbol{u_{1j} = a_{1j}, \quad j=1,2,\dots,n} \]

即\(\boldsymbol{U}\)的第1行与\(\boldsymbol{A}\)的第1行完全相等。

（2）\(\boldsymbol{L}\)的第1列推导

当\(j=1\)时，\(\boldsymbol{U}\)的第1列只有\(u_{11} \neq 0\)，其余\(u_{k1}=0 \ (k>1)\)，因此对\(i \geq 2\)：

\[a_{i1} = \sum_{k=1}^n l_{ik}u_{k1} = l_{i1}u_{11} \]

由充要条件，1阶顺序主子式\(D_1=u_{11}=a_{11} \neq 0\)，因此可以移项得到\(\boldsymbol{L}\)第1列的计算公式：

\[\boldsymbol{l_{i1} = \frac{a_{i1}}{u_{11}}, \quad i=2,3,\dots,n} \]

步骤2：递推计算\(\boldsymbol{U}\)的第\(r\)行、\(\boldsymbol{L}\)的第\(r\)列（\(r=2,3,\dots,n\)）

我们采用数学归纳法的思路：假设已经计算出\(\boldsymbol{U}\)的前\(r-1\)行、\(\boldsymbol{L}\)的前\(r-1\)列的所有元素，现在推导第\(r\)行的\(\boldsymbol{U}\)和第\(r\)列的\(\boldsymbol{L}\)。

（1）\(\boldsymbol{U}\)的第\(r\)行推导

\(\boldsymbol{U}\)是上三角矩阵，因此第\(r\)行的元素仅当\(j \geq r\)时有非零值。对\(j \geq r\)，展开矩阵乘法：

\[a_{rj} = \sum_{k=1}^n l_{rk}u_{kj} \]

由于\(\boldsymbol{L}\)是单位下三角矩阵，\(k>r\)时\(l_{rk}=0\)，因此求和上限可缩减为\(r\)：

\[a_{rj} = \sum_{k=1}^r l_{rk}u_{kj} \]

将\(k=r\)的项单独拆分出来（\(l_{rr}=1\)）：

\[a_{rj} = \sum_{k=1}^{r-1} l_{rk}u_{kj} + l_{rr}u_{rj} = \sum_{k=1}^{r-1} l_{rk}u_{kj} + u_{rj} \]

移项后直接得到\(\boldsymbol{U}\)第\(r\)行的计算公式：

\[\boldsymbol{u_{rj} = a_{rj} - \sum_{k=1}^{r-1} l_{rk}u_{kj}, \quad j=r,r+1,\dots,n} \]

说明：公式右侧的\(l_{rk}\)（\(k<r\)）是\(\boldsymbol{L}\)第\(r\)行前\(r-1\)个元素，已在之前的步骤中计算完成；\(u_{kj}\)（\(k<r\)）是\(\boldsymbol{U}\)前\(r-1\)行的元素，也已计算完成，因此该式无未知量，可直接计算。

（2）\(\boldsymbol{L}\)的第\(r\)列推导

\(\boldsymbol{L}\)是单位下三角矩阵，因此第\(r\)列的元素仅当\(i > r\)时有非零值。对\(i > r\)，展开矩阵乘法：

\[a_{ir} = \sum_{k=1}^n l_{ik}u_{kr} \]

由于\(\boldsymbol{U}\)是上三角矩阵，\(k>r\)时\(u_{kr}=0\)，因此求和上限可缩减为\(r\)：

\[a_{ir} = \sum_{k=1}^r l_{ik}u_{kr} \]

将\(k=r\)的项单独拆分出来：

\[a_{ir} = \sum_{k=1}^{r-1} l_{ik}u_{kr} + l_{ir}u_{rr} \]

由充要条件，\(r\)阶顺序主子式\(D_r = D_{r-1} \cdot u_{rr} \neq 0\)，因此\(u_{rr} \neq 0\)，移项后得到\(\boldsymbol{L}\)第\(r\)列的计算公式：

\[\boldsymbol{l_{ir} = \frac{1}{u_{rr}} \left( a_{ir} - \sum_{k=1}^{r-1} l_{ik}u_{kr} \right), \quad i=r+1,r+2,\dots,n} \]

说明：公式右侧的\(l_{ik}\)（\(k<r\)）、\(u_{kr}\)（\(k<r\)）均已在之前的步骤中计算完成，\(u_{rr}\)是刚算出的\(\boldsymbol{U}\)第\(r\)行的对角线元素，因此该式可直接计算。

步骤3：分解的计算流程总结

1. 初始化：计算\(\boldsymbol{U}\)第1行\(u_{1j}=a_{1j}\)，\(\boldsymbol{L}\)第1列\(l_{i1}=a_{i1}/u_{11}\)；
2. 递推循环：对\(r=2,3,\dots,n\)，依次执行：
   a. 计算\(\boldsymbol{U}\)的第\(r\)行：\(u_{rj} = a_{rj} - \sum_{k=1}^{r-1} l_{rk}u_{kj} \ (j \geq r)\)；
   b. 计算\(\boldsymbol{L}\)的第\(r\)列：\(l_{ir} = \frac{1}{u_{rr}} \left( a_{ir} - \sum_{k=1}^{r-1} l_{ik}u_{kr} \right) \ (i > r)\)；
3. 循环结束，得到完整的\(\boldsymbol{L}\)和\(\boldsymbol{U}\)。

---

三、三角分解后方程组的求解公式推导

完成\(\boldsymbol{A}=\boldsymbol{L}\boldsymbol{U}\)分解后，方程组\(\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}\)转化为\(\boldsymbol{L}\boldsymbol{y}=\boldsymbol{b}\)和\(\boldsymbol{U}\boldsymbol{x}=\boldsymbol{y}\)两个三角方程组，我们分别推导其求解公式。

1. 单位下三角方程组\(\boldsymbol{L}\boldsymbol{y}=\boldsymbol{b}\)的前向替换公式

\(\boldsymbol{L}\boldsymbol{y}=\boldsymbol{b}\)展开为：

\[\begin{cases} y_1 = b_1 \\ l_{21}y_1 + y_2 = b_2 \\ l_{31}y_1 + l_{32}y_2 + y_3 = b_3 \\ \quad \vdots \\ l_{n1}y_1 + l_{n2}y_2 + \dots + l_{n,n-1}y_{n-1} + y_n = b_n \end{cases} \]

可以看到，第一个方程直接给出\(y_1\)，后续每个方程仅包含一个新的未知量\(y_i\)，其余未知量均已在之前的步骤中算出。

对第\(i\)个方程（\(i \geq 2\)），整理得：

\[y_i + \sum_{k=1}^{i-1} l_{ik}y_k = b_i \]

移项得到前向替换的通用公式：

\[\boldsymbol{ \begin{cases} y_1 = b_1 \\ y_i = b_i - \sum_{k=1}^{i-1} l_{ik}y_k, \quad i=2,3,\dots,n \end{cases} } \]

2. 上三角方程组\(\boldsymbol{U}\boldsymbol{x}=\boldsymbol{y}\)的后向替换公式

\(\boldsymbol{U}\boldsymbol{x}=\boldsymbol{y}\)展开为：

\[\begin{cases} u_{11}x_1 + u_{12}x_2 + \dots + u_{1n}x_n = y_1 \\ u_{22}x_2 + \dots + u_{2n}x_n = y_2 \\ \quad \vdots \\ u_{n-1,n-1}x_{n-1} + u_{n-1,n}x_n = y_{n-1} \\ u_{nn}x_n = y_n \end{cases} \]

可以看到，最后一个方程直接给出\(x_n\)，从后往前每个方程仅包含一个新的未知量\(x_i\)，其余未知量均已算出。

对第\(i\)个方程（\(i \leq n-1\)），整理得：

\[u_{ii}x_i + \sum_{k=i+1}^n u_{ik}x_k = y_i \]

由\(\boldsymbol{A}\)非奇异，\(u_{ii} \neq 0\)，移项得到后向替换的通用公式：

\[\boldsymbol{ \begin{cases} x_n = \frac{y_n}{u_{nn}} \\ x_i = \frac{1}{u_{ii}} \left( y_i - \sum_{k=i+1}^n u_{ik}x_k \right), \quad i=n-1,n-2,\dots,1 \end{cases} } \]

---

四、杜利特尔分解唯一性的严格证明

前面我们给出了分解存在唯一的充要条件，现在给出完整的证明，帮助大家理解方法的理论根基。

定理

n阶方阵\(\boldsymbol{A}\)存在唯一的杜利特尔分解\(\boldsymbol{A}=\boldsymbol{L}\boldsymbol{U}\)（\(\boldsymbol{L}\)单位下三角，\(\boldsymbol{U}\)上三角）的充要条件是\(\boldsymbol{A}\)的前\(n-1\)阶顺序主子式全不为零。

1. 充分性证明（条件满足→分解存在且唯一）

采用数学归纳法证明：

• 基例（n=1）：1阶矩阵\(\boldsymbol{A}=[a_{11}]\)，显然可分解为\(\boldsymbol{A}=[1][a_{11}]\)，\(\boldsymbol{L}=[1]\)是单位下三角，\(\boldsymbol{U}=[a_{11}]\)是上三角，分解存在且唯一。

• 归纳假设：假设对\(n=k-1\)阶矩阵，若其前\(k-2\)阶顺序主子式全不为零，则存在唯一的杜利特尔分解。

• 归纳步骤（n=k）：对k阶矩阵\(\boldsymbol{A}\)，将其分块为：

  \[\boldsymbol{A} = \begin{pmatrix} \boldsymbol{A}_{k-1} & \boldsymbol{c} \\ \boldsymbol{r}^T & a_{kk} \end{pmatrix} \]

  其中\(\boldsymbol{A}_{k-1}\)是\(\boldsymbol{A}\)的\(k-1\)阶顺序主子矩阵，\(\boldsymbol{c}\)是\(k-1\)维列向量，\(\boldsymbol{r}^T\)是\(k-1\)维行向量。

  我们要找的\(\boldsymbol{L}\)和\(\boldsymbol{U}\)可对应分块为：

  \[\boldsymbol{L} = \begin{pmatrix} \boldsymbol{L}_{k-1} & \boldsymbol{0} \\ \boldsymbol{m}^T & 1 \end{pmatrix}, \quad \boldsymbol{U} = \begin{pmatrix} \boldsymbol{U}_{k-1} & \boldsymbol{v} \\ \boldsymbol{0} & u_{kk} \end{pmatrix} \]

  其中\(\boldsymbol{L}_{k-1}\)是\(k-1\)阶单位下三角矩阵，\(\boldsymbol{U}_{k-1}\)是\(k-1\)阶上三角矩阵。

  由\(\boldsymbol{A}=\boldsymbol{L}\boldsymbol{U}\)，展开分块乘法得到四个等式：

  1. \(\boldsymbol{L}_{k-1}\boldsymbol{U}_{k-1} = \boldsymbol{A}_{k-1}\)：由归纳假设，\(\boldsymbol{A}_{k-1}\)的前\(k-2\)阶顺序主子式全不为零，因此存在唯一的\(\boldsymbol{L}_{k-1}\)和\(\boldsymbol{U}_{k-1}\)满足该式；
  2. \(\boldsymbol{L}_{k-1}\boldsymbol{v} = \boldsymbol{c}\)：\(\boldsymbol{L}_{k-1}\)是单位下三角矩阵，行列式为1，非奇异，因此存在唯一的\(\boldsymbol{v}=\boldsymbol{L}_{k-1}^{-1}\boldsymbol{c}\)；
  3. \(\boldsymbol{m}^T\boldsymbol{U}_{k-1} = \boldsymbol{r}^T\)：由条件，\(\boldsymbol{A}_{k-1}\)的行列式（即\(k-1\)阶顺序主子式）\(D_{k-1}=\det(\boldsymbol{U}_{k-1}) \neq 0\)，因此\(\boldsymbol{U}_{k-1}\)非奇异，存在唯一的\(\boldsymbol{m}^T=\boldsymbol{r}^T\boldsymbol{U}_{k-1}^{-1}\)；
  4. \(\boldsymbol{m}^T\boldsymbol{v} + u_{kk} = a_{kk}\)：移项得唯一的\(u_{kk}=a_{kk}-\boldsymbol{m}^T\boldsymbol{v}\)。

  综上，k阶矩阵\(\boldsymbol{A}\)的杜利特尔分解存在且唯一，由数学归纳法，充分性得证。

2. 必要性证明（分解存在且唯一→条件满足）

若\(\boldsymbol{A}\)存在唯一的杜利特尔分解\(\boldsymbol{A}=\boldsymbol{L}\boldsymbol{U}\)，对任意\(k=1,2,\dots,n-1\)，取\(\boldsymbol{A}\)的k阶顺序主子矩阵\(\boldsymbol{A}_k\)，对应\(\boldsymbol{L}\)和\(\boldsymbol{U}\)的k阶顺序主子矩阵\(\boldsymbol{L}_k\)和\(\boldsymbol{U}_k\)，则有：

\[\boldsymbol{A}_k = \boldsymbol{L}_k\boldsymbol{U}_k \]

其中\(\boldsymbol{L}_k\)是单位下三角矩阵，\(\det(\boldsymbol{L}_k)=1\)；\(\boldsymbol{U}_k\)是上三角矩阵，\(\det(\boldsymbol{U}_k)=\prod_{i=1}^k u_{ii}\)。

因此\(\boldsymbol{A}\)的k阶顺序主子式：

\[D_k = \det(\boldsymbol{A}_k) = \det(\boldsymbol{L}_k)\det(\boldsymbol{U}_k) = \prod_{i=1}^k u_{ii} \]

由于分解唯一，\(\boldsymbol{U}\)的前\(n-1\)个对角线元素\(u_{11},u_{22},\dots,u_{n-1,n-1}\)必须全不为零（若存在\(u_{kk}=0\)，\(k \leq n-1\)，则\(\boldsymbol{U}_k\)奇异，分解将不唯一），因此\(D_k \neq 0\)，\(k=1,2,\dots,n-1\)，必要性得证。

---

五、算例详解（对应教材例5.4）

我们用完整的步骤求解教材中的例题，验证上述公式的正确性。

例题：用直接三角分解法求解线性方程组

\[\begin{pmatrix} 1 & 2 & 3 \\ 2 & 5 & 2 \\ 3 & 1 & 5 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 14 \\ 18 \\ 20 \end{pmatrix} \]

步骤1：验证分解条件

计算前\(n-1=2\)阶顺序主子式：

• 1阶顺序主子式\(D_1=|1|=1 \neq 0\)；
• 2阶顺序主子式\(D_2=\begin{vmatrix}1&2\\2&5\end{vmatrix}=5-4=1 \neq 0\)；
  满足分解条件，存在唯一的杜利特尔分解。

步骤2：计算\(\boldsymbol{L}\)和\(\boldsymbol{U}\)

1. 第1步：计算\(\boldsymbol{U}\)第1行、\(\boldsymbol{L}\)第1列

   • \(\boldsymbol{U}\)第1行：\(u_{11}=1, u_{12}=2, u_{13}=3\)；
   • \(\boldsymbol{L}\)第1列：\(l_{21}=a_{21}/u_{11}=2/1=2\)，\(l_{31}=a_{31}/u_{11}=3/1=3\)；

2. 第2步（r=2）：计算\(\boldsymbol{U}\)第2行、\(\boldsymbol{L}\)第2列

   • \(\boldsymbol{U}\)第2行：
     \(u_{22}=a_{22}-l_{21}u_{12}=5-2×2=1\)；
     \(u_{23}=a_{23}-l_{21}u_{13}=2-2×3=-4\)；
   • \(\boldsymbol{L}\)第2列：
     \(l_{32}=\frac{1}{u_{22}}(a_{32}-l_{31}u_{12})=\frac{1}{1}(1-3×2)=-5\)；

3. 第3步（r=3）：计算\(\boldsymbol{U}\)第3行

   • \(u_{33}=a_{33}-(l_{31}u_{13}+l_{32}u_{23})=5-(3×3 + (-5)×(-4))=5-29=-24\)；

最终得到分解结果：

\[\boldsymbol{L}=\begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & -5 & 1 \end{pmatrix}, \quad \boldsymbol{U}=\begin{pmatrix} 1 & 2 & 3 \\ 0 & 1 & -4 \\ 0 & 0 & -24 \end{pmatrix} \]

验证：\(\boldsymbol{L}\boldsymbol{U}=\boldsymbol{A}\)，分解正确。

步骤3：前向替换求解\(\boldsymbol{L}\boldsymbol{y}=\boldsymbol{b}\)

\(\boldsymbol{b}=(14,18,20)^T\)，代入前代公式：

• \(y_1 = b_1 =14\)；
• \(y_2 = b_2 - l_{21}y_1 = 18-2×14=-10\)；
• \(y_3 = b_3 - (l_{31}y_1 + l_{32}y_2) = 20 - (3×14 + (-5)×(-10))=20-92=-72\)；
  得到\(\boldsymbol{y}=(14,-10,-72)^T\)。

步骤4：后向替换求解\(\boldsymbol{U}\boldsymbol{x}=\boldsymbol{y}\)

代入回代公式：

• \(x_3 = y_3/u_{33} = (-72)/(-24)=3\)；
• \(x_2 = (y_2 - u_{23}x_3)/u_{22} = (-10 - (-4)×3)/1=2\)；
• \(x_1 = (y_1 - (u_{12}x_2 + u_{13}x_3))/u_{11} = (14 - (2×2+3×3))/1=1\)；

最终得到方程组的解：\(\boldsymbol{x}=(1,2,3)^T\)，代入原方程验证，等式成立。

---

六、方法的核心特性与工程实现要点

1. 计算量：杜利特尔分解的乘除法计算量约为\(\frac{n^3}{3}\)，与高斯消去法完全一致；但对于多右端项方程组，每新增一个右端项，仅需增加\(n^2\)次乘除法（前代+回代），远低于高斯消去法的重复\(\frac{n^3}{3}\)计算，工程优势显著。
2. 存储优化：由于\(\boldsymbol{L}\)的对角线元素全为1，无需额外存储，因此\(\boldsymbol{L}\)的下三角元素和\(\boldsymbol{U}\)的上三角元素，可直接覆盖存储在原矩阵\(\boldsymbol{A}\)的存储空间中，无需额外开辟内存，适合大规模矩阵计算。
3. 数值稳定性：不选主元的直接三角分解法，数值稳定性与顺序高斯消去法一致，当出现小主元时，会产生严重的舍入误差，因此实际工程计算中，通常采用列选主元的三角分解法，保证数值稳定性。
4. 适用场景：特别适合求解大规模、同系数矩阵、多右端项的线性方程组，如有限元计算、电路仿真、参数优化等工程领域。

---

七、核心知识点归纳总结表

模块	核心内容	详细公式/说明
方法本质	高斯消去法的矩阵形式，将\(\boldsymbol{A}\)分解为单位下三角矩阵\(\boldsymbol{L}\)和上三角矩阵\(\boldsymbol{U}\)的乘积	\(\boldsymbol{A}=\boldsymbol{L}\boldsymbol{U}\)，将\(\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}\)转化为\(\boldsymbol{L}\boldsymbol{y}=\boldsymbol{b}\)和\(\boldsymbol{U}\boldsymbol{x}=\boldsymbol{y}\)
分解存在唯一的充要条件	保证分解可计算、结果唯一的前提	\(\boldsymbol{A}\)的前\(n-1\)阶顺序主子式全不为零，即\(D_k \neq 0, \ k=1,2,\dots,n-1\)
杜利特尔分解递推公式	计算\(\boldsymbol{L}\)和\(\boldsymbol{U}\)的核心公式	1. 初始行/列：\(u_{1j}=a_{1j} \ (j=1,\dots,n)\)；\(l_{i1}=a_{i1}/u_{11} \ (i=2,\dots,n)\) 2. \(\boldsymbol{U}\)第\(r\)行：\(u_{rj}=a_{rj}-\sum_{k=1}^{r-1}l_{rk}u_{kj} \ (j=r,\dots,n)\) 3. \(\boldsymbol{L}\)第\(r\)列：\(l_{ir}=\frac{1}{u_{rr}}\left(a_{ir}-\sum_{k=1}^{r-1}l_{ik}u_{kr}\right) \ (i=r+1,\dots,n)\)
前代求解公式	单位下三角方程组\(\boldsymbol{L}\boldsymbol{y}=\boldsymbol{b}\)的求解	\(\begin{cases} y_1 = b_1 \\ y_i = b_i - \sum_{k=1}^{i-1}l_{ik}y_k \ (i=2,\dots,n) \end{cases}\)
回代求解公式	上三角方程组\(\boldsymbol{U}\boldsymbol{x}=\boldsymbol{y}\)的求解	\(\begin{cases} x_n = y_n/u_{nn} \\ x_i = \frac{1}{u_{ii}}\left(y_i - \sum_{k=i+1}^n u_{ik}x_k\right) \ (i=n-1,\dots,1) \end{cases}\)
计算量	算法的时间复杂度	1. 分解过程：约\(\frac{n^3}{3}\)次乘除法 2. 单右端项前代+回代：约\(n^2\)次乘除法
核心优势	工程应用的核心价值	1. 一次分解可重复用于多右端项方程组，批量计算效率极高 2. 存储可复用，无需额外内存 3. 公式规整，易于编程实现
局限性	方法的适用边界	1. 不选主元时，小主元会导致严重舍入误差，数值稳定性不足 2. 顺序主子式为零时，无法直接分解，需选主元处理

---

以上就是直接三角分解法的完整讲解，从理论根基到公式推导，再到算例验证和工程特性，形成了完整的知识体系。这个方法是数值线性代数的核心基础，后续的平方根法、追赶法等特殊矩阵的分解方法，都是基于这个核心框架延伸而来的。

列选主元的三角分解法（部分选主元LU分解）完整讲解

承接上一讲的不选主元直接三角分解，本次我们系统讲解工程中真正实用的列选主元三角分解法，它解决了不选主元方法的核心缺陷，是工业级数值计算中求解线性方程组的标准算法。我会从提出背景、核心原理、公式推导、完整算法、算例验证、拓展应用（求逆矩阵）、数值特性全链条展开，最后用表格做系统归纳。

---

一、方法提出的核心背景：解决不选主元分解的致命缺陷

不选主元的杜利特尔分解（直接三角分解），在理论上是高斯消去法的矩阵形式，但在实际数值计算中存在两个不可忽视的致命问题：

1. 计算中断风险：当分解过程中某一步的主元\(u_{rr}=0\)时，分解公式会出现除以0的操作，计算直接中断。即使矩阵\(\boldsymbol{A}\)是非奇异（可逆）的，只要其前\(n-1\)阶顺序主子式有一个为0，就无法完成直接分解，例如\(\boldsymbol{A}=\begin{pmatrix}0&1\\1&0\end{pmatrix}\)，可逆但无法直接做LU分解。
2. 舍入误差严重放大：当主元\(u_{rr}\)的绝对值极小时，计算\(l_{ir}=s_i/u_{rr}\)时，相当于用极小的数做除数，会把\(s_i\)的舍入误差急剧放大，最终导致解的严重失真，甚至完全错误。

核心解决思路与理论支撑

对于n阶非奇异矩阵\(\boldsymbol{A}\)，一定存在排列矩阵\(\boldsymbol{P}\)（行交换矩阵），使得行交换后的矩阵\(\boldsymbol{PA}\)可以做唯一的杜利特尔分解：

\[\boldsymbol{PA}=\boldsymbol{LU} \]

其中：

• \(\boldsymbol{L}\)是单位下三角矩阵，且所有非对角线元素满足\(|l_{ij}|\leq1\)，从根源上避免小除数的误差放大；
• \(\boldsymbol{U}\)是上三角矩阵；
• 该方法与列主元高斯消去法完全等价，但具备“一次分解、多次求解”的效率优势。

---

二、选主元三角分解的核心公式与详细推导

我们以杜利特尔分解为基础，推导列选主元的分解逻辑，核心是在计算主元前先做选主元+行交换，保证主元绝对值最大。

1. 第r步分解的基础回顾

在第r步分解时，我们已经计算出\(\boldsymbol{U}\)的前\(r-1\)行、\(\boldsymbol{L}\)的前\(r-1\)列，接下来需要计算\(\boldsymbol{U}\)的第r行和\(\boldsymbol{L}\)的第r列。
不选主元时，\(\boldsymbol{L}\)第r列的计算公式为：

\[l_{ir} = \frac{a_{ir} - \sum_{k=1}^{r-1} l_{ik}u_{kr}}{u_{rr}}, \quad i=r+1,\dots,n \]

其中分母\(u_{rr}\)就是主元，也是误差放大的核心来源。

2. 中间量\(s_i\)的定义（选主元的预处理）

在第r步，我们先对第r列、第r行到第n行的所有元素，提前完成求和项的计算，定义中间量：

\[\boldsymbol{s_i = a_{ir} - \sum_{k=1}^{r-1} l_{ik}u_{kr}, \quad i=r,r+1,\dots,n} \]

这个\(s_i\)的物理意义：

• 当\(i=r\)时，\(s_r\)就是原本要计算的主元\(u_{rr}\)；
• 当\(i>r\)时，\(s_i\)就是\(l_{ir}\)的分子部分，原本需要除以\(u_{rr}\)得到\(l_{ir}\)。

提前计算\(s_i\)后，我们就可以在确定主元前，先选出当前列的最优主元。

3. 列选主元的规则

我们在所有\(s_i\)（\(i=r\)到\(n\)）中，选出绝对值最大的元素，记为：

\[|s_{i_r}| = \max_{r\leq i\leq n} |s_i| \]

其中\(i_r\)就是我们选定的主元行号。

4. 行交换操作

交换矩阵\(\boldsymbol{A}\)的第\(r\)行和第\(i_r\)行的所有元素，同时用整型数组\(Ip(r)\)记录本次交换的行号（\(Ip(r)=i_r\)），后续处理右端项和求逆矩阵时会用到。
行交换后，原本的\(s_{i_r}\)被换到\((r,r)\)位置，成为新的主元\(s_r\)，保证主元的绝对值是当前列最大的。

5. 分解计算（行交换后）

行交换完成后，即可安全计算\(\boldsymbol{U}\)的第r行和\(\boldsymbol{L}\)的第r列：

1. 主元赋值：\(u_{rr}=s_r\)（即交换后的\(a_{rr}\)）；
2. \(\boldsymbol{L}\)的第r列：\(l_{ir} = \frac{s_i}{u_{rr}}, \quad i=r+1,\dots,n\)
   由于我们选了绝对值最大的\(s_r\)作为主元，因此天然满足\(|l_{ir}|\leq1\)，完全避免了小除数导致的误差放大；
3. \(\boldsymbol{U}\)的第r行：\(u_{ri} = a_{ri} - \sum_{k=1}^{r-1} l_{rk}u_{ki}, \quad i=r+1,\dots,n\)
   与不选主元的公式完全一致，行交换后\(l_{rk}\)、\(u_{ki}\)均已计算完成，可直接求解。

工程实现优化：\(\boldsymbol{L}\)的下三角元素（不含对角线，对角线恒为1）可直接覆盖存储在\(\boldsymbol{A}\)的下三角区域，\(\boldsymbol{U}\)的上三角元素可覆盖存储在\(\boldsymbol{A}\)的上三角区域，无需额外开辟内存。

---

三、选主元三角分解的完整算法详解

我们将算法拆分为PA=LU分解阶段和方程组求解阶段，每一步明确操作、公式和意义，完全对应教材中的算法5.2。

算法：列选主元三角分解法求解\(\boldsymbol{Ax=b}\)

已知：n阶非奇异矩阵\(\boldsymbol{A}\)，右端项\(\boldsymbol{b}\)；
输出：方程组的解\(\boldsymbol{x}\)（存储在\(\boldsymbol{b}\)中），\(\boldsymbol{PA=LU}\)的分解结果（存储在\(\boldsymbol{A}\)中），排列矩阵\(\boldsymbol{P}\)（记录在\(Ip\)数组中）。

第一阶段：\(\boldsymbol{PA=LU}\)选主元分解（核心）

对\(r=1,2,\dots,n\)，依次执行以下4步：

1. 计算中间量\(s_i\)
   对\(i=r,r+1,\dots,n\)，计算：

   \[a_{ir} \leftarrow s_i = a_{ir} - \sum_{k=1}^{r-1} a_{ik}a_{kr} \]

   （工程实现中，\(a_{ik}\)存储\(l_{ik}\)，\(a_{kr}\)存储\(u_{kr}\)，直接覆盖存储，无需额外内存）
   意义：提前完成主元列的求和计算，为选主元做准备。

2. 列选主元
   在\(i=r\)到\(n\)的范围内，找到绝对值最大的\(a_{ir}\)对应的行号\(i_r\)：

   \[|a_{i_r r}| = \max_{r\leq i\leq n} |a_{ir}| \]

   记录主元行号：\(Ip(r) \leftarrow i_r\)
   意义：选出当前列的最大主元，保证数值稳定性。

3. 行交换操作
   交换\(\boldsymbol{A}\)的第\(r\)行和第\(i_r\)行的所有元素：

   \[a_{rj} \leftrightarrow a_{i_r j}, \quad j=1,2,\dots,n \]

   意义：将最大主元换到\((r,r)\)位置，作为分解的主元。

4. 计算\(\boldsymbol{L}\)的第r列和\(\boldsymbol{U}\)的第r行
   a. 计算\(\boldsymbol{L}\)的第r列（\(r\neq n\)时执行，最后一行无\(\boldsymbol{L}\)的非对角线元素）：
   对\(i=r+1,\dots,n\)，计算：
   $$a_{ir} \leftarrow l_{ir} = \frac{a_{ir}}{a_{rr}}$$
   天然满足\(|l_{ir}|\leq1\)，控制元素大小。
   b. 计算\(\boldsymbol{U}\)的第r行：
   对\(i=r+1,\dots,n\)，计算：
   $$a_{ri} \leftarrow u_{ri} = a_{ri} - \sum_{k=1}^{r-1} a_{rk}a_{ki}$$
   意义：完成第r步的分解，得到\(\boldsymbol{U}\)的第r行和\(\boldsymbol{L}\)的第r列。

分解完成后：

• \(\boldsymbol{A}\)的上三角部分（含对角线）存储\(\boldsymbol{U}\)的全部元素；
• \(\boldsymbol{A}\)的下三角部分（不含对角线）存储\(\boldsymbol{L}\)的全部非对角线元素；
• \(Ip\)数组记录了每一步的行交换信息，对应排列矩阵\(\boldsymbol{P}\)。

第二阶段：求解方程组\(\boldsymbol{Ly=Pb}\)、\(\boldsymbol{Ux=y}\)

分解完成后，\(\boldsymbol{Ax=b}\)等价于\(\boldsymbol{PAx=Pb}\)，即\(\boldsymbol{LUx=Pb}\)，分三步求解：

1. 对右端项\(\boldsymbol{b}\)做行交换，得到\(\boldsymbol{Pb}\)
   对\(i=1,2,\dots,n-1\)，依次执行：
   a. \(t \leftarrow Ip(i)\)；
   b. 若\(i\neq t\)，交换\(\boldsymbol{b}\)的第\(i\)个和第\(t\)个元素：\(b_i \leftrightarrow b_t\)。
   意义：将排列矩阵\(\boldsymbol{P}\)作用到\(\boldsymbol{b}\)上，与分解后的\(\boldsymbol{PA}\)匹配。

2. 前向替换求解\(\boldsymbol{Ly=Pb}\)
   对\(i=2,3,\dots,n\)，计算：

   \[b_i \leftarrow b_i - \sum_{k=1}^{i-1} a_{ik}b_k \]

   （\(a_{ik}\)为\(l_{ik}\)）
   计算完成后，\(\boldsymbol{b}\)中存储中间解\(\boldsymbol{y}\)。

3. 后向替换求解\(\boldsymbol{Ux=y}\)
   a. 计算最后一个元素：\(b_n \leftarrow \frac{b_n}{a_{nn}}\)；
   b. 对\(i=n-1,n-2,\dots,1\)，倒序计算：
   $$b_i \leftarrow \frac{1}{a_{ii}} \left( b_i - \sum_{k=i+1}^n a_{ik}b_k \right)$$
   计算完成后，\(\boldsymbol{b}\)中存储方程组的最终解\(\boldsymbol{x}\)。

---

四、拓展应用：用选主元三角分解计算矩阵的逆矩阵

矩阵求逆是工程中的高频需求，选主元三角分解是大规模矩阵求逆的标准方法，计算效率远高于伴随矩阵法。

1. 核心公式推导

由\(\boldsymbol{PA=LU}\)，两边同时右乘\(\boldsymbol{A^{-1}}\)，得：

\[\boldsymbol{PAA^{-1}=LUA^{-1}} \implies \boldsymbol{P=LUA^{-1}} \]

两边依次左乘\(\boldsymbol{U^{-1}}\)、\(\boldsymbol{L^{-1}}\)，得逆矩阵的核心计算公式：

\[\boldsymbol{A^{-1}=U^{-1}L^{-1}P} \]

2. 求逆的详细步骤

1. 对\(\boldsymbol{A}\)做列选主元三角分解，得到\(\boldsymbol{PA=LU}\)，存储\(\boldsymbol{L}\)、\(\boldsymbol{U}\)于\(\boldsymbol{A}\)中，\(Ip\)数组记录行交换信息；
2. 计算上三角矩阵\(\boldsymbol{U}\)的逆矩阵\(\boldsymbol{U^{-1}}\)（上三角矩阵的逆仍为上三角矩阵，可通过回代法高效计算）；
3. 计算单位下三角矩阵\(\boldsymbol{L}\)的逆矩阵\(\boldsymbol{L^{-1}}\)，再计算矩阵乘积\(\boldsymbol{M=U^{-1}L^{-1}}\)；
4. 对\(\boldsymbol{M}\)的列做交换（与\(Ip\)数组的行交换规则逆序对应），最终得到\(\boldsymbol{A^{-1}}\)。

3. 计算量说明

该方法求逆的总乘除法计算量约为\(n^3\)次，是大规模矩阵求逆的最优算法之一，LAPACK、NumPy等主流数值库的矩阵求逆均基于该方法实现。

---

五、算法数值特性与等价性说明

1. 数值稳定性
   列选主元三角分解通过每一步选择当前列绝对值最大的主元，保证了\(\boldsymbol{L}\)的所有非对角线元素\(|l_{ij}|\leq1\)，从根源上避免了“大数除以小数”导致的舍入误差放大，是数值稳定的算法，在绝大多数工程场景下都能得到精度可靠的解。

2. 与列主元高斯消去法的等价性
   可严格证明，列选主元三角分解与列主元高斯消去法完全等价：

   • 列主元高斯消去法的消元乘子，就是\(\boldsymbol{L}\)的非对角线元素；
   • 消元过程的行交换，对应排列矩阵\(\boldsymbol{P}\)；
   • 消元得到的上三角矩阵，就是分解中的\(\boldsymbol{U}\)。
     两者的计算量、数值稳定性完全一致，核心区别在于：三角分解法将“消元”与“回代”分离，一次分解可重复用于多个右端项的方程组求解，批量计算效率远高于高斯消去法。

---

六、算例演示（完整算法流程）

我们用一个含小主元的典型案例，验证选主元方法的优势（不选主元会出现严重误差放大）。
例题：用列选主元三角分解法求解方程组

\[\begin{pmatrix} 0.001 & 2.000 & 3.000 \\ 1.000 & 3.712 & 4.623 \\ 2.000 & 1.072 & 5.643 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 1.000 \\ 2.000 \\ 3.000 \end{pmatrix} \]

步骤1：分解阶段

初始化：

\[\boldsymbol{A} = \begin{pmatrix} 0.001 & 2.000 & 3.000 \\ 1.000 & 3.712 & 4.623 \\ 2.000 & 1.072 & 5.643 \end{pmatrix}, \quad Ip=[0,0,0] \]

r=1（第一步分解）

1. 计算\(s_i\)：\(r=1\)，求和项为0，\(s_i=a_{i1}\)，即\(s_1=0.001, s_2=1.000, s_3=2.000\)；
2. 选主元：\(\max|s_i|=2.000\)，对应行号\(i_r=3\)，\(Ip(1)=3\)；
3. 行交换：交换第1行和第3行，\(\boldsymbol{A}\)变为：
   \[\boldsymbol{A} = \begin{pmatrix} 2.000 & 1.072 & 5.643 \\ 1.000 & 3.712 & 4.623 \\ 0.001 & 2.000 & 3.000 \end{pmatrix} \]

4. 计算\(\boldsymbol{L}\)第1列、\(\boldsymbol{U}\)第1行：
   • \(l_{21}=1.000/2.000=0.5\)，\(l_{31}=0.001/2.000=0.0005\)；
   • \(\boldsymbol{U}\)第1行与\(\boldsymbol{A}\)第1行一致；
     此时\(\boldsymbol{A}\)变为：
   \[\boldsymbol{A} = \begin{pmatrix} 2.000 & 1.072 & 5.643 \\ 0.5 & 3.712 & 4.623 \\ 0.0005 & 2.000 & 3.000 \end{pmatrix} \]

r=2（第二步分解）

1. 计算\(s_i\)：
   \(s_2=3.712 - 0.5\times1.072=3.176\)，\(s_3=2.000 - 0.0005\times1.072=1.999464\)；
2. 选主元：\(\max|s_i|=3.176\)，对应行号\(i_r=2\)，\(Ip(2)=2\)，无需行交换；
3. 计算\(\boldsymbol{L}\)第2列、\(\boldsymbol{U}\)第2行：
   • \(l_{32}=1.999464/3.176\approx0.629554\)；
   • \(u_{23}=4.623 - 0.5\times5.643=1.8015\)；
     此时\(\boldsymbol{A}\)变为：
   \[\boldsymbol{A} = \begin{pmatrix} 2.000 & 1.072 & 5.643 \\ 0.5 & 3.176 & 1.8015 \\ 0.0005 & 0.629554 & 3.000 \end{pmatrix} \]

r=3（第三步分解）

1. 计算\(s_3=3.000 - (0.0005\times5.643 + 0.629554\times1.8015)\approx1.8630285\)；
2. 选主元：仅\(i=3\)，\(Ip(3)=3\)，无需行交换；
3. 主元赋值：\(u_{33}=1.8630285\)，分解完成。

最终分解结果：

\[\boldsymbol{L}=\begin{pmatrix}1&0&0\\0.5&1&0\\0.0005&0.629554&1\end{pmatrix}, \quad \boldsymbol{U}=\begin{pmatrix}2.000&1.072&5.643\\0&3.176&1.8015\\0&0&1.8630285\end{pmatrix} \]

步骤2：求解阶段

右端项\(\boldsymbol{b}=(1.000,2.000,3.000)^T\)

1. 行交换：\(Ip(1)=3\)，交换\(b_1\)和\(b_3\)，得\(\boldsymbol{b}=(3.000,2.000,1.000)^T\)；
2. 前向替换：\(y_1=3.000\)，\(y_2=0.5\)，\(y_3\approx0.683723\)；
3. 后向替换：\(x_3\approx0.3670\)，\(x_2\approx-0.0507\)，\(x_1\approx0.4916\)。

最终解：\(\boldsymbol{x}\approx(0.4916, -0.0507, 0.3670)^T\)，代入原方程验证精度极高，而不选主元的方法解的误差会超过100%，充分体现了选主元的优势。

---

七、核心知识点归纳总结表

模块	不选主元直接三角分解	列选主元三角分解
核心分解式	\(\boldsymbol{A=LU}\)	\(\boldsymbol{PA=LU}\)（\(\boldsymbol{P}\)为排列矩阵）
适用条件	\(\boldsymbol{A}\)的前\(n-1\)阶顺序主子式全不为0	仅要求\(\boldsymbol{A}\)非奇异（可逆），无顺序主子式要求
主元选择	固定使用\((r,r)\)位置的元素作为主元	每一步选第r列\(r\sim n\)行中绝对值最大的元素作为主元
\(\boldsymbol{L}\)元素特性	非对角线元素绝对值可能远大于1，误差易放大	非对角线元素满足$
数值稳定性	数值不稳定，小主元会导致解严重失真	数值稳定，是工业级计算的标准算法
计算中断风险	主元为0时计算直接中断	非奇异矩阵不会出现计算中断
计算量	分解阶段约\(n^3/3\)次乘除法，单右端项求解约\(n^2\)次	与不选主元方法完全一致，选主元仅增加极小的比较开销
核心优势	公式简单，易于理解	适用范围广、数值稳定，一次分解可批量求解多右端项方程组
工程应用	仅用于理论教学，无实际工程应用	有限元、CFD、电路仿真、优化求解等所有工程数值计算场景
拓展能力	仅能求解顺序主子式非零的方程组	可求解任意非奇异方程组，可高效计算矩阵的逆矩阵

---

平方根法（Cholesky分解）与改进平方根法完整讲解

本次我们讲解针对对称正定矩阵的专属高效解法——平方根法与改进平方根法，这是有限元分析、结构力学、优化计算等工程领域求解线性方程组的核心算法，是三角分解法在对称正定矩阵上的极致优化。我会从理论基础、公式推导、算法流程、数值特性、算例验证全链条展开，最后用表格完成核心知识点的系统归纳。

---

一、方法背景与核心前提

1. 应用场景

在有限元法求解结构力学、热传导、电磁场等工程问题时，最终得到的线性方程组\(\boldsymbol{Ax=b}\)的系数矩阵\(\boldsymbol{A}\)，绝大多数都具有对称正定的性质。针对这类特殊矩阵，通用的LU分解法没有利用对称性，存在计算和存储的冗余，而平方根法正是为对称正定矩阵量身定制的高效解法。

2. 核心前提：对称正定矩阵的核心性质

n阶矩阵\(\boldsymbol{A}\)为对称正定矩阵，满足：

1. 对称性：\(\boldsymbol{A}=\boldsymbol{A}^T\)；
2. 正定性：对任意非零n维向量\(\boldsymbol{x}\)，都有\(\boldsymbol{x}^T\boldsymbol{Ax}>0\)；
3. 等价性质：\(\boldsymbol{A}\)的所有顺序主子式\(D_k>0\)（\(k=1,2,\dots,n\)），所有特征值均为正数。

这一性质是我们所有分解和推导的理论基础，保证了分解的存在性、唯一性和数值稳定性。

---

二、对称矩阵的LDLᵀ分解（平方根法的理论基础）

我们先从通用对称矩阵的三角分解出发，再延伸到对称正定矩阵的Cholesky分解。

1. 分解定理推导

定理5.9（对称矩阵的三角分解定理）：设\(\boldsymbol{A}\)为n阶对称矩阵，且\(\boldsymbol{A}\)的所有顺序主子式均不为零，则\(\boldsymbol{A}\)可唯一分解为：

\[\boldsymbol{A=LDL^T} \]

其中，\(\boldsymbol{L}\)为单位下三角矩阵（对角线元素全为1），\(\boldsymbol{D}\)为对角矩阵。

详细推导过程

1. 由LU分解的前提，\(\boldsymbol{A}\)的所有顺序主子式不为零，因此存在唯一的杜利特尔分解：
   \[\boldsymbol{A=LU} \]
   其中\(\boldsymbol{L}\)为单位下三角矩阵，\(\boldsymbol{U}\)为上三角矩阵。
2. 对\(\boldsymbol{U}\)做对角拆分，将对角线元素提取出来，分解为对角矩阵和单位上三角矩阵的乘积：
   \[\boldsymbol{U} = \begin{pmatrix} u_{11} & & & \\ & u_{22} & & \\ & & \ddots & \\ & & & u_{nn} \end{pmatrix} \begin{pmatrix} 1 & \frac{u_{12}}{u_{11}} & \dots & \frac{u_{1n}}{u_{11}} \\ & 1 & \dots & \frac{u_{2n}}{u_{22}} \\ & & \ddots & \vdots \\ & & & 1 \end{pmatrix} = \boldsymbol{DU_0} \]
   其中\(\boldsymbol{D}\)为对角矩阵，\(\boldsymbol{U_0}\)为单位上三角矩阵。
3. 代入LU分解式，得到：
   \[\boldsymbol{A=LDU_0} \]

4. 利用\(\boldsymbol{A}\)的对称性\(\boldsymbol{A=A^T}\)，对等式两边转置：
   \[\boldsymbol{A^T=(LDU_0)^T=U_0^T D^T L^T=U_0^T D L^T} \]
   （对角矩阵转置等于自身，即\(\boldsymbol{D^T=D}\)）
5. 由分解的唯一性，杜利特尔分解是唯一的，因此对比\(\boldsymbol{A=LDU_0}\)和\(\boldsymbol{A=U_0^T D L^T}\)，可得：
   \[\boldsymbol{U_0^T=L} \implies \boldsymbol{U_0=L^T} \]

6. 代回原式，最终得到对称矩阵的唯一分解式：
   \[\boldsymbol{A=LDL^T} \]

---

三、平方根法（Cholesky分解）

针对对称正定矩阵，我们可以对LDLᵀ分解做进一步简化，得到平方根法。

1. 分解定理

定理5.10（对称正定矩阵的Cholesky分解）：如果\(\boldsymbol{A}\)为n阶对称正定矩阵，则存在一个实的非奇异下三角矩阵\(\boldsymbol{L}\)，使得：

\[\boldsymbol{A=LL^T} \]

当限定\(\boldsymbol{L}\)的对角元素为正数时，这种分解是唯一的。

详细推导

1. 由对称正定矩阵的性质，\(\boldsymbol{A}\)的所有顺序主子式\(D_k>0\)，因此满足LDLᵀ分解的条件，存在唯一分解\(\boldsymbol{A=LDL^T}\)，其中\(\boldsymbol{L}\)为单位下三角矩阵，\(\boldsymbol{D}\)为对角矩阵。
2. 证明\(\boldsymbol{D}\)的对角元素全为正数：
   由顺序主子式与对角元的关系，\(d_1=D_1>0\)，\(d_i=\frac{D_i}{D_{i-1}}>0\)（\(i=2,3,\dots,n\)），因此\(\boldsymbol{D}\)的所有对角元\(d_i>0\)。
3. 对\(\boldsymbol{D}\)做开方分解：
   \[\boldsymbol{D} = \begin{pmatrix} d_1 & & \\ & \ddots & \\ & & d_n \end{pmatrix} = \begin{pmatrix} \sqrt{d_1} & & \\ & \ddots & \\ & & \sqrt{d_n} \end{pmatrix} \begin{pmatrix} \sqrt{d_1} & & \\ & \ddots & \\ & & \sqrt{d_n} \end{pmatrix} = \boldsymbol{D}^{\frac{1}{2}} \boldsymbol{D}^{\frac{1}{2}} \]

4. 代入LDLᵀ分解式：
   \[\boldsymbol{A=L D^{\frac{1}{2}} D^{\frac{1}{2}} L^T = (L D^{\frac{1}{2}}) (L D^{\frac{1}{2}})^T} \]

5. 令\(\boldsymbol{L_1 = L D^{\frac{1}{2}}}\)，则\(\boldsymbol{L_1}\)是对角元为正的下三角矩阵，最终得到：
   \[\boldsymbol{A=L_1 L_1^T} \]
   这就是Cholesky分解，也叫平方根分解，对应的求解方法就是平方根法。

2. Cholesky分解的递推公式详细推导

我们通过矩阵乘法的定义，直接推导下三角矩阵\(\boldsymbol{L}\)的元素计算公式。

设n阶对称正定矩阵\(\boldsymbol{A}=(a_{ij})_{n\times n}\)，下三角矩阵\(\boldsymbol{L}=(l_{ij})_{n\times n}\)，满足\(l_{ij}=0\)（\(i<j\)），\(l_{ii}>0\)，且\(\boldsymbol{A=LL^T}\)。

由矩阵乘法定义，对任意\(i\geq j\)（仅需计算下三角，上三角由对称性可得），有：

\[a_{ij} = \sum_{k=1}^n l_{ik} l_{jk} \]

由于\(\boldsymbol{L}\)是下三角矩阵，当\(k>i\)时\(l_{ik}=0\)，当\(k>j\)时\(l_{jk}=0\)，因此求和上限可缩减为\(\min(i,j)=j\)，即：

\[a_{ij} = \sum_{k=1}^j l_{ik} l_{jk} \]

将\(k=j\)的项单独拆分出来，得到：

\[a_{ij} = \sum_{k=1}^{j-1} l_{ik} l_{jk} + l_{ij} l_{jj} \]

我们按列递推（\(j=1,2,\dots,n\)），分两种情况得到计算公式：

（1）对角线元素\(l_{jj}\)的计算公式

当\(i=j\)时，上式变为：

\[a_{jj} = \sum_{k=1}^{j-1} l_{jk}^2 + l_{jj}^2 \]

移项开方（限定\(l_{jj}>0\)），得到：

\[\boldsymbol{l_{jj} = \sqrt{a_{jj} - \sum_{k=1}^{j-1} l_{jk}^2}, \quad j=1,2,\dots,n} \]

这就是教材中的公式(5.18)第1式。

（2）非对角线元素\(l_{ij}\)（\(i>j\)）的计算公式

当\(i>j\)时，对\(a_{ij} = \sum_{k=1}^{j-1} l_{ik} l_{jk} + l_{ij} l_{jj}\)移项，得到：

\[\boldsymbol{l_{ij} = \frac{1}{l_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} l_{ik} l_{jk} \right), \quad i=j+1,j+2,\dots,n} \]

这就是教材中的公式(5.18)第2式。

3. 平方根法求解方程组的完整步骤

完成\(\boldsymbol{A=LL^T}\)分解后，方程组\(\boldsymbol{Ax=b}\)等价于两个三角方程组：

1. 下三角方程组\(\boldsymbol{Ly=b}\)，通过前向替换求解中间变量\(\boldsymbol{y}\)；
2. 上三角方程组\(\boldsymbol{L^T x=y}\)，通过后向替换求解最终解\(\boldsymbol{x}\)。

（1）前向替换求解\(\boldsymbol{Ly=b}\)

展开下三角方程组，对\(i=1,2,\dots,n\)，推导得：

\[\boldsymbol{y_i = \frac{1}{l_{ii}} \left( b_i - \sum_{k=1}^{i-1} l_{ik} y_k \right)} \]

其中\(y_1 = b_1 / l_{11}\)。

（2）后向替换求解\(\boldsymbol{L^T x=y}\)

\(\boldsymbol{L^T}\)是上三角矩阵，展开后倒序计算，对\(i=n,n-1,\dots,1\)，推导得：

\[\boldsymbol{x_i = \frac{1}{l_{ii}} \left( y_i - \sum_{k=i+1}^n l_{ki} x_k \right)} \]

其中\(x_n = y_n / l_{nn}\)。

4. 平方根法的核心特性

1. 数值稳定性极佳
   由分解公式可得\(a_{jj} = \sum_{k=1}^j l_{jk}^2\)，因此\(l_{jk}^2 \leq a_{jj} \leq \max_{1\leq i\leq n} a_{ii}\)，即\(\boldsymbol{L}\)的所有元素的绝对值都不会超过\(\boldsymbol{A}\)对角线元素的最大值，元素数量级不会增长，舍入误差不会被放大，是天然数值稳定的算法，无需选主元。
2. 计算量极小
   平方根法的总乘除法计算量约为\(\frac{n^3}{6}\)，仅为通用LU分解法的一半，计算效率大幅提升。
3. 存储量节省
   利用对称性，仅需存储\(\boldsymbol{A}\)的下三角部分，共\(\frac{n(n+1)}{2}\)个元素，无需存储整个矩阵，内存占用仅为通用方法的一半左右，工程上可通过一维数组压缩存储。
4. 唯一局限性
   计算对角线元素时需要做开方运算，开方运算会增加少量计算耗时，且存在微小的浮点精度损失。

---

四、改进平方根法（LDLᵀ分解法）

为了避免平方根法的开方运算，我们直接使用对称矩阵的LDLᵀ分解，得到改进平方根法，这是工程中求解对称正定方程组的首选算法。

1. 核心思想

直接基于\(\boldsymbol{A=LDL^T}\)分解求解方程组，其中\(\boldsymbol{L}\)为单位下三角矩阵，\(\boldsymbol{D}\)为对角矩阵，全程无需开方运算，计算量与平方根法相当，且保留了数值稳定、计算高效的全部优势。

2. LDLᵀ分解的递推公式推导

设\(\boldsymbol{A=LDL^T}\)，\(\boldsymbol{L}=(l_{ij})\)为单位下三角矩阵（\(l_{ii}=1\)，\(l_{ij}=0\)，\(i<j\)），\(\boldsymbol{D}=\text{diag}(d_1,d_2,\dots,d_n)\)为对角矩阵。

由矩阵乘法定义，对\(i\geq j\)，有：

\[a_{ij} = \sum_{k=1}^n (\boldsymbol{LD})_{ik} (\boldsymbol{L^T})_{kj} = \sum_{k=1}^j l_{ik} d_k l_{jk} \]

将\(k=j\)的项单独拆分（\(l_{jj}=1\)），得到：

\[a_{ij} = \sum_{k=1}^{j-1} l_{ik} d_k l_{jk} + l_{ij} d_j \]

我们按行递推（\(i=1,2,\dots,n\)），得到计算公式：

（1）非对角线元素\(l_{ij}\)（\(j<i\)）的计算公式

对\(j=1,2,\dots,i-1\)，移项得：

\[l_{ij} = \frac{1}{d_j} \left( a_{ij} - \sum_{k=1}^{j-1} l_{ik} d_k l_{jk} \right) \]

（2）对角矩阵元素\(d_i\)的计算公式

当\(j=i\)时，\(a_{ii} = \sum_{k=1}^{i-1} l_{ik}^2 d_k + d_i\)，移项得：

\[d_i = a_{ii} - \sum_{k=1}^{i-1} l_{ik}^2 d_k \]

3. 计算优化：避免重复计算

上述公式中，\(l_{ik} d_k\)会被重复计算，为了提升效率，引入中间变量\(t_{ij}=l_{ij} d_j\)，将公式改写为按行计算的优化形式：

1. 初始化：\(d_1 = a_{11}\)
2. 对\(i=2,3,\dots,n\)，依次执行：
   ① 计算中间量：\(\boldsymbol{t_{ij} = a_{ij} - \sum_{k=1}^{j-1} t_{ik} l_{jk}, \quad j=1,2,\dots,i-1}\)
   ② 计算\(\boldsymbol{L}\)的元素：\(\boldsymbol{l_{ij} = \frac{t_{ij}}{d_j}, \quad j=1,2,\dots,i-1}\)
   ③ 计算\(\boldsymbol{D}\)的对角元：\(\boldsymbol{d_i = a_{ii} - \sum_{k=1}^{i-1} t_{ik} l_{ik}}\)

该优化形式避免了重复计算，进一步提升了计算效率，是工程实现的标准形式。

4. 改进平方根法求解方程组的完整步骤

完成\(\boldsymbol{A=LDL^T}\)分解后，方程组\(\boldsymbol{Ax=b}\)等价于：

\[\boldsymbol{LDL^T x = b} \]

拆分为两个三角方程组分步求解：

1. 单位下三角方程组\(\boldsymbol{Ly=b}\)，前向替换求解\(\boldsymbol{y}\)；
2. 对角+上三角方程组\(\boldsymbol{DL^T x=y}\)，先做对角缩放，再后向替换求解\(\boldsymbol{x}\)。

（1）前向替换求解\(\boldsymbol{Ly=b}\)

\(\boldsymbol{L}\)是单位下三角矩阵，对角线元素为1，因此公式为：

\[\boldsymbol{ \begin{cases} y_1 = b_1 \\ y_i = b_i - \sum_{k=1}^{i-1} l_{ik} y_k, \quad i=2,3,\dots,n \end{cases} }\]

（2）后向替换求解\(\boldsymbol{DL^T x=y}\)

先令\(\boldsymbol{z=D^{-1}y}\)，即\(z_i = y_i / d_i\)，再求解\(\boldsymbol{L^T x=z}\)，最终公式为：

\[\boldsymbol{ \begin{cases} x_n = y_n / d_n \\ x_i = \frac{y_i}{d_i} - \sum_{k=i+1}^n l_{ki} x_k, \quad i=n-1,n-2,\dots,1 \end{cases} }\]

5. 改进平方根法的核心优势

1. 全程无开方运算，避免了开方带来的计算耗时和精度损失；
2. 计算量与平方根法相当，约为\(\frac{n^3}{6}\)，仅为通用LU分解的一半；
3. 保留了平方根法数值稳定、无需选主元、存储节省的全部优势；
4. 公式规整，易于编程实现，是工业级软件求解对称正定方程组的标准算法。

---

五、算例验证

我们用一个典型的对称正定矩阵，分别用平方根法和改进平方根法求解，验证算法的正确性。

例题：求解对称正定方程组

\[\begin{pmatrix} 4 & 2 & -2 \\ 2 & 2 & -3 \\ -2 & -3 & 14 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 10 \\ 5 \\ 4 \end{pmatrix} \]

方法1：平方根法（Cholesky分解）

步骤1：计算\(\boldsymbol{L}\)的元素

1. \(j=1\)：
   \(l_{11}=\sqrt{a_{11}}=\sqrt{4}=2\)；
   \(l_{21}=a_{21}/l_{11}=2/2=1\)；
   \(l_{31}=a_{31}/l_{11}=-2/2=-1\)；
2. \(j=2\)：
   \(l_{22}=\sqrt{a_{22}-l_{21}^2}=\sqrt{2-1}=1\)；
   \(l_{32}=(a_{32}-l_{31}l_{21})/l_{22}=(-3 - (-1)\times1)/1=-2\)；
3. \(j=3\)：
   \(l_{33}=\sqrt{a_{33}-l_{31}^2-l_{32}^2}=\sqrt{14-1-4}=3\)；

最终得到：

\[\boldsymbol{L}=\begin{pmatrix} 2 & 0 & 0 \\ 1 & 1 & 0 \\ -1 & -2 & 3 \end{pmatrix} \]

步骤2：求解\(\boldsymbol{Ly=b}\)

\(\boldsymbol{b}=(10,5,4)^T\)，前代得：
\(y_1=10/2=5\)；
\(y_2=(5-1\times5)/1=0\)；
\(y_3=(4 - (-1)\times5 - (-2)\times0)/3=3\)；
\(\boldsymbol{y}=(5,0,3)^T\)。

步骤3：求解\(\boldsymbol{L^T x=y}\)

回代得：
\(x_3=3/3=1\)；
\(x_2=(0 - (-2)\times1)/1=2\)；
\(x_1=(5 - 1\times2 - (-1)\times1)/2=2\)；

最终解：\(\boldsymbol{x}=(2,2,1)^T\)，代入原方程验证成立。

方法2：改进平方根法（LDLᵀ分解）

步骤1：计算\(\boldsymbol{L}\)和\(\boldsymbol{D}\)

1. \(i=1\)：\(d_1=a_{11}=4\)；
2. \(i=2\)：
   \(t_{21}=a_{21}=2\)；
   \(l_{21}=t_{21}/d_1=2/4=0.5\)；
   \(d_2=a_{22}-t_{21}l_{21}=2-2\times0.5=1\)；
3. \(i=3\)：
   \(t_{31}=a_{31}=-2\)；\(l_{31}=t_{31}/d_1=-2/4=-0.5\)；
   \(t_{32}=a_{32}-t_{31}l_{21}=-3 - (-2)\times0.5=-2\)；
   \(l_{32}=t_{32}/d_2=-2/1=-2\)；
   \(d_3=a_{33}-t_{31}l_{31}-t_{32}l_{32}=14 - (-2)\times(-0.5) - (-2)\times(-2)=9\)；

最终得到：

\[\boldsymbol{L}=\begin{pmatrix} 1 & 0 & 0 \\ 0.5 & 1 & 0 \\ -0.5 & -2 & 1 \end{pmatrix}, \quad \boldsymbol{D}=\begin{pmatrix} 4 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 9 \end{pmatrix} \]

步骤2：求解\(\boldsymbol{Ly=b}\)

前代得：\(y_1=10\)，\(y_2=0\)，\(y_3=9\)。

步骤3：求解\(\boldsymbol{DL^T x=y}\)

回代得：\(x_3=9/9=1\)，\(x_2=0/1 - (-2)\times1=2\)，\(x_1=10/4 - 0.5\times2 - (-0.5)\times1=2\)，与平方根法结果完全一致。

---

六、核心知识点归纳总结表

方法	通用LU分解	平方根法（Cholesky分解\(\boldsymbol{LL^T}\)）	改进平方根法（\(\boldsymbol{LDL^T}\)分解）
适用矩阵	顺序主子式非零的方阵	对称正定矩阵	对称正定/对称非奇异矩阵
核心分解式	\(\boldsymbol{A=LU}\)	\(\boldsymbol{A=LL^T}\)	\(\boldsymbol{A=LDL^T}\)
矩阵特性	\(\boldsymbol{L}\)单位下三角，\(\boldsymbol{U}\)上三角	\(\boldsymbol{L}\)下三角，对角元为正	\(\boldsymbol{L}\)单位下三角，\(\boldsymbol{D}\)对角矩阵
计算量	约\(\frac{n^3}{3}\)次乘除法	约\(\frac{n^3}{6}\)次乘除法，含n次开方	约\(\frac{n^3}{6}\)次乘除法，无开方运算
存储量	需存储整个n阶矩阵，\(n^2\)个元素	仅需存储下三角，\(\frac{n(n+1)}{2}\)个元素	仅需存储下三角，\(\frac{n(n+1)}{2}\)个元素
数值稳定性	不稳定，需选主元	天然稳定，无需选主元	天然稳定，无需选主元
核心优势	适用范围广，通用型强	计算效率高，存储节省，数值稳定	无开方运算，效率最高，工程实用性最强
局限性	计算和存储冗余，对称矩阵未利用对称性	需开方运算，有微小精度损失	仅适用于对称矩阵，通用性弱于LU分解
核心应用场景	通用非对称线性方程组求解	小规模对称正定方程组求解	有限元、结构力学等大规模对称正定方程组工程求解

---

追赶法（Thomas算法）完整讲解

追赶法（也叫托马斯Thomas算法）是专门针对三对角线性方程组的高效定制化解法，是高斯消去法与三角分解法在三对角稀疏矩阵上的极致优化，广泛应用于常微分方程边值问题求解、热传导方程数值计算、三次样条插值、船体数学放样等工程与科学计算场景。本次讲解将从理论基础、矩阵分解推导、算法流程、稳定性证明、数值特性、算例验证全链条展开，最终完成核心知识点的系统归纳。

---

一、三对角方程组与核心理论前提

1. 三对角方程组的标准形式

我们求解的线性方程组\(\boldsymbol{Ax=f}\)，其系数矩阵\(\boldsymbol{A}\)为三对角矩阵：仅主对角线、主对角线下一条次对角线、主对角线上一条次对角线有非零元素，其余元素全为0，标准形式为：

\[\begin{pmatrix} b_1 & c_1 & & & \\ a_2 & b_2 & c_2 & & \\ & \ddots & \ddots & \ddots & \\ & & a_{n-1} & b_{n-1} & c_{n-1} \\ & & & a_n & b_n \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_{n-1} \\ x_n \end{pmatrix} = \begin{pmatrix} f_1 \\ f_2 \\ \vdots \\ f_{n-1} \\ f_n \end{pmatrix} \]

简记为\(\boldsymbol{Ax=f}\)，满足：当\(|i-j|>1\)时，\(a_{ij}=0\)。

2. 方程组有唯一解的前提条件

为保证方程组有唯一解、算法数值稳定，要求三对角矩阵\(\boldsymbol{A}\)满足对角占优条件：
① \(|b_1| > |c_1| > 0\)；
② \(|b_i| \geq |a_i| + |c_i|\)，且\(a_i,c_i \neq 0\)，\(i=2,3,\dots,n-1\)；
③ \(|b_n| > |a_n| > 0\)。

为严谨证明解的存在唯一性，我们先引入两个核心定义：

定义1：对角占优矩阵

设n阶矩阵\(\boldsymbol{A}=(a_{ij})_{n\times n}\)：

1. 严格对角占优矩阵：若对所有\(i=1,2,\dots,n\)，都有
   \[|a_{ii}| > \sum_{\substack{j=1 \\ j\neq i}}^n |a_{ij}| \]
   即主对角线元素的绝对值，严格大于该行其余所有元素绝对值之和。
2. 弱对角占优矩阵：若对所有\(i=1,2,\dots,n\)，都有
   \[|a_{ii}| \geq \sum_{\substack{j=1 \\ j\neq i}}^n |a_{ij}| \]
   且至少存在一个\(i\)，使得上述不等式严格成立。

定义2：可约与不可约矩阵

设n阶矩阵\(\boldsymbol{A}(n\geq2)\)，若存在置换矩阵\(\boldsymbol{P}\)，使得

\[\boldsymbol{P^TAP} = \begin{pmatrix} \boldsymbol{A_{11}} & \boldsymbol{A_{12}} \\ \boldsymbol{0} & \boldsymbol{A_{22}} \end{pmatrix} \]

其中\(\boldsymbol{A_{11}}\)为r阶方阵，\(\boldsymbol{A_{22}}\)为n-r阶方阵（\(1\leq r <n\)），则称\(\boldsymbol{A}\)为可约矩阵；否则称\(\boldsymbol{A}\)为不可约矩阵。
通俗理解：可约矩阵可通过行列重排拆分为两个低阶方程组求解；而所有\(a_i,c_i\neq0\)的三对角矩阵是典型的不可约矩阵，无法拆分。

定理5.11（对角占优定理）

若\(\boldsymbol{A}\)为严格对角占优矩阵，或不可约的弱对角占优矩阵，则\(\boldsymbol{A}\)为非奇异矩阵（可逆）。

详细证明（反证法，针对严格对角占优矩阵）

假设\(\boldsymbol{A}\)是奇异矩阵，则齐次方程组\(\boldsymbol{Ax=0}\)存在非零解\(\boldsymbol{x}=(x_1,x_2,\dots,x_n)^T\)。
设\(|x_k| = \max_{1\leq i\leq n} |x_i| \neq 0\)，即\(x_k\)是解向量中绝对值最大的元素。
对齐次方程组的第k个方程：

\[\sum_{j=1}^n a_{kj}x_j = 0 \]

移项得：

\[a_{kk}x_k = -\sum_{\substack{j=1 \\ j\neq k}}^n a_{kj}x_j \]

两边取绝对值，由三角不等式：

\[|a_{kk}| |x_k| \leq \sum_{\substack{j=1 \\ j\neq k}}^n |a_{kj}| |x_j| \]

由于\(|x_j| \leq |x_k|\)对所有j成立，代入得：

\[|a_{kk}| |x_k| \leq |x_k| \sum_{\substack{j=1 \\ j\neq k}}^n |a_{kj}| \]

两边约去非零的\(|x_k|\)，得到：

\[|a_{kk}| \leq \sum_{\substack{j=1 \\ j\neq k}}^n |a_{kj}| \]

这与\(\boldsymbol{A}\)是严格对角占优矩阵的定义矛盾，因此假设不成立，\(\boldsymbol{A}\)非奇异。

三对角方程组解的唯一性结论

满足对角占优条件①②③的三对角矩阵\(\boldsymbol{A}\)，是不可约的弱对角占优矩阵，因此非奇异，对应的三对角方程组\(\boldsymbol{Ax=f}\)存在唯一解，为追赶法提供了理论基础。

---

二、三对角矩阵的LU分解与追赶法公式推导

我们利用三对角矩阵的稀疏特性，对\(\boldsymbol{A}\)做定制化的LU分解，避免通用LU分解的冗余计算，推导追赶法的核心公式。

1. 定制化LU分解形式

针对三对角矩阵的结构，我们将\(\boldsymbol{A}\)分解为：

\[\boldsymbol{A=LU} \]

其中：

• \(\boldsymbol{L}\)是下双对角矩阵（仅主对角线和下次对角线有非零元素）：
  \[\boldsymbol{L} = \begin{pmatrix} \alpha_1 & & & & \\ a_2 & \alpha_2 & & & \\ & \ddots & \ddots & & \\ & & a_{n-1} & \alpha_{n-1} & \\ & & & a_n & \alpha_n \end{pmatrix}\]

• \(\boldsymbol{U}\)是单位上双对角矩阵（主对角线全为1，仅上次对角线有非零元素）：
  \[\boldsymbol{U} = \begin{pmatrix} 1 & \beta_1 & & & \\ & 1 & \beta_2 & & \\ & & \ddots & \ddots & \\ & & & 1 & \beta_{n-1} \\ & & & & 1 \end{pmatrix}\]

这种分解形式完全匹配三对角矩阵的稀疏性，分解后\(\boldsymbol{L}\)和\(\boldsymbol{U}\)也仅保留三条对角线的非零元素，无冗余计算。

2. 分解系数的递推公式推导

通过矩阵乘法\(\boldsymbol{A=LU}\)，对比等式两边对应位置的元素，推导\(\alpha_i\)和\(\beta_i\)的计算公式。

（1）第1行元素对比

\(\boldsymbol{A}\)的第1行：\([b_1, c_1, 0, \dots, 0]\)
\(\boldsymbol{LU}\)的第1行：\([\alpha_1 \times 1, \alpha_1 \times \beta_1, 0, \dots, 0]\)
对应相等得：

\[b_1 = \alpha_1, \quad c_1 = \alpha_1 \beta_1 \]

因此得到初始值：

\[\boldsymbol{\alpha_1 = b_1, \quad \beta_1 = \frac{c_1}{b_1}} \]

（2）第i行元素对比（\(2\leq i\leq n-1\)）

\(\boldsymbol{A}\)的第i行非零元素：\(a_i, b_i, c_i\)
\(\boldsymbol{LU}\)的第i行展开后：

• 第i-1列元素：\(a_i \times 1 = a_i\)，与\(\boldsymbol{A}\)的\(a_i\)完全匹配，无需额外计算；
• 第i列元素：\(a_i \times \beta_{i-1} + \alpha_i \times 1 = b_i\)；
• 第i+1列元素：\(\alpha_i \times \beta_i = c_i\)。

整理得到递推公式：

\[\boldsymbol{\alpha_i = b_i - a_i \beta_{i-1}, \quad \beta_i = \frac{c_i}{\alpha_i}, \quad i=2,3,\dots,n-1} \]

（3）第n行元素对比

\(\boldsymbol{A}\)的第n行非零元素：\(a_n, b_n\)
\(\boldsymbol{LU}\)的第n行第n列元素：\(a_n \times \beta_{n-1} + \alpha_n = b_n\)
整理得到：

\[\boldsymbol{\alpha_n = b_n - a_n \beta_{n-1}} \]

3. 分解可行性与稳定性证明（归纳法）

我们证明：满足对角占优条件的三对角矩阵，所有\(\alpha_i \neq 0\)，分解不会出现除以零的情况，且\(0<|\beta_i|<1\)，保证数值稳定。

• 基例（i=1）：由条件①\(|b_1|>|c_1|>0\)，因此\(\alpha_1=b_1\neq0\)，且\(|\beta_1|=|c_1/b_1|<1\)，同时\(|\alpha_1|=|b_1|>|c_1|>0\)。
• 归纳假设：假设对i-1，有\(|\alpha_{i-1}|>|c_{i-1}|>0\)，且\(0<|\beta_{i-1}|<1\)。
• 归纳步骤（i）：
  由三角不等式和对角占优条件②：
  \[|\alpha_i| = |b_i - a_i \beta_{i-1}| \geq |b_i| - |a_i| |\beta_{i-1}| \]
  由归纳假设\(|\beta_{i-1}|<1\)，因此：
  \[|\alpha_i| > |b_i| - |a_i| \geq |c_i| > 0 \]
  因此\(\alpha_i \neq 0\)，且\(|\beta_i|=|c_i/\alpha_i|<1\)，归纳成立。

最终得到定理5.12的核心结论：

1. \(0<|\beta_i|<1\)，\(i=1,2,\dots,n-1\)；
2. \(0<|c_i| \leq |b_i|-|a_i| < |\alpha_i| < |b_i|+|a_i|\)，\(i=2,\dots,n-1\)；\(0<|b_n|-|a_n|<|\alpha_n|<|b_n|+|a_n|\)。

该结论说明：分解过程中的中间量\(\alpha_i,\beta_i\)均有界，不会出现数值爆炸，因此追赶法是数值稳定的算法，无需选主元。

---

三、追赶法求解方程组的完整算法流程

完成\(\boldsymbol{A=LU}\)分解后，方程组\(\boldsymbol{Ax=f}\)等价于两个双对角三角方程组：

1. 下双对角方程组\(\boldsymbol{Ly=f}\)，通过前向替换（追的过程）求解中间变量\(\boldsymbol{y}\)；
2. 单位上双对角方程组\(\boldsymbol{Ux=y}\)，通过后向替换（赶的过程）求解最终解\(\boldsymbol{x}\)。

完整算法分为3个核心步骤：

步骤1：计算分解系数\(\{\beta_i\}\)（追的过程第一部分）

从前往后依次计算：

1. 初始值：\(\beta_1 = \frac{c_1}{b_1}\)；
2. 递推计算：\(\beta_i = \frac{c_i}{b_i - a_i \beta_{i-1}}\)，\(i=2,3,\dots,n-1\)。

步骤2：前向替换求解\(\boldsymbol{Ly=f}\)（追的过程第二部分）

从前往后依次计算\(\boldsymbol{y}\)的元素：

1. 初始值：\(y_1 = \frac{f_1}{b_1}\)；
2. 递推计算：\(y_i = \frac{f_i - a_i y_{i-1}}{b_i - a_i \beta_{i-1}}\)，\(i=2,3,\dots,n\)。

步骤3：后向替换求解\(\boldsymbol{Ux=y}\)（赶的过程）

从后往前依次计算\(\boldsymbol{x}\)的元素：

1. 初始值：\(x_n = y_n\)；
2. 递推计算：\(x_i = y_i - \beta_i x_{i+1}\)，\(i=n-1,n-2,\dots,1\)。

算法名称由来：步骤1和步骤2是从第1个元素到第n个元素的前向迭代，称为“追”的过程；步骤3是从第n个元素到第1个元素的后向迭代，称为“赶”的过程，因此该算法被形象地称为“追赶法”。

---

四、追赶法的核心数值特性与工程优势

1. 计算量极小
   追赶法的总乘除法计算量仅为\(\boldsymbol{5n-4}\)次，远低于通用LU分解的\(\frac{n^3}{3}\)次和高斯消去法的\(\frac{n^3}{3}\)次。当n很大时（如n=10000），追赶法的计算量仅为通用方法的几十万分之一，效率提升极为显著。
   对于同系数矩阵、多右端项的方程组，每新增一个右端项，仅需增加\(3n-2\)次乘除运算，批量计算效率极高。

2. 存储量极省
   无需存储整个n阶矩阵，仅需3个一维数组分别存储三对角线元素\(\{a_i\},\{b_i\},\{c_i\}\)，再用2个一维数组存储中间量\(\{\beta_i\}\)和\(\{y_i\}\)即可，内存占用仅为\(O(n)\)，适合求解超大规模的三对角方程组（如偏微分方程离散后的百万阶方程组）。

3. 数值稳定性极佳
   由定理5.12的结论，分解过程中的中间量\(\beta_i\)满足\(|\beta_i|<1\)，所有中间量均有界，舍入误差不会被放大，是天然数值稳定的算法，无需选主元，进一步简化了计算流程。

4. 编程实现极简
   算法仅为简单的循环递推，无复杂的矩阵操作，代码量极少，易于工程实现和调试，是工业级数值软件中求解三对角方程组的标准算法。

---

五、算例验证

我们用一个典型的三对角方程组，完整演示追赶法的计算流程。

例题：用追赶法求解三对角方程组

\[\begin{pmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \\ 1 \end{pmatrix} \]

解：
首先确定三对角线元素：
\(a_2=-1, a_3=-1\)；\(b_1=2, b_2=2, b_3=2\)；\(c_1=-1, c_2=-1\)；右端项\(f=(1,0,1)^T\)。
验证对角占优条件：\(|b_1|=2>|-1|=|c_1|\)，\(|b_2|=2=|-1|+|-1|\)，\(|b_3|=2>|-1|=|a_3|\)，满足条件，可使用追赶法。

步骤1：计算\(\{\beta_i\}\)

• \(\beta_1 = c_1/b_1 = -1/2 = -0.5\)；
• \(\beta_2 = c_2/(b_2 - a_2 \beta_1) = -1/(2 - (-1)\times(-0.5)) = -2/3\)；

步骤2：前向替换求解\(\boldsymbol{y}\)

• \(y_1 = f_1/b_1 = 1/2 = 0.5\)；
• \(y_2 = (f_2 - a_2 y_1)/(b_2 - a_2 \beta_1) = (0 - (-1)\times0.5)/1.5 = 1/3\)；
• \(y_3 = (f_3 - a_3 y_2)/(b_3 - a_3 \beta_2) = (1 - (-1)\times(1/3))/(2 - (-1)\times(-2/3)) = 1\)；

步骤3：后向替换求解\(\boldsymbol{x}\)

• \(x_3 = y_3 = 1\)；
• \(x_2 = y_2 - \beta_2 x_3 = 1/3 - (-2/3)\times1 = 1\)；
• \(x_1 = y_1 - \beta_1 x_2 = 0.5 - (-0.5)\times1 = 1\)；

最终解：\(\boldsymbol{x}=(1,1,1)^T\)，代入原方程组验证，等式完全成立。

---

六、核心知识点归纳总结表

方法	通用LU分解	高斯消去法	追赶法（Thomas算法）
适用矩阵	顺序主子式非零的n阶方阵	非奇异n阶方阵	对角占优的三对角矩阵
核心形式	\(\boldsymbol{A=LU}\)，L单位下三角，U上三角	增广矩阵行变换消元	定制化\(\boldsymbol{A=LU}\)，L下双对角，U单位上双对角
计算量	约\(\frac{n^3}{3}\)次乘除法	约\(\frac{n^3}{3}\)次乘除法	仅\(5n-4\)次乘除法，线性时间复杂度
存储量	\(O(n^2)\)，需存储整个矩阵	\(O(n^2)\)，需存储增广矩阵	\(O(n)\)，仅需3个一维数组存储三对角线
数值稳定性	不稳定，需选主元	不稳定，需选主元	天然稳定，无需选主元，中间量有界
编程复杂度	中等，需双重循环处理矩阵	中等，需行变换和消元操作	极低，仅需简单的单循环递推
核心应用场景	通用非对称线性方程组	小规模通用线性方程组	三次样条插值、常微分方程边值问题、热传导等大规模三对角方程组求解