1. 为什么我要手写迭代法而不是直接调库

数值分析这门课里，线性方程组的求解是绕不开的一道坎。当年我第一次接触高斯消元的时候觉得挺美，一个消元回代就出结果，多干净。直到后来做一个热传导仿真的小项目，网格稍微划细一点，系数矩阵直接膨胀到几千阶，高斯消元跑一次要等好几分钟，而且内存占用肉眼可见地往上飙。那时候我才真正理解为什么工程上大量使用迭代解法——它不需要存储整个矩阵的LU分解结果，迭代过程中矩阵可以保持稀疏，内存占用小，而且很多时候收敛速度完全够用。

这篇文章要聊的就是用C语言从零实现三种最经典的线性方程组迭代解法：雅可比迭代法、高斯-赛德尔迭代法和超松弛迭代法。我会把每种方法的数学原理、C语言实现细节、收敛条件判断、以及实际跑的时候会遇到的各种坑都掰开揉碎讲清楚。适合正在学数值分析的学生、需要自己写求解器的工程师、以及任何对迭代法感兴趣但一直没动手写过的朋友。代码我会给完整的可编译版本，参数怎么选、什么时候该换方法、精度怎么控制，这些实战经验我都会一并交代。

先说清楚定位：这不是一篇数学推导文，我不会花大篇幅证明收敛性定理。我要做的是让你看完之后能直接打开编辑器，把代码敲进去跑起来，并且知道每一行在干什么、为什么这么写、出了问题去哪里找原因。

2. 三种迭代解法的核心思路与选型逻辑

2.1 迭代法的基本思想：从猜测到逼近

线性方程组写成矩阵形式就是Ax = b，其中A是系数矩阵，x是未知数向量，b是右端项。直接法（比如高斯消元）是一步到位算出精确解（在浮点精度范围内），而迭代法的思路完全不同：先随便猜一个初始解x⁽⁰⁾，然后构造一个迭代格式，让每一步得到的x⁽ᵏ⁺¹⁾都比上一步更接近真实解，直到相邻两次迭代的差值小到可以接受为止。

这个思路用生活化的类比来说就是：你在一个完全黑暗的房间里找一盏灯，你不知道灯具体在哪，但你能感觉到光的方向。每次你朝光更亮的方向挪一步，挪着挪着就到了灯旁边。迭代法就是这样，每次修正一点点，最终逼近真实解。

把A拆成三部分：A = D - L - U，其中D是对角矩阵（只保留A的对角线元素），L是严格下三角部分取负号，U是严格上三角部分取负号。这个拆分是理解所有迭代法的关键。基于这个拆分，可以构造出不同的迭代格式，对应不同的方法。

2.2 雅可比迭代：最直观但收敛最慢

雅可比迭代的公式非常直观。对于第i个方程，把除了x_i以外的所有项移到右边，得到：

x_i⁽ᵏ⁺¹⁾ = (b_i - Σ_{j≠i} a_ij · x_j⁽ᵏ⁾) / a_ii

注意右边用的全是上一步的旧值x⁽ᵏ⁾，这就是雅可比迭代的特征——同步更新。所有分量都算完之后才统一替换。写成矩阵形式就是x⁽ᵏ⁺¹⁾ = D⁻¹(L+U)x⁽ᵏ⁾ + D⁻¹b。

雅可比迭代的优点是实现简单、天然适合并行计算（每个分量互不依赖），缺点是收敛速度通常最慢，而且对矩阵性质要求比较苛刻。它收敛的充分条件是矩阵A严格对角占优，即每一行对角元素的绝对值大于该行其他元素绝对值之和。

2.3 高斯-赛德尔迭代：用上新值加速收敛

高斯-赛德尔迭代在雅可比的基础上做了一个看似微小但影响巨大的改动：计算x_i⁽ᵏ⁺¹⁾的时候，如果x_1到x_{i-1}的新值已经算出来了，就直接用新值而不是旧值。公式变成：

x_i⁽ᵏ⁺¹⁾ = (b_i - Σ_{j<i} a_ij · x_j⁽ᵏ⁺¹⁾ - Σ_{j>i} a_ij · x_j⁽ᵏ⁾) / a_ii

这个改动的好处是收敛速度通常比雅可比快一倍左右，而且不需要额外的存储空间（直接在原数组上更新即可）。代价是失去了并行性，因为每个分量依赖前一个分量的新值。矩阵形式是x⁽ᵏ⁺¹⁾ = (D-L)⁻¹Ux⁽ᵏ⁾ + (D-L)⁻¹b。

2.4 超松弛迭代：加个松弛因子再快一点

超松弛迭代（SOR）在高斯-赛德尔的基础上引入了一个松弛因子ω：

x_i⁽ᵏ⁺¹⁾ = (1-ω)x_i⁽ᵏ⁾ + ω·x_i^{GS}

其中x_i^{GS}是高斯-赛德尔算出来的值。当ω=1时退化为高斯-赛德尔；当0<ω<1时叫低松弛，用于处理某些不收敛的情况；当1<ω<2时叫超松弛，用于加速收敛。理论上存在一个最优的ω值能让收敛速度最大化，对于特定的矩阵类型（比如三对角矩阵）甚至有解析公式可以算出来。

选型逻辑很清晰：如果矩阵严格对角占优且规模不大，雅可比就够了；如果想要更快的收敛且不介意串行执行，用高斯-赛德尔；如果对收敛速度有更高要求且愿意花时间调参，上SOR。实际工程中高斯-赛德尔和SOR用得最多。

3. C语言实现的核心细节与避坑指南

3.1 数据结构设计：二维数组还是稀疏存储

先解决一个最基本的问题：矩阵怎么存。对于教学和小规模问题，直接用二维数组double A[N][N]最省事。但如果你要处理几千阶的稀疏矩阵，二维数组就太浪费了——大部分元素都是零，存它干嘛。

我一般会分两种情况处理：

• 稠密矩阵：直接用二维数组，代码简单，缓存友好。N在几百以内都没问题。
• 稀疏矩阵：用CSR（压缩稀疏行）格式存储，只存非零元素的值、列索引和行偏移。迭代法的一大优势就是可以充分利用稀疏性，每次迭代只遍历非零元素。

为了不让代码过于复杂，这篇文章的示例代码用二维数组，但我会在关键位置说明如果换成CSR格式该怎么改。

#include <stdio.h> #include <stdlib.h> #include <math.h> #include <string.h> #define N 100 // 方程组阶数 #define MAX_ITER 10000 // 最大迭代次数 #define TOL 1e-10 // 收敛容差 // 全局矩阵和向量 double A[N][N]; double b[N], x[N], x_new[N];

这里有几个细节值得说。TOL取1e-10是比较严格的标准，实际工程中1e-6到1e-8就够用了。MAX_ITER设10000是防止不收敛时死循环，这个必须设，我见过太多人忘了设最大迭代次数，结果程序跑了一晚上还在转。

3.2 雅可比迭代的C语言实现

雅可比迭代的核心就是两层循环：外层控制迭代次数，内层遍历每个未知数。注意内层必须用旧值计算，所以需要一个额外的数组x_new来暂存新值。

int jacobi(int n, double A[N][N], double *b, double *x, int max_iter, double tol) { double *x_new = (double *)malloc(n * sizeof(double)); if (!x_new) return -1; int iter; for (iter = 0; iter < max_iter; iter++) { double max_diff = 0.0; // 对每个未知数计算新值 for (int i = 0; i < n; i++) { double sum = 0.0; for (int j = 0; j < n; j++) { if (j != i) { sum += A[i][j] * x[j]; } } x_new[i] = (b[i] - sum) / A[i][i]; } // 计算最大变化量并更新解向量 for (int i = 0; i < n; i++) { double diff = fabs(x_new[i] - x[i]); if (diff > max_diff) max_diff = diff; x[i] = x_new[i]; } // 判断收敛 if (max_diff < tol) { free(x_new); return iter + 1; // 返回迭代次数 } } free(x_new); return -1; // 未收敛 }

这段代码里有个容易忽略的点：收敛判断用的是最大变化量而不是欧几里得范数。用最大范数（无穷范数）的好处是计算量小，不需要开方，而且对于大多数问题来说判断效果和2-范数差不多。如果你对精度要求特别高，可以改用2-范数，但每次迭代要多算一次开方。

另一个细节是A[i][i]不能为零。如果对角线上有零元素，程序直接除零崩溃。实际使用前应该检查一下，如果对角元为零，要么做行交换，要么换方法。

3.3 高斯-赛德尔迭代的实现差异

高斯-赛德尔和雅可比的代码结构几乎一样，唯一的区别是不需要x_new数组，直接在x上更新：

int gauss_seidel(int n, double A[N][N], double *b, double *x, int max_iter, double tol) { int iter; for (iter = 0; iter < max_iter; iter++) { double max_diff = 0.0; for (int i = 0; i < n; i++) { double sum = 0.0; for (int j = 0; j < n; j++) { if (j != i) { sum += A[i][j] * x[j]; // 注意：j<i时x[j]已是新值 } } double x_new = (b[i] - sum) / A[i][i]; double diff = fabs(x_new - x[i]); if (diff > max_diff) max_diff = diff; x[i] = x_new; // 立即更新 } if (max_diff < tol) return iter + 1; } return -1; }

看起来只改了一点点，但收敛行为完全不同。高斯-赛德尔之所以更快，是因为它利用了最新信息。你可以这样理解：雅可比是每步所有人都用昨天的信息做决策，高斯-赛德尔是每个人做决策时已经知道了前面人的最新决策，信息更及时，决策更准。

3.4 SOR迭代：松弛因子的引入

SOR的代码在高斯-赛德尔基础上加一行加权平均：

int sor(int n, double A[N][N], double *b, double *x, double omega, int max_iter, double tol) { int iter; for (iter = 0; iter < max_iter; iter++) { double max_diff = 0.0; for (int i = 0; i < n; i++) { double sum = 0.0; for (int j = 0; j < n; j++) { if (j != i) { sum += A[i][j] * x[j]; } } double x_gs = (b[i] - sum) / A[i][i]; double x_new = (1.0 - omega) * x[i] + omega * x_gs; double diff = fabs(x_new - x[i]); if (diff > max_diff) max_diff = diff; x[i] = x_new; } if (max_diff < tol) return iter + 1; } return -1; }

ω的取值是SOR的灵魂。理论最优值ω_opt = 2 / (1 + sqrt(1 - ρ²))，其中ρ是雅可比迭代矩阵的谱半径。但实际中ρ很难算，所以工程上常用试凑法：从1.0开始，以0.1为步长往上试，找到迭代次数最少的那个值。我试过的大部分问题，ω在1.2到1.8之间效果最好。

注意：ω必须严格小于2，否则迭代一定发散。ω=2时公式退化为x_new = 2*x_gs - x_old，这是外推，不稳定。

4. 完整可运行的测试程序与参数调优

4.1 构造测试用例：对角占优矩阵

要测试迭代法，得先有一个保证收敛的矩阵。最简单的方法是构造严格对角占优矩阵：随机生成非对角元素，然后让对角元素等于该行其他元素绝对值之和加一个正数。

void generate_test_matrix(int n, double A[N][N], double *b, double *x_exact) { srand(42); // 固定种子，保证可复现 // 生成精确解（全1向量方便验证） for (int i = 0; i < n; i++) { x_exact[i] = 1.0; } // 生成对角占优矩阵 for (int i = 0; i < n; i++) { double row_sum = 0.0; for (int j = 0; j < n; j++) { if (j != i) { A[i][j] = (double)(rand() % 100) / 100.0 - 0.5; row_sum += fabs(A[i][j]); } } A[i][i] = row_sum + 1.0 + (double)(rand() % 100) / 100.0; } // 计算b = A * x_exact for (int i = 0; i < n; i++) { b[i] = 0.0; for (int j = 0; j < n; j++) { b[i] += A[i][j] * x_exact[j]; } } }

这个构造方法的好处是精确解已知（全是1），跑完迭代后可以直接对比误差。srand(42)固定随机种子是为了让每次运行结果一致，方便调试。

4.2 主函数与性能对比

int main() { double x_exact[N]; double x_jac[N], x_gs[N], x_sor[N]; generate_test_matrix(N, A, b, x_exact); // 初始猜测全零 memset(x_jac, 0, sizeof(double) * N); memset(x_gs, 0, sizeof(double) * N); memset(x_sor, 0, sizeof(double) * N); // 雅可比 int iter_jac = jacobi(N, A, b, x_jac, MAX_ITER, TOL); printf("Jacobi: iter=%d, error=%.2e\n", iter_jac, compute_max_error(x_jac, x_exact, N)); // 高斯-赛德尔 int iter_gs = gauss_seidel(N, A, b, x_gs, MAX_ITER, TOL); printf("Gauss-Seidel: iter=%d, error=%.2e\n", iter_gs, compute_max_error(x_gs, x_exact, N)); // SOR，试几个omega double omegas[] = {1.0, 1.2, 1.5, 1.8}; for (int k = 0; k < 4; k++) { memset(x_sor, 0, sizeof(double) * N); int iter_sor = sor(N, A, b, x_sor, omegas[k], MAX_ITER, TOL); printf("SOR(omega=%.1f): iter=%d, error=%.2e\n", omegas[k], iter_sor, compute_max_error(x_sor, x_exact, N)); } return 0; } double compute_max_error(double *x, double *x_exact, int n) { double max_err = 0.0; for (int i = 0; i < n; i++) { double err = fabs(x[i] - x_exact[i]); if (err > max_err) max_err = err; } return max_err; }

我实测跑N=100的结果大概是这样的：雅可比需要800多次迭代，高斯-赛德尔400次左右，SOR在ω=1.5时只要150次左右。这个差距随着N增大还会拉大。所以如果你的问题规模比较大，花点时间调ω是值得的。

4.3 收敛性判断与异常处理

实际使用中最大的坑就是不收敛。迭代法不是万能的，矩阵不满足收敛条件时它会一直震荡或者发散。除了设最大迭代次数，还应该监控残差的变化趋势：

// 在迭代循环内加入残差监控 double residual = 0.0; for (int i = 0; i < n; i++) { double r = b[i]; for (int j = 0; j < n; j++) { r -= A[i][j] * x[j]; } residual += r * r; } residual = sqrt(residual); // 如果残差开始增大，说明可能发散 if (iter > 10 && residual > prev_residual * 1.5) { printf("Warning: residual increasing, may diverge at iter %d\n", iter); break; } prev_residual = residual;

这个监控很实用。我遇到过好几次矩阵只是“弱对角占优”的情况，前几十步残差在降，后面突然开始涨，如果没有这个监控，程序会一直跑到最大迭代次数才停，浪费时间。

5. 常见问题排查与实战经验

5.1 迭代不收敛怎么办

这是被问得最多的问题。排查顺序我一般是这样：

对角占优是最关键的。如果矩阵不满足，雅可比和高斯-赛德尔都可能不收敛。这时候可以考虑对原方程做行交换，让对角元尽可能大。或者用SOR加低松弛（ω<1）来“稳住”迭代。

5.2 精度上不去的原因

有时候迭代明明收敛了，但误差就是降不到1e-10以下。常见原因有三个：

第一是浮点累积误差。每次迭代都有舍入误差，迭代次数多了误差会累积。这时候把TOL放宽到1e-8通常就够了，没必要追求机器精度。

第二是矩阵条件数太大。病态矩阵对扰动极其敏感，解本身就不稳定。这种情况迭代法救不了，得用预处理或者换直接法。

第三是收敛判据选得不好。用最大变化量判断时，如果解向量某些分量变化很小但误差很大，可能提前误判收敛。可以改用残差范数判断：||b - Ax|| < tol。这个判据更可靠，但每次迭代要多算一次矩阵向量乘法。

5.3 性能优化技巧

如果N很大，迭代法的性能瓶颈在矩阵向量乘法上。几个实用的优化手段：

• 利用稀疏性：如果矩阵大部分是零，用CSR格式存储，每次迭代只遍历非零元素。我做过一个对比，N=5000的稀疏矩阵，CSR格式比二维数组快20倍以上。
• 循环展开：内层循环手动展开4次，减少循环开销。现代编译器一般会自动做，但手动展开有时效果更好。
• OpenMP并行：雅可比迭代天然并行，加一行#pragma omp parallel for就能利用多核。高斯-赛德尔不能直接并行，但可以用红黑排序（red-black ordering）来并行化。
• 编译优化：用-O2或-O3编译，开启向量化。-ffast-math可以进一步加速但会牺牲一些精度，慎用。

提示：优化之前先用-pg编译跑一遍gprof，找到真正的热点再动手。我见过有人花半天优化一个只占5%时间的函数，纯属浪费时间。

5.4 初始猜测的选择

初始猜测对迭代次数影响很大。全零初值是最省事的，但如果真解离零很远，前几十次迭代基本在“赶路”。如果对问题有先验知识，比如知道解大概在什么范围，用那个范围的中心作为初值能省不少迭代。

一个实用的技巧是：先用低精度（TOL=1e-3）跑一遍，把结果作为高精度（TOL=1e-10）的初值。这样总时间通常比直接跑高精度要短，因为低精度阶段很快就收敛了，高精度阶段从很近的地方出发，迭代次数少。

6. 从代码到工程：我的几点体会

写完这三种迭代法并且在实际项目中用过之后，有几个体会特别深。

第一，迭代法的价值在于“够用就好”。很多工程问题不需要1e-15的精度，1e-6完全够用。在这个精度下，迭代法往往比直接法快得多，尤其是矩阵规模大的时候。我做过一个热传导问题，N=2000，高斯消元要跑好几分钟，高斯-赛德尔迭代30秒就收敛到1e-6了。

第二，收敛性比速度重要。一个收敛慢但稳定的方法，比一个快但偶尔发散的方法靠谱得多。这也是为什么高斯-赛德尔比雅可比更常用——它收敛更快，而且收敛条件更宽松（雅可比要求严格对角占优，高斯-赛德尔只要求对角占优或不可约弱对角占优）。

第三，调参是门手艺。SOR的ω值、收敛容差、最大迭代次数，这些参数没有万能值，得根据具体问题试。我的习惯是先跑一个小规模（N=50）的测试，把参数调好，再放大到实际规模。小规模跑一次几毫秒，大规模跑一次可能几分钟，先在小规模上试错成本低得多。

最后分享一个调试技巧：把每次迭代的解向量和残差都写到文件里，用Python或者gnuplot画出来。收敛过程一目了然，是单调下降还是震荡、收敛速度是线性还是超线性，看图比看数字直观得多。这个习惯帮我省了很多瞎猜的时间。