简介这是一套面向数值分析与计算方法学习者的zip压缩包以牛顿迭代法和雅可比迭代算法为核心重点解决非线性方程求根与大型稀疏线性方程组迭代求解问题同时收入直接迭代法、拉格朗日插值、牛顿插值等经典数值方法可作为计算物理、工程计算等课程的课后实践与复习资料。包内共44个文件压缩包体积仅2.36MB主要文件类型包括C/C源代码.cpp、.c、可直接运行的exe程序以及Visual Studio工程文件.dsp、.dsw、.opt、.plg、调试符号文件.pdb、.idb、.ncb和少量文本文件.txt既能运行验证也便于对照源码修改调试。目前已有548人学习下载代码注释自然清晰适合初学者将算法公式与程序实现一一对应也可供高年级学生快速复现牛顿、雅可比、拉格朗日插值、直接迭代等数值实验借助包内不同算例还可横向比较各迭代格式的收敛速度与适用条件从而系统梳理计算方法导引的核心思路。1. 计算方法导引里的两个经典迭代牛顿迭代法与雅可比迭代算法做《计算方法导引》这门课的人多半是从牛顿迭代法和雅可比迭代算法这两个名字开始真正接触迭代的。前者用来解非线性方程后者用来解线性方程组一个靠切线逼近一个把矩阵拆开轮流更新表面看没什么关系骨子里却是同一件事把求精确解的问题替换成“反复做同一个四则运算直到误差小到可以接受”。这个替换看起来简单坑却全藏在参数里——初值选在哪、迭代上限设多少、停机容差给多大直接决定结果是收敛、振荡还是干脆溢出。这篇文章把两个迭代从公式推到最小实现再沿收敛判据、失效修复、松弛加速一路落到可运行的 Python 代码上适合正在学数值计算或者需要在工程里手写迭代器、不能调库的人。2. 牛顿迭代法从泰勒展开到最小实现收敛判据怎么设2.1 切线法的几何直觉与迭代公式对单变量方程 f(x)0牛顿迭代法的常规推导是在当前近似解 x_k 附近做一阶泰勒展开f(x_{k1}) ≈ f(x_k) f(x_k)(x_{k1} - x_k)令右侧等于 0解出 x_{k1}得到迭代公式x_{k1} x_k - f(x_k) / f(x_k)几何上这一步等于过点 (x_k, f(x_k)) 做切线把切线与 x 轴的交点当作下一个近似解所以也叫切线法。这个公式成立的前提是 f(x_k) 不为 0如果在某步导数恰好接近 0下一步会被甩到很远后面会专门讨论这种情况。2.2 最小实现以 x^3 - 2x - 5 0 为例经典教材里常拿这个三次方程做演示它只有一个实根导数也好算。最小实现如下def newton(f, df, x0, tol1e-8, max_iter50): x x0 for k in range(max_iter): fx f(x) dfx df(x) if abs(dfx) 1e-12: raise ValueError(f第 {k} 步导数为 {dfx}牛顿法失效) x_new x - fx / dfx if abs(x_new - x) tol: return x_new, k 1 x x_new raise RuntimeError(f超过最大迭代次数 {max_iter}未收敛) if __name__ __main__: f lambda x: x**3 - 2*x - 5 df lambda x: 3*x**2 - 2 root, steps newton(f, df, x02.0) print(f根: {root:.10f}, 迭代次数: {steps})这段代码里x0 是初值tol 是停机阈值max_iter 是保护性上限。第 8 行用相邻两次迭代的差来判断收敛而不是用 f(x) 的绝对值原因是牛顿法在根附近残差与步长同时趋零但残差判断在 f(x) 本身数值很小但离根还很远的病态函数下会提前误判步长判断更接近“解不再移动”的直观含义。注意第 10 行的 raise 语句工程实现里必须给迭代加“安全带”否则死循环比发散更麻烦。以 x02.0 启动牛顿法大约 5 步就能把根压到 1e-8 以内比二分法快一个数量级。代价是每一步都要计算 f 和 f 两个函数值后续会看到这个代价在方程组场景会被放大。2.3 停机参数怎么配tol 与 max_iter 的经验值下面的表是单变量牛顿迭代较通用的参数区间可作为起点再按实际问题调整参数典型值说明x0根附近或变号区间内离根太远时极易振荡见第 3 章tol1e-8 ~ 1e-10低于 1e-14 会撞上双精度浮点表示极限max_iter20 ~ 50再往上通常不是收敛慢而是根本不收敛导数阈值1e-12低于该值直接报错防止除零和越界跳变tol 不建议小于 1e-14。64 位浮点的机器精度约 2.2e-16当相邻步差缩到接近 1e-15 时x_new 与 x 的差已经落在舍入误差范围里继续迭代只会原地抖动收敛判据设得比浮点精度还严结果就是一路跑满 max_iter 然后抛异常。如果确实需要高精度应该换高精度算术库而不是反复收紧 tol。3. 牛顿迭代法的失效场景初值发散、多重根降阶与阻尼修正3.1 初值选择如何决定收敛与否同一个方程 x^3 - 2x - 5 0把 x0 从 2.0 换成 0.0迭代序列会先跳到 -2.5随后在远离根的地方振荡发散。这不是算法写错了而是牛顿法只是局部收敛只有初值落在根的收敛半径内二次收敛才会发生。工程上最稳妥的做法是先做全局扫描找变号区间再用牛顿法精化。变号区间只需 30 行以内的朴素代码在区间 [-10, 10] 上按步长 0.1 采样检查 f(x_i) 与 f(x_{i1}) 异号的位置作为牛顿法的启动区间。这个扫描成本很低但能把“发散”这类问题拦截在迭代之前。对高维问题没有便宜的全局面扫描只能靠多初值并行启动加残差筛选。失效类型现象兜底手段初值远离根振荡或直接溢出全局扫描取变号区间多重根收敛变慢误差每步减半乘重数修正或阻尼导数接近零跨出有效范围阻尼牛顿或切回二分3.2 多重根附近收敛降阶带重数修正的牛顿迭代法牛顿法在单根附近有二次收敛速度但根的重数 m 大于 1 时收敛会退化为线性且每次迭代误差约乘 (m-1)/m。以 f(x) (x-1)^2 为例根 x1 的重数为 2迭代格式退化为 x_{k1} (x_k 1)/2误差每步减半肉眼可见地慢。修正方案是给导数项乘上重数 mx_{k1} x_k - m * f(x_k) / f(x_k)def newton_multi(f, df, m, x0, tol1e-10, max_iter50): x x0 for k in range(max_iter): x_new x - m * f(x) / df(x) if abs(x_new - x) tol: return x_new, k 1 x x_new raise RuntimeError(f超过最大迭代次数 {max_iter}) f lambda x: (x - 1.0) ** 2 df lambda x: 2.0 * (x - 1.0) print(newton_multi(f, df, m2, x00.5))这里的 m 通常要预先知道比如多项式因式分解或对函数结构做分析。如果 m 未知可以在根附近用 f(x)/f(x) 的比值估计不过工程上更常见的做法是用下面的阻尼措施而不是强行猜重数。另外提醒一句牛顿法碰到重根时残差下降会变慢但 x 的位移其实并没有完全停止此时用步长判据会比残差判据更可靠。3.3 阻尼牛顿与二分兜底把发散拉回来实际计算中初值很难保证一定落在收敛半径内。常见的兜底策略有两种。第一种是阻尼牛顿在每次更新时引入步长因子 alphax_{k1} x_k - alpha * f(x_k)/f(x_k)先取 alpha1如果 f 的绝对值没有下降就把 alpha 减半重试alpha 减到某个下界仍不下降就认为在当前点无法推进。第二种是退回二分法把扫描阶段得到的变号区间保留下来当牛顿迭代跳出区间时切回二分压几轮再交给牛顿。阻尼牛顿与二分法可以组合成混合求解器这也是 SciPy 里 brentq 类算法的工程动机。单变量问题没必要硬撑牛顿法左右横跳的代价很小真正需要精心调参数的是后面要讲的线性方程组迭代那里没有二分法可以兜底。4. 雅可比迭代算法并行友好的线性方程组求解器4.1 分量迭代格式与矩阵形式是怎么来的对于线性方程组 Axb雅可比迭代算法的出发点是把第 i 个方程单独拎出来用其他分量的旧值解出 x_ix_i^{(k1)} (b_i - Σ_{j≠i} a_ij x_j^{(k)}) / a_ii注意所有 x_i 的更新都只用第 k 步的值第 k1 步的新值在同一个 pass 内互不依赖因此雅可比迭代也叫同步迭代。整轮更新写成矩阵形式x^{(k1)} M x^{(k)} c其中 M I - D^{-1}Ac D^{-1}b这里 D 是取 A 对角线组成的对角矩阵。M 的非对角元为 -a_ij/a_ii对角元为 0。有的教材把 A 拆成 D-L-U再把迭代矩阵写成 D^{-1}(LU)符号约定不同写代码时容易对不上号直接用 M I - D^{-1}A 可以绕开这个坑。4.2 最小实现与对角占优检查为了后面复用这里用 NumPy 实现并配套一个严格对角占优检查函数import numpy as np def jacobi(A, b, x0None, tol1e-8, max_iter1000): A np.asarray(A, dtypefloat) b np.asarray(b, dtypefloat) n b.size x np.zeros(n) if x0 is None else np.array(x0, dtypefloat) x_next np.zeros_like(x) for k in range(max_iter): for i in range(n): s A[i] x - A[i, i] * x[i] x_next[i] (b[i] - s) / A[i, i] diff np.max(np.abs(x_next - x)) x[:] x_next if diff tol: return x, k 1 raise RuntimeError(f超过最大迭代次数 {max_iter}) def is_diagonally_dominant(A): A np.asarray(A, dtypefloat) for i in range(A.shape[0]): if abs(A[i, i]) np.sum(np.abs(A[i, :])) - abs(A[i, i]): return False return True A np.array([[4, 1, 1], [1, 4, 1], [1, 1, 4]], dtypefloat) b np.array([6, 6, 6], dtypefloat) print(is_diagonally_dominant(A)) print(jacobi(A, b, tol1e-10))代码里s A[i] x - A[i, i] * x[i]这行是核心先算整行点积再减掉自己的那一项正好对应分量公式里的求和。用两份数组 x 与 x_next 是必须的如果像高斯-赛德尔那样原地覆盖算法的性质就变了。is_diagonally_dominant 里的写法没有先取绝对值求和再减对角而是直接用总和减对角绝对值避免额外分配数组。4.3 收敛判定充分条件与充要条件分别看什么严格对角占优定义为对每一行|a_ii| Σ_{j≠i} |a_ij|。雅可比迭代在这个条件下必然收敛但它只是充分条件不满足时算法仍然可能收敛。真正的充要条件是迭代矩阵 M I - D^{-1}A 的谱半径 ρ(M) 1谱半径是 M 的特征值模的最大值。判定方法如下判定条件结论类型实际操作严格对角占优充分条件逐行扫描成本 O(n²)ρ(M) 1充要条件对中小矩阵用 np.linalg.eigvalsa_ii 0直接失效需要先行交换或做预处理用np.linalg.eigvals(np.eye(n) - np.linalg.inv(D) A)可以直接算谱半径但矩阵规模一大特征值计算本身就比迭代贵得多所以工程上通常先做对角占优检查再靠迭代过程的残差曲线判断。另一个容易被忽略的问题是a_ii 为零时分母直接失效即便不为零但绝对值很小也容易放大舍入误差此时先做行交换让大对角元上移迭代效果会明显改善。提示不要用残差作为雅可比迭代唯一的停机判据。矩阵病态时残差可能已经很小但解还有明显误差用相邻迭代步的差与残差同时判断更稳妥。5. 雅可比迭代的变体G-S、SOR 与 JOR 的松弛参数怎么调5.1 一行代码从雅可比改到高斯-赛德尔把 jacobi 函数里第 i 行的求和项从“全部用旧值”改成“j i 时用新值j i 时用旧值”就是高斯-赛德尔迭代x_i^{(k1)} (b_i - Σ_{ji} a_ij x_j^{(k1)} - Σ_{ji} a_ij x_j^{(k)}) / a_ii实现上只需把内层循环里的x[j]改成x_next[j]当 j i即可。这一步改动让每个分量能立即使用本轮的更新结果通常迭代次数能比雅可比少一半左右但代价是完全失去了并行性后一个分量的计算依赖前一个分量。对稀疏矩阵高斯-赛德尔的内存访问模式也更好尤其是三对角系统收敛速度差距明显。下面给一个直接可对比的变体def gauss_seidel(A, b, x0None, tol1e-8, max_iter1000): A np.asarray(A, dtypefloat) b np.asarray(b, dtypefloat) n b.size x np.zeros(n) if x0 is None else np.array(x0, dtypefloat) for k in range(max_iter): x_old x.copy() for i in range(n): s A[i] x - A[i, i] * x[i] x[i] (b[i] - s) / A[i, i] if np.max(np.abs(x - x_old)) tol: return x, k 1 raise RuntimeError(f超过最大迭代次数 {max_iter})这里直接原地更新 x再靠 x_old 判断停机代码比 jacobi 更短。需要留意的是高斯-赛德尔的收敛条件与雅可比并不等价存在雅可比收敛而高斯-赛德尔发散的例子。换算法之前先判断 A 是否对称正定对称正定矩阵保证高斯-赛德尔收敛。5.2 SOR 与 JOR松弛因子 ω 的作用和调法雅可比和高斯-赛德尔都对应一个带松弛的推广。JOR 是雅可比加上松弛x^{(k1)} (1-ω) x^{(k)} ω D^{-1}(b (D - A) x^{(k)})SOR 则在高斯-赛德尔迭代式上做同样的加权x_i^{(k1)} (1-ω) x_i^{(k)} ω (b_i - Σ_{ji} a_ij x_j^{(k1)} - Σ_{ji} a_ij x_j^{(k)}) / a_ii两者的 ω 都必须落在 (0, 2) 区间内。ω 1 称为欠松弛用于抑制振荡ω 1 称为超松弛用于加速收敛。对具有相容次序的一类矩阵SOR 的最优松弛因子可以解析地写成 ω_opt 2/(1√(1-ρ(J)²))其中 ρ(J) 是雅可比迭代矩阵 J I - D^{-1}A 的谱半径。但工程中往往不知道 ρ(J)常见做法是从 ω1.2 开始观察残差曲线的下降斜率逐次上调 0.1 试出拐点。下表是几个变体的选型对照变体更新基准并行性典型收敛速度雅可比JORω1全部旧值强慢JORω≠1旧值加松弛强略好于雅可比高斯-赛德尔混合新旧值无约为雅可比两倍SOR混合新值加松弛无最优时远快于 G-S5.3 工程选型什么时候坚持用雅可比而不是 G-S看到这里容易产生一个疑问既然 G-S 和 SOR 收敛更快为什么还要学雅可比答案常在并行场景。雅可比每次迭代所有分量独立计算可以按行切分到不同线程或 GPU block 上几乎不需要线程间通信G-S 每更新一个分量就依赖前一个分量的结果在分布式内存机器上一步同步都省不了。对 GPU 上的巨大稀疏系统雅可比配合 JOR 是目前工程落地的主路线即使迭代次数多单步吞吐量也足以补回来。实际取舍通常遵循这样的经验矩阵规模小、单机跑、收敛困难优先 G-S 或 SOR矩阵规模大、有 GPU、tol 要求不极端优先雅可比加 JOR。还有一类情况必须用雅可比求解器需要完全确定的副作用顺序比如做定点数模拟或硬件仿真时同步更新便于回归测试比对。6. 不精确牛顿迭代法内层用雅可比迭代解线性方程组6.1 双层迭代框架外层牛顿内层雅可比把牛顿迭代法从单变量推广到非线性方程组 F(x)0每一步需要解线性方程组 J(x_k) Δx_k -F(x_k)其中 J(x_k) 是雅可比矩阵Jacobian matrix。注意这里的“雅可比矩阵”和前文的“雅可比迭代算法”是同一个英文词源中文翻译也相同但一个指导数矩阵一个指线性迭代器讨论代码时很容易混淆。外层每走一步内层就可以用第 4 章的 jacobi 函数求解增量 Δximport numpy as np def newton_jacobi(F, J, x0, tol_out1e-10, tol_inner1e-6): x np.array(x0, dtypefloat) for k in range(20): dx, _ jacobi(J(x), -F(x), x0np.zeros_like(x), toltol_inner, max_iter200) x x dx if np.max(np.abs(dx)) tol_out: return x, k 1 raise RuntimeError(外层牛顿迭代未收敛) def F(x): return np.array([x[0]**2 x[1]**2 - 1, x[0] - x[1]]) def Jacobian(x): return np.array([[2*x[0], 2*x[1]], [1.0, -1.0]]) print(newton_jacobi(F, Jacobian, x0np.array([0.5, 0.5])))这个框架下外层只在乎 dx 能不能让解往前走一步内层雅可比迭代的精度是由 tol_inner 控制的。6.2 内外层容差怎么搭配才不浪费算力内层迭代收得越紧外层每步的 dx 越准确但内层迭代步数也越多。经验是外层容差决定最终精度内层容差只需比外层松一到两个数量级即可。把 tol_inner 从 1e-4 收紧到 1e-10外层牛顿步数往往不变变的是每个外层步里的内层迭代步数后者会成倍上升。常见做法是固定 tol_inner 1e-6而不是设置成 tol_out 的几分之一因为外层对 dx 的要求本质上与最终解的精度解耦真正影响外层步数的是内层迭代能否给出一个足够好的下降方向。6.3 验证方法观察内外层步数的变化曲线把上面的 newton_jacobi 稍加修改返回内层平均迭代步数然后对同一初值分别用 tol_inner 1e-4、1e-6、1e-8 各跑一遍。实际观察到的结果通常是外层步数完全相同比如都是 4 步收敛内层平均步数却从个位数涨到几十。这个实验说明不精确牛顿法的瓶颈在内外层容差的配合上而不是内层解方程的精度。碰到外层反复震荡时优先检查 Jacobian 矩阵本身是否准确再考虑收紧内层容差这个顺序几乎总能节省调试时间。本文还有配套的精品资源点击获取