Lesson 63: LU分解求解器(带部分主元)
练习任务
难度：难
实现一个带部分主元（partial pivoting）的 LU 分解线性方程组求解器，求解固定的 3×3 系统：
A = [[2, 1, 1],
[4, 3, 3], b = [5, 11, 29] 解: x = [2, -2, 3]
[8, 7, 9]]你需要完成七个核心函数：
print_matrix(label, m)— 格式化打印 N×N 矩阵print_vector(label, v)— 格式化打印 N 维向量pivot(A, p, k)— 部分主元选择与行交换lu_decomp(A, L, U, p)— PA = LU 分解forward_subst(L, b, p, y)— 前向替换 Ly = Pbback_subst(U, y, x)— 后向替换 Ux = ysolve(A, b, x)+main()— 主流程编排与验证
程序输出每一步的中间结果（L 矩阵、U 矩阵、y 向量、x 解向量），最后验证 Ax = b。
预期输出（核心段）：
=== LU Decomposition Solver with Partial Pivoting ===
Step 1: PA = LU decomposition
Pivot: swap row 0 <-> row 2
Pivot: swap row 1 <-> row 2
Permutation vector p: [2 0 1]
L matrix:
1.0000 0.0000 0.0000
0.2500 1.0000 0.0000
0.5000 0.6667 1.0000
U matrix:
8.0000 7.0000 9.0000
0.0000 -0.7500 -1.2500
0.0000 0.0000 -0.6667
Step 2: Forward substitution Ly = Pb
y: [ 29.0000 -2.2500 -2.0000]
Step 3: Backward substitution Ux = y
x (solution): [ 2.0000 -2.0000 3.0000]
=== Verification: Ax ===
Ax: [ 5.0000 11.0000 29.0000]
b: [ 5.0000 11.0000 29.0000]提示：LU 分解的核心在于将"求解 Ax=b"拆成三步——先分解 PA=LU（O(n³)），再解两个三角系统 Ly=Pb（前向）和 Ux=y（后向）（各 O(n²)）。如果 A 不变而 b 变化多次，只需分解只需做一次。
核心知识点
- LU 分解的矩阵视角 — A 被分解为单位下三角 L 和上三角 U 的乘积；本质是高斯消元的矩阵形式存储
- 部分主元（Partial Pivoting） — 每步选择当前列中绝对值最大的元素作为主元，抑制浮点误差放大，保证 |乘数| ≤ 1
- 前向替换 Ly = Pb — 利用 L 是单位下三角（对角线全 1），从第一行开始直接逐行计算 y，无需除法
- 后向替换 Ux = y — 利用 U 是上三角，从最后一行向前逐行求解，必须除以对角元 U[i][i]
- 排列向量 p 的作用 — p 记录行交换历史：PA = LU 中的 P 由 p 隐式表示，前向替换时用它重排 b
- 浮点精度与数值稳定性 — 部分主元保证所有乘数 l_ik ∈ [-1,1]，阻止舍入误差在后续步骤中被放大
- 三角求解 O(n²) 的效率优势 — 一般线性系统 O(n³)，三角系统只需 O(n²)；多次求解右端项时优势显著
代码框架
c
#include <math.h>
#include <stdio.h>
#define N 3
/* ─── TODO 1: 打印 N×N 矩阵 ─── */
void print_matrix(const char *label, double m[N][N]) {
// ① printf("%s:\n", label)
// ② for i = 0..N-1：
// for j = 0..N-1：printf("%8.4f ", m[i][j])
// printf("\n")
}
/* ─── TODO 2: 打印 N 维向量 ─── */
void print_vector(const char *label, double v[N]) {
// printf("%s: [", label)
// for i = 0..N-1：printf("%8.4f%c", v[i], i<N-1 ? ' ' : ?)
// printf("]\n")
}
/* ─── TODO 3: 部分主元 ─── */
void pivot(double A[N][N], int p[N], int k) {
// ① 在 [k, N-1] 行中找 |A[i][k]| 最大的行号 max_row
// 用 fabs()（需 #include <math.h> 和链接 -lm）
// ② 如果 max_row != k：
// 交换 A 的第 k 行和第 max_row 行（所有 N 列）
// 交换 p[k] 和 p[max_row]
// printf(" Pivot: swap row %d <-> row %d\n", k, max_row)
}
/* ─── TODO 4: PA = LU 分解 ─── */
void lu_decomp(double A[N][N], double L[N][N], double U[N][N], int p[N]) {
// ① 初始化排列：p[i] = i for i = 0..N-1
//
// ② 对 k = 0..N-1：
// a. pivot(A, p, k)
// b. 对 i = k+1..N-1：
// A[i][k] /= A[k][k] // 乘数 l_ik
// 对 j = k+1..N-1：
// A[i][j] -= A[i][k] * A[k][j] // Schur 补
//
// ③ 提取 L 和 U：
// if i > j: L[i][j]=A[i][j], U[i][j]=0
// if i == j: L[i][j]=1, U[i][j]=A[i][j]
// if i < j: L[i][j]=0, U[i][j]=A[i][j]
}
/* ───── TODO 5: 前向替换 Ly = Pb ─── */
void forward_subst(double L[N][N], double b[N], int p[N], double y[N]) {
// ① 用 p 重排 b：Pb[i] = b[p[i]]
// ② y[i] = Pb[i] - Σ(j=0 to i-1) L[i][j] * y[j]
// （L[i][i] = 1，无需除法）
}
/* ─── TODO 6: 后向替换 Ux = y ─── */
void back_subst(double U[N][N], double y[N], double x[N]) {
// 从 i = N-1 到 0：
// x[i] = (y[i] - Σ(j=i+1 to N-1) U[i][j] * x[j]) / U[i][i]
// ⚠ 最后除以 U[i][i] 不能忘！
}
/* ─── TODO 7: 主流程 ─── */
void solve(double A_orig[N][N], double b[N], double x[N]) {
// ① 复制 A_orig 到 local A（分解会修改 A）
// ② 打印标题
// ③ 调用 lu_decomp → 打印 L, U, p
// ④ 调用 forward_subst → 打印 y
// ⑤ 调用 back_subst → 打印 x
// ⑥ 验证：用 A_orig 计算 Ax，与 b 比较
}
int main(void) {
// double A[3][3] = {{2,1,1},{4,3,3},{8,7,9}}
// double b[3] = {5, 11, 29}
// double x[3]
// solve(A, b, x)
// return 0
}阅读骨架后，尝试自己填充 // ① 到 // ⑥ 标记的部分。关键挑战：lu_decomp 中"先 pivot 再消去"的嵌套顺序是否正确？提取 L/U 时对角线如何分配？前向替换中排列向量 p 如何使用？后向替换最后除以 U[i][i] 是否遗漏？
TIP
先在纸上用 A = [[2,1,1],[4,3,3],[8,7,9]] 手动完成 PA=LU 分解全过程。代入三重消元——k=0 选主元交换 row0↔row2 → 消去 row1,row2 → k=1 再选主元交换 row1↔row2 → 消去 row2。验算 L×U 是否等于 PA。
深度讲解
1. LU 分解——高斯消元的矩阵形式
1.1 从高斯消元到 LU
高斯消元大家都会：用某行的倍数去消另一行。但每次消元操作实际上是一个矩阵乘法。把所有消元操作的逆累积起来，就得到了 L。
高斯消元的过程对应三个初等矩阵：
消去第 1 列（k=0）:
E₁ = [[ 1, 0, 0], E₁⁻¹ = [[ 1, 0, 0],
[-2, 1, 0], [ 2, 1, 0],
[-4, 0, 1]] [ 4, 0, 1]]
消去第 2 列（k=1）:
E₂ = [[ 1, 0, 0], E₂⁻¹ = [[ 1, 0, 0],
[ 0, 1, 0], [ 0, 1, 0],
[ 0, -3, 1]] [ 0, 3, 1]]
则: E₂·E₁·A = U → A = E₁⁻¹·E₂⁻¹·U = L·U
L = E₁⁻¹·E₂⁻¹ = [[1,0,0], ← 消元乘数 l_ik 恰好填入 L！
[2,1,0],
[4,3,1]]c
/* LU 分解本质：记录高斯消元中每一步的乘数 */IMPORTANT
LU 分解不是新算法——它就是高斯消元。区别在于：高斯消元做完就忘（只留下 U），LU 把每一步的乘数保存在 L 中，使得后续求解新右端项 b 时可以直接复用 L 和 U。
1.2 In-Place 存储与分解过程
A 在分解后同时存储 L 和 U（L 对角线上的 1 不存储）：
A (分解后): 提取规则:
┌─────────────────────┐ i > j → L[i][j]=A[i][j], U[i][j]=0
│ u₀₀ u₀₁ u₀₂ │ i == j → L[i][j]=1, U[i][j]=A[i][j]
│ l₁₀ u₁₁ u₁₂ │ i < j → L[i][j]=0, U[i][j]=A[i][j]
│ l₂₀ l₂₁ u₂₂ │
└─────────────────────┘以本题 A = [[2,1,1],[4,3,3],[8,7,9]] 的完整分解过程：
k=0: 列 0 选主元 — |2|=2, |4|=4, |8|=8 → max=8 at row2, swap row0↔row2
A=[[8,7,9],[4,3,3],[2,1,1]], p=[2,1,0]
消去: l₁₀=4/8=0.5, l₂₀=2/8=0.25
A=[[8,7,9],[0.5,-0.5,-1.5],[0.25,-0.75,-1.25]]
k=1: 列 1 选主元 — |-0.5|=0.5, |-0.75|=0.75 → max=0.75 at row2, swap row1↔row2
A=[[8,7,9],[0.25,-0.75,-1.25],[0.5,-0.5,-1.5]], p=[2,0,1]
消去: l₂₁=-0.5/(-0.75)=0.6667
A=[[8,7,9],[0.25,-0.75,-1.25],[0.5,0.6667,-0.6667]]
提取:
L=[[1,0,0],[0.25,1,0],[0.5,0.6667,1]], U=[[8,7,9],[0,-0.75,-1.25],[0,0,-0.6667]]
验证: L×U = [[8,7,9],[2,1,1],[4,3,3]] = PA ✓WARNING
提取 L/U 时最常见的错误是把对角线元素同时放入 L 和 U。记住：L 的对角线固定为 1（单位下三角），U 的对角线 = A[i][i]。两者必须严格区分。
2. 部分主元——数值稳定性的基石
2.1 没有主元会怎样？
不选主元时，小元素出现在对角线上会导致乘数巨大，浮点误差被灾难性放大：
示例: ε = 1e-16（接近机器精度的微小正数）
[ε 1][x₁] = [1]
[1 1][x₂] [2]
不选主元: l₂₁=1/ε≈1e16 → (1-1e16)x₂≈2-1e16
双精度下: 1-1e16≈-1e16, 2-1e16≈-1e16 → x₂≈1.0
x₁=(1-1)/1e-16=0 ← 完全错误！（真实解 x₁=x₂≈1）
选主元: swap → l₂₁=ε≈1e-16 → x₂≈1.0, x₁≈1.0 ← 正确！
这就是"灾难性抵消"——相近大数相减，有效数字全部丢失。NOTE
部分主元保证所有乘数 $|l_{ik}| \le 1$，从而抑制误差在后续步骤中的传播。这是 LU 分解成为工业标准的关键——不是因为它更快，而是因为它稳定。
2.2 部分主元的实现细节
c
void pivot(double A[N][N], int p[N], int k) {
int max_row = k;
double max_val = fabs(A[k][k]); // ⚠ 初始化为 |A[k][k]|，不是 0！
/* 在 [k, N-1] 行中找绝对值最大的 */
for (int i = k + 1; i < N; i++) {
double abs_val = fabs(A[i][k]);
if (abs_val > max_val) {
max_val = abs_val;
max_row = i;
}
}
if (max_row != k) {
/* 交换 A 的整行——所有 N 列都需要交换！ */
for (int j = 0; j < N; j++) {
double tmp = A[k][j];
A[k][j] = A[max_row][j];
A[max_row][j] = tmp;
}
/* 同步更新排列向量 */
int tmp_p = p[k];
p[k] = p[max_row];
p[max_row] = tmp_p;
printf(" Pivot: swap row %d <-> row %d\n", k, max_row);
}
}CAUTION
两个关键陷阱：（1）max_val 必须初始化为 fabs(A[k][k])，如果初始化成 0，当所有待选元素为负时 max_val 永远不会被更新；（2）交换时必须是整行（所有 N 列都交换），不能只交换第 k 列及之后——因为之前列中已经存储了 L 的乘数，它们也需要跟着行一起转移。
2.3 全主元 vs 部分主元
| 策略 | 搜索范围 | 复杂度 | 使用场景 |
|---|---|---|---|
| 部分主元 | 第 k 列 [k:N-1] 行 | O(N²) | 绝大多数矩阵（工业标准） |
| 全主元 | 右下[k:N-1]×[k:N-1] | O(N³) | 病态矩阵、极限精度需求 |
部分主元在实践中几乎总是足够。
3. 前向替换——L 是单位下三角的特权
3.1 前向替换的核心公式
Ly = Pb，L 是单位下三角（对角线全为 1）：
y[0] = Pb[0] ← 只依赖 Pb[0]
y[1] = Pb[1] - L[1][0]·y[0] ← 依赖 y[0]
y[2] = Pb[2] - L[2][0]·y[0] - L[2][1]·y[1] ← 依赖 y[0], y[1]
一般形式: y[i] = Pb[i] - Σ(j=0 to i-1) L[i][j]·y[j]c
void forward_subst(double L[N][N], double b[N], int p[N], double y[N]) {
/* 用排列向量重排 b */
double Pb[N];
for (int i = 0; i < N; i++)
Pb[i] = b[p[i]];
/* 前向替换 */
for (int i = 0; i < N; i++) {
y[i] = Pb[i];
for (int j = 0; j < i; j++)
y[i] -= L[i][j] * y[j];
/* L[i][i] == 1，所以不需要除法 */
}
}3.2 本题的计算过程
p = [2, 0, 1] → Pb = [b[2], b[0], b[1]] = [29, 5, 11]
y₀ = 29
y₁ = 5 - 0.25 × 29 = 5 - 7.25 = -2.25
y₂ = 11 - 0.5 × 29 - 0.6667 × (-2.25)
= 11 - 14.5 + 1.5 = -2.0
y = [29, -2.25, -2.0]CAUTION
初学者最容易忘记的一步：先用 p 重排 b，再做前向替换。不能直接在原始 b 上计算——因为 A 的行已经被交换了（PA=LU），b 也必须按同样的方式重排。
4. 后向替换——不能忘记的除法
4.1 后向替换的核心公式
Ux = y，U 是上三角（对角元非零）：
x[2] = y[2] / U[2][2] ← 只依赖 y[2]
x[1] = (y[1] - U[1][2]·x[2]) / U[1][1] ← 依赖 x[2]
x[0] = (y[0] - U[0][1]·x[1] - U[0][2]·x[2]) / U[0][0] ← 依赖 x[1], x[2]
一般形式: x[i] = (y[i] - Σ(j=i+1 to N-1) U[i][j]·x[j]) / U[i][i]c
void back_subst(double U[N][N], double y[N], double x[N]) {
for (int i = N - 1; i >= 0; i--) {
x[i] = y[i];
for (int j = i + 1; j < N; j++)
x[i] -= U[i][j] * x[j];
x[i] /= U[i][i]; // ⚠ 除以对角元——这是最容易被遗忘的一步！
}
}本题计算结果：x₂=-2.0/(-0.6667)=3.0, x₁=(-2.25+3.75)/(-0.75)=-2.0, x₀=(29+14-27)/8=2.0 → x=[2,-2,3]
前向 vs 后向对比：前向顶部开始、依赖前面、不需除法（L[i][i]=1）；后向底部开始、依赖后面、必须除 U[i][i]。
NOTE
前向替换不需要除法的原因：L 是单位下三角，对角线固定为 1。后向替换必须除法的原因：U 的对角线元素一般不为 1，且由于选主元它们也不为 0（保证了除法安全）。
5. 浮点精度与数值稳定性
浮点运算有三类误差来源，部分主元是抑制它们的关键：
- 舍入误差：有限位数无法精确表示某些实数（如 1/3），用 double（64 位，约 15-17 位十进制有效数字）可基本覆盖
- 抵消误差：相近大数相减损失有效数字，部分主元将大元素放在对角线上减少抵消
- 放大误差：大乘数放大之前步骤的误差，部分主元保证所有 |l_ik| <= 1，误差最多保持原样不增长
c
double compute_residual(double A[N][N], double x[N], double b[N]) {
double residual = 0.0;
for (int i = 0; i < N; i++) {
double sum = 0.0;
for (int j = 0; j < N; j++) sum += A[i][j] * x[j];
double diff = sum - b[i];
residual += diff * diff;
}
return residual; // 本题残差应在 1e-30 量级
}float（32位，6-7位有效）vs double（64位，15-17位有效）：数值计算中几乎总是用 double。
6. 排列向量与排列矩阵
排列矩阵 P 是正交矩阵（P-1 = P^T），每行每列恰好一个 1。用长度为 N 的排列向量 p 即可隐式表示它：
p = [2, 0, 1] 对应 P:
| 0 0 1 | (第0行: 1在第2列)
P = | 1 0 0 | (第1行: 1在第0列)
| 0 1 0 | (第2行: 1在第1列)
作用: (Pv)[i] = v[p[i]], 即 Pb[i] = b[p[i]]排列向量的演变：
初始: p = [0, 1, 2] (恒等排列)
k=0: swap row0↔row2 → p = [2, 1, 0]
k=1: swap row1↔row2 → p = [2, 0, 1] (最终)p 的语义：p[i] = "当前第 i 行中存的是原矩阵的第几行"。
行列式计算：det(A) = (-1)^(#交换次数) * Prod U[i][i]。本题：2次交换×8×(-0.75)×(-0.6667)=4.0，与直接计算一致。
参考解答
TODO 1: print_matrix — 格式化打印矩阵
c
#include <math.h>
#include <stdio.h>
#define N 3
void print_matrix(const char *label, double m[N][N]) {
printf("%s:\n", label);
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
printf("%8.4f ", m[i][j]);
}
printf("\n");
}
}要点：格式 %8.4f 保证 8 位宽度含 4 位小数，每行一个换行。
TODO 2: print_vector — 格式化打印向量
c
void print_vector(const char *label, double v[N]) {
printf("%s: [", label);
for (int i = 0; i < N; i++) {
printf("%8.4f", v[i]);
if (i < N - 1) printf(" ");
}
printf("]\n");
}要点：最后一个元素后不加空格，括号外换行。
TODO 3: pivot — 部分主元选择与行交换
c
void pivot(double A[N][N], int p[N], int k) {
int max_row = k;
double max_val = fabs(A[k][k]); /* 初始化为 |A[k][k]| */
/* 在 [k, N-1] 行中寻找绝对值最大的 */
for (int i = k + 1; i < N; i++) {
double abs_val = fabs(A[i][k]);
if (abs_val > max_val) {
max_val = abs_val;
max_row = i;
}
}
if (max_row != k) {
/* 如果 A 的整行（所有 N 列） */
for (int j = 0; j < N; j++) {
double tmp = A[k][j];
A[k][j] = A[max_row][j];
A[max_row][j] = tmp;
}
/* 同步更新排列向量 */
int tmp_p = p[k];
p[k] = p[max_row];
p[max_row] = tmp_p;
printf(" Pivot: swap row %d <-> row %d\n", k, max_row);
}
}要点：交换整行而非只交换部分列；max_val 初始化为 fabs(A[k][k])，不能是 0（当所有元素为负时，0 不会被更新）。
TODO 4: lu_decomp — PA = LU 分解
c
void lu_decomp(double A[N][N], double L[N][N], double U[N][N], int p[N]) {
/* 步骤 1: 初始化排列为恒等 */
for (int i = 0; i < N; i++)
p[i] = i;
/* 步骤 2: 逐列分解 */
for (int k = 0; k < N; k++) {
pivot(A, p, k); /* 选主元 */
for (int i = k + 1; i < N; i++) {
A[i][k] /= A[k][k]; /* 乘数 l_ik */
for (int j = k + 1; j < N; j++)
A[i][j] -= A[i][k] * A[k][j]; /* Schur 补更新 */
}
}
/* 步骤 3: 从修改后的 A 中提取 L 和 U */
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
if (i > j) { L[i][j] = A[i][j]; U[i][j] = 0.0; }
else if (i == j) { L[i][j] = 1.0; U[i][j] = A[i][j]; }
else { L[i][j] = 0.0; U[i][j] = A[i][j]; }
}
}
}要点：先 pivot 再消去——这个顺序是关键。提取时 L 对角线固定为 1，U 对角线取自 A[i][i]。
TODO 5: forward_subst — 前向替换 Ly = Pb
c
void forward_subst(double L[N][N], double b[N], int p[N], double y[N]) {
double Pb[N];
for (int i = 0; i < N; i++)
Pb[i] = b[p[i]]; /* 用排列向量重排 b */
for (int i = 0; i < N; i++) {
y[i] = Pb[i];
for (int j = 0; j < i; j++)
y[i] -= L[i][j] * y[j];
/* L[i][i] == 1，无需除法 */
}
}要点：必须先 Pb[i] = b[p[i]] 将 b 按排列重排；L 的对角线为 1 所以不需要除法。
TODO 6: back_subst — 后向替换 Ux = y
c
void back_subst(double U[N][N], double y[N], double x[N]) {
for (int i = N - 1; i >= 0; i--) {
x[i] = y[i];
for (int j = i + 1; j < N; j++)
x[i] -= U[i][j] * x[j];
x[i] /= U[i][i]; /* 除以对角元——不能忘！ */
}
}要点：从底部向上求解；最后除以 U[i][i] 这一步最容易被遗忘。
TODO 7: solve + main — 主流程编排
c
#include <math.h>
#include <stdio.h>
#define N 3
/* 前文 print_matrix, print_vector, pivot, lu_decomp,
forward_subst, back_subst 的完整实现 */
void solve(double A_orig[N][N], double b[N], double x[N]) {
/* 复制 A_orig 到局部 A（分解会修改 A） */
double A[N][N];
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
A[i][j] = A_orig[i][j];
double L[N][N], U[N][N], y[N];
int p[N];
printf("=== LU Decomposition Solver with Partial Pivoting ===\n\n");
/* 步骤 1: LU 分解 */
printf("Step 1: PA = LU decomposition\n");
lu_decomp(A, L, U, p);
printf("\nPermutation vector p: [%d %d %d]\n\n",
p[0], p[1], p[2]);
print_matrix("L matrix", L);
printf("\n");
print_matrix("U matrix", U);
printf("\n");
/* 步骤 2: 前向替换 */
printf("Step 2: Forward substitution Ly = Pb\n");
forward_subst(L, b, p, y);
print_vector("y", y);
printf("\n");
/* 步骤 3: 后向替换 */
printf("Step 3: Backward substitution Ux = y\n");
back_subst(U, y, x);
print_vector("x (solution)", x);
printf("\n");
/* 验证 Ax = b */
printf("=== Verification: Ax ===\n");
double Ax[N];
for (int i = 0; i < N; i++) {
Ax[i] = 0.0;
for (int j = 0; j < N; j++)
Ax[i] += A_orig[i][j] * x[j];
}
print_vector("Ax", Ax);
print_vector("b", b);
}
int main(void) {
double A[N][N] = {{2, 1, 1},
{4, 3, 3},
{8, 7, 9}};
double b[N] = {5, 11, 29};
double x[N];
solve(A, b, x);
return 0;
}核心逻辑解析：
- 复制原矩阵：
solve内部拷贝一份A_orig到A，因为lu_decomp会原地修改它。验证时仍需原矩阵计算 Ax。 - 三步求解：PA=LU 分解（O(N³)）→ 前向替换（O(N²)）→ 后向替换（O(N²)）。当 N 很大或需要多次求解不同 b 时，分解只需做一次。
- 排列向量的使用：
p在pivot中被更新，在forward_subst中用于重排b。这是 PA=LU 中 P 的实现方式。 - 验证环节：用原始
A_orig和求出的x计算 Ax，与原始b对比——这是验证解正确性的标准方法。
对照检查：
pivot中交换了整行吗？lu_decomp中"先 pivot 再消去"的顺序对吗？提取 L/U 时对角线分配正确吗？前向替换前先重排了 b 吗？后向替换最后除以U[i][i]了吗？验证环节用了A_orig而非被修改的A吗？
课堂讨论
- 如果不选主元，用本题的矩阵直接做高斯消元——分解结果还正确吗？如果不正确，在哪一步出错？
- 部分主元搜索的是整列的最大绝对值。如果将该策略改为"搜索第一个非零元素"（first non-zero），会出现什么问题？
- 前向替换中 L 的对角线是 1，不需要除法。如果没有部分主元，U 的对角线可能为零——这意味着什么？后向替换能处理吗？
- 排列向量 p 记录了行交换历史。PA = LU 中，P 是排列矩阵。p = [2, 0, 1] 对应的 3×3 P 矩阵是什么？验证 P × b = Pb。
- 本题固定 N=3。如果改为 N=100 的通用求解器，哪些函数需要修改？哪些可以保持不变？
- 对于对称正定矩阵（如 [[4, 1, 1], [1, 3, 1], [1, 1, 2]]），可以用更高效的 Cholesky 分解（A = LLᵀ，无需选主元）。为什么本题的矩阵不能用 Cholesky 分解？
讨论问题
Q1: 不选主元会怎样？
本题 A = [[2,1,1],[4,3,3],[8,7,9]] 不选主元时：
k=0: 主元 A[0][0]=2, 消去:
l₁₀=4/2=2: row₁ -= 2·row₀ → [0, 1, 1]
l₂₀=8/2=4: row₂ -= 4·row₀ → [0, 3, 5]
k=1: 主元 A[1][1]=1, 消去:
l₂₁=3/1=3: row₂ -= 3·row₁ → [0, 0, 2]
L = [[1,0,0], U = [[2,1,1],
[2,1,0], [0,1,1],
[4,3,1]] [0,0,2]]
L×U = [[2,1,1], 验证: 4=4 ✓, 3=3 ✓
[4,3,3], 8=8 ✓, 7=7 ✓, 9=9 ✓
[8,7,9]] ≈ A ✓
本题恰好正确！因为主元 2 和 1 都不算小，
且各行元素本身已经是类似数量级。不选主元正确与否取决于矩阵如何。当各主元足够大时正确，但一旦遇到小主元就灾难性失败。这是"数值算法"与"精确数学"的核心区别——在浮点运算下，正确性可能取决于输入数据的条件不是算法逻辑。
Q2: 选第一个非零元素会怎样？
"第一个非零"策略在两个层面上失败：
问题一: 非零但极小
A = [[1e-16, 1], 第一个非零: A[0][0]=1e-16
[1, 1]]
乘数 l₂₁ = 1/1e-16 = 1e16 → 灾难性抵消，同不选主元
问题二: 非零但和其他元素差很多数量级
A = [[0.001, 1000],
[1, 1]]
不选主元: 0.001 作为主元, 乘数 l₂₁=1000 → 可能不安全
选主元: 1 > 0.001 → swap → 乘数 l₂₁=0.001 → 安全
结论: "非零"不等于"安全"。"最大绝对值"才是数值稳定的选择。Q3: U 对角元为零的含义
如果 U 的对角元为 0，意味着矩阵 A 是奇异的（行列式为 0），方程组要么无解要么有无穷多解。
U = [[2, 1, 1],
[0, 0, 3], ← U[1][1] = 0 — 矩阵奇异！
[0, 0, 1]]
后向替换遇到 x₁ 时:
x₁ = (y₁ - 3·x₂) / 0 → 除以零！
部分主元可以避免此问题（只要矩阵非奇异），
因为 pivot 总是选择绝对值最大的元素作为主元。Q4: 排列矩阵 P 的构造
p = [2, 0, 1] → P[i][j] = 1 当 j == p[i] 时，否则 0
┌ ┐
P = │ 0 0 1 │ (第 0 行: p[0]=2)
│ 1 0 0 │ (第 1 行: p[1]=0)
│ 0 1 0 │ (第 2 行: p[2]=1)
└ ┘
b = [5, 11, 29]ᵀ
P × b = [b[2], b[0], b[1]]ᵀ = [29, 5, 11]ᵀ
C 实现: Pb[i] = b[p[i]]Q5: 通用 N 的改造
c
/* 将 #define N 3 改为 #define N 100 */
/* 需要修改的函数: */
/* ① print_matrix: 无需修改——两个循环都到 N */
/* ② print_vector: 无需修改——循环到 N */
/* al pivot: 无需修改——循环条件都是 i < N, j < N */
/* ④ lu_decomp: 无需修改——三重循环都到 N */
/* ⑤ forward_subst: 无需修改——循环到 N */
/* ⑥ back_subst: 无需修改——循环到 N */
/* ⑦ solve: 只需修改矩阵/向量的初始化方式 */
/* 不再用字面量初始化 A 和 b，需要从输入读取 */
/* 同时需要注意:
* - 对于大 N，A[N][N] 可能超出栈空间（N>500 时）
* 应改用 malloc 在堆上分配
* - 编译时常量 N 可通过 -DN=100 传入
*/核心发现：算法的结构与 N 无关——这正是 LU 分解的通用性。
Q6: 为什么不用 Cholesky？
Cholesky 分解 A = LLᵀ 仅适用于对称正定矩阵。本题的 A：
A = [[2, 1, 1],
[4, 3, 3],
[8, 7, 9]]
不满足对称: A[1][0]=4 ≠ A[0][1]=1 ✗
不满足正定: 不需要进一步检查
Cholesky 的优精度:
- 计算量约为 LU 的一半
- 无需选主元（对称正定保证了数值稳定）
Cholesky 的限制:
- 仅适用于于对称正定矩阵
- 遇到非对称矩阵时必须用 LU
```c [cholesky_vs_lu.c]
/* Cholesky 分解示意（仅对称正定矩阵可用）:
* for i in 0..N-1:
* for j in 0..i:
* sum = A[i][j]
* for k in 0..j-1: sum -= L[i][k] * L[j][k]
* if i==j: L[i][j] = sqrt(sum)
* else: L[i][j] = sum / L[j][j]
*
* LU 分解（通用）:
* 适用于任何非奇异方阵，需选主元保证数值稳定
*/课后练习
不选主元的 LU 分解。编写
lu_decomp_no_pivot(A, L, U)，不做 pivot，总是用 A[k][k] 作为主元。用以下两组矩阵测试：（a）本课的 A；（b）A = [[1e-16, 1], [1, 1]]。对比两种结果的精度差异。知识点提示：对（a）结果可能恰好正确，但对（b）应观察到显著的精度退化——验证 |Ax - b|。浮点运算下，正确性取决于输入数据的条件，而非算法逻辑本身。
参考解答
c#include <math.h> #include <stdio.h> #define N 3 void lu_decomp_no_pivot(double A[N][N], double L[N][N], double U[N][N]) { for (int k = 0; k < N; k++) { for (int i = k + 1; i < N; i++) { A[i][k] /= A[k][k]; for (int j = k + 1; j < N; j++) A[i][j] -= A[i][k] * A[k][j]; } } /* 提取 L/U（同上） */ for (int i = 0; i < N; i++) for (int j = 0; j < N; j++) { if (i > j) { L[i][j] = A[i][j]; U[i][j] = 0; } else if (i == j){ L[i][j] = 1; U[i][j] = A[i][j]; } else { L[i][j] = 0; U[i][j] = A[i][j]; } } } int main(void) { double A[N][N] = {{2,1,1},{4,3,3},{8,7,9}}; double L[N][N], U[N][N]; lu_decomp_no_pivot(A, L, U); printf("No-pivot LU for course matrix:\n"); /* 打印 L, U, 验证 L*U ≈ A */ return 0; }通用规模的 LU 求解器。修改所有函数，使 N 可作为运行时参数（通过
#define N 100或函数参数传递）。从标准输入读取 N、A、b，输出解 x。知识点提示：对于大 N，栈空间的
double A[N][N]可能会溢出。改用double **A = malloc(N * sizeof(double*))并逐行分配。注意这改变了函数签名中矩阵参数的类型（从double m[N][N]到double **m）。参考解答
c#include <math.h> #include <stdio.h> #include <stdlib.h> /* 分配 N×N 矩阵 */ double **alloc_matrix(int n) { double **m = malloc(n * sizeof(double *)); for (int i = 0; i < n; i++) m[i] = malloc(n * sizeof(double)); return m; } void free_matrix(double **m, int n) { for (int i = 0; i < n; i++) free(m[i]); free(m); } void pivot(double **A, int *p, int n, int k) { int max_row = k; double max_val = fabs(A[k][k]); for (int i = k + 1; i < n; i++) { if (fabs(A[i][k]) > max_val) { max_val = fabs(A[i][k]); max_row = i; } } if (max_row != k) { for (int j = 0; j < n; j++) { double tmp = A[k][j]; A[k][j] = A[max_row][j]; A[max_row][j] = tmp; } int tmp = p[k]; p[k] = p[max_row]; p[max_row] = tmp; } } /* lu_decomp, forward_subst, back_subst 类似地增加 int n 参数 */残差计算。编写
double residual(double A[N][N], double x[N], double b[N])，计算 ‖Ax - b‖₂（欧几里得范数）。对本课的解和故意扰动 0.01 的"错误解"，分别输出残差。观察正确解的残差数量级（应接近机器精度 ~1e-16）。知识点提示：残差 ‖Ax - b‖ 是衡量解质量的直接度量——逼近 0 越好。用
sqrt(Σ (Ax[i]-b[i])²)计算。参考解答
c#include <math.h> #include <stdio.h> #define N 3 double residual(double A[N][N], double x[N], double b[N]) { double sum = 0.0; for (int i = 0; i < N; i++) { double ax = 0.0; for (int j = 0; j < N; j++) ax += A[i][j] * x[j]; double diff = ax - b[i]; sum += diff * diff; } return sqrt(sum); } int main(void) { double A[N][N] = {{2,1,1},{4,3,3},{8,7,9}}; double b[N] = {5, 11, 29}; double x_correct[N] = {2.0, -2.0, 3.0}; double x_wrong[N] = {2.01, -2.0, 3.0}; printf("Correct solution residual: %e\n", residual(A, x_correct, b)); printf("Perturbed solution residual: %e\n", residual(A, x_wrong, b)); return 0; } /* 预期: * Correct residual: ~1e-15 或更小 * Perturbed residual: ~0.1 左右 */排列矩阵显式构造。编写函数
void build_permutation(int p[N], double P[N][N])，从排列向量 p 显式构造 N×N 排列矩阵 P。用 P 验证 PA ≈ LU （即 P·A ≈ L·U）。知识点提示：
P[i][j] = (j == p[i]) ? 1.0 : 0.0。矩阵乘法用三重循环计算 P × A 和 L × U，逐元素比较。参考解答
c#define N 3 void build_permutation(int p[N], double P[N][N]) { for (int i = 0; i < N; i++) for (int j = 0; j < N; j++) P[i][j] = (j == p[i]) ? 1.0 : 0.0; } void mat_mul(double A[N][N], double B[N][N], double C[N][N]) { for (int i = 0; i < N; i++) for (int j = 0; j < N; j++) { C[i][j] = 0.0; for (int k = 0; k < N; k++) C[i][j] += A[i][k] * B[k][j]; } } /* 使用: double P[N][N], LU[N][N], PA[N][N]; build_permutation(p, P); mat_mul(P, A_orig, PA); mat_mul(L, U, LU); // 比较 PA 和 LU — 应几乎相等 */行列式计算器。扩展 lu_decomp，增加计算行列式的功能。利用 p 中的交换次数和 U 的对角线乘积：det(A) = (-1)^(交换次数) × ∏ U[i][i]。
知识点提示：在 pivot 中维护交换计数器。分解完成后，符号 × U 对角元乘积即为行列式。验证：用本课 A 应得 4.0。
参考解答
cdouble determinant_from_lu(int swaps, double U[N][N]) { int sign = (swaps % 2 == 0) ? 1 : -1; double det = sign; for (int i = 0; i < N; i++) det *= U[i][i]; return det; } /* 在 lu_decomp 中增加 int swaps = 0; 每次 pivot 交换时将 swaps++ */
参考资料
- Trefethen & Bau, Numerical Linear Algebra, SIAM 1997. Chapters 20-21 — LU 分解与选主元的权威讲解
- Golub & Van Loan, Matrix Computations, 4th ed. Chapter 3 — 线性系统求解的全面覆盖
- Wilkinson, J.H., The Algebraic Eigenvalue Problem, 1965 — 部分主元分析的奠基之作
- LAPACK 文档:
dgetrf(LU 分解),dgetrs(用 LU 求解) — 工业标准实现 - MATLAB:
[L,U,P] = lu(A)— 交互式验证 LU 分解的便捷工具
"The purpose of computing is insight, not numbers." — Richard Hamming