简介面向哈工大硕士《数值分析》课程上机实验这套代码包围绕研究生数值分析常见实验整理涵盖非线性方程组迭代求解、线性代数方程组高斯消元、最小二乘数据拟合和龙贝格数值积分四大模块适合相关专业研究生、高年级本科生及自学者参考。压缩包共8个文件、仅567KB内含4个MATLAB .m源码、2份Word文档、1个Visio流程图和1个README说明.m文件可直接运行分别实现四种核心算法Word文档记录代码思路与实验结果流程图辅助理解控制流。目前已有129人学习下载。读者可对照源码逐步理解数值方法的编程细节也可以借助流程图梳理程序执行顺序并通过实验数据验证结果精度资源允许自行修改适合迁移到课程设计、算法复现或二次开发场景。具体内容包括迭代法逼近非线性方程根、高斯消元求解线性方程组、最小二乘拟合数据曲线以及龙贝格积分近似定积分完整展示了从数学原理到MATLAB实现的转化路径。1. 数值分析上机实验为什么比理论推导更难条件、边界与浮点误差研究生阶段的数值分析课上最难的不是理解教材里的公式推导而是把这些算法写成能跑出合理结果的代码。牛顿法从给定的初值出发直接溢出高斯消去遇到主元接近零被中断插值多项式在端点处剧烈震荡这些现象在板书里几乎看不到却是上机实验最常见的现状。数值分析上机实验要解决的就是从公式到可执行代码之间那条容易被忽略的鸿沟算法能不能收敛、误差怎么控制、参数怎么设、出了问题从哪里查。这份材料面向两类人一类是要交课程报告的学生需要把教材方法跑通并解释结果另一类是在自己工程里复现这些算法的开发者更关心接口怎么设计、换一种矩阵要改哪里。源码和说明书放在一起目的就是让这两类人拿到后都能动手改而不是只当一个黑盒。2. 数值分析上机实验的算法骨架求根、线性方程组、插值与RK4不管课程实验怎么分组有四个主题出现频率最高非线性方程求根、线性方程组求解、插值与数值积分、常微分方程初值问题。这一章按这四个方向给出代码骨架重点不是把算法背一遍而是说明哪些参数必须暴露出来、哪些细节教材没写但上机不处理就会翻车。2.1 非线性方程求根迭代终止条件是数值分析上机实验代码里最先要定的事教材里牛顿法通常写成“重复计算直到前后两次近似值之差小于ε”但只靠这一个条件在上机时会出两类问题函数值已经接近机器精度而位移量还没到容差程序空转初值选得不好时相邻两步差值很小实际却离根越来越远。因此代码里的终止条件要和教材里的收敛性讨论分开处理。我一般会同时检查位移量、函数值和最大迭代次数任何一个条件先触发都返回当前状态而不是让程序死循环。def newton_root(f, df, x0, tol1e-8, max_iter50, reportFalse): 牛顿法求 f(x)0 的单根。 f : 目标函数 df : 导函数 x0 : 初始迭代点 tol : 位移量绝对容差 max_iter : 最大迭代步数 report : 为 True 时返回 (root, iter_count, history) x x0 history [x] for k in range(max_iter): fx f(x) dfx df(x) # 导数接近 0 时继续迭代只会让位移量爆炸 if abs(dfx) 1e-15: raise ZeroDivisionError(fdf({x}) too small) x_new x - fx / dfx # 位移量和函数值同时满足才认为收敛 if abs(x_new - x) tol and abs(f(x_new)) tol: if report: return x_new, k 1, history return x_new, k 1 x x_new history.append(x) raise RuntimeError(fnot converge in {max_iter} iterations)代码里的终止条件被拆成abs(x_new - x) tol与abs(f(x_new)) tol两个判断前者管步长后者管函数取值两个条件同时成立时才返回。max_iter是最后的保险防止初值不合适时无限循环。实际使用时tol要根据实验目的来调只要求交报告1e-6足够要验证高阶收敛性取1e-10以上也合理如果矩阵或函数本身病态把tol放宽到1e-4反而能得到一个可解释的结果。导数接近零时抛异常是故意为之这通常不是程序错误而是初值落在驻点附近换成弦截法或重新选初值即可。2.2 线性方程组列主元消去是数值分析上机实验必须手动实现的改动顺序高斯消去在小规模教学矩阵上表现良好但遇到主元绝对值远小于同列其他元素的矩阵时消元因子会变得很大后续回代解的误差被成倍放大。列主元消去每次把当前列绝对值最大的行换到主元位置一行改动就能把大部分算法错误挡在外面。以下实现不依赖 scipy适合作为实验课源码的起点import numpy as np def gauss_pivot(A, b): 列主元高斯消去求解 Axb。 A : n×n 系数矩阵 b : n 维右端向量 返回解向量 x n A.shape[0] # 组装增广矩阵并确保参与运算的是浮点数 M np.hstack([A.astype(float), b.reshape(-1, 1)]) for col in range(n): # 在当前列内找绝对值最大的元素偏移 col 得到全局行号 piv col np.argmax(np.abs(M[col:, col])) if abs(M[piv, col]) 1e-15: raise ValueError(pivot near zero) if piv ! col: M[[col, piv], :] M[[piv, col], :] for row in range(col 1, n): factor M[row, col] / M[col, col] M[row, col:] - factor * M[col, col:] # 回代求解 x np.zeros(n) for i in range(n - 1, -1, -1): x[i] (M[i, -1] - M[i, i 1:] x[i 1:]) / M[i, i] return xnp.argmax(np.abs(M[col:, col]))返回的是从col开始的相对偏移所以必须加上col才能得到正确的行号。M[[col, piv], :]利用 NumPy 的花式索引做整行交换比临时变量更简洁。回代部分用M[i, i1:] x[i1:]完成已知解与系数行的点积避免再开一层循环。以下表格总结了四种常见实现的选择依据方案主元处理适用情况主要实现成本顺序高斯消去直接用对角元主元有保证的小规模方阵代码最短列主元高斯消去当前列绝对最大值一般实验矩阵多一次 argmax 和行交换全主元消去全局最大值病态严重的矩阵交换逻辑复杂实验较少用LU 分解与列主元结合多次求解不同右端项需要额外存储 L、U2.3 插值与逼近拉格朗日基函数和高次插值的Runge震荡教材上写拉格朗日插值和牛顿插值在数学上等价实验代码里两者却有明显差别。拉格朗日基函数每算一个点都要重新累乘一遍连乘项计算开销大而且当插值节点较多时高次多项式在端点附近的振幅会被非线性放大这就是 Runge 现象。牛顿插值可以递推构造差商表增加一个节点时只需要补一列不必推翻重算因此更适合上机实验里的反复试算。def div_diff_table(x, y): 构造牛顿插值的差商表返回第一行系数。 x : 插值节点 y : 对应函数值 n len(x) coef np.array(y, dtypefloat) # 用 y 作为差商表的第一列 result [coef[0]] for j in range(1, n): # 从当前位置向后更新完成差商的逐层递推 coef[j:] (coef[j:] - coef[j - 1]) / (x[j:] - x[j - 1]) result.append(coef[j]) return result这段代码的核心是coef[j:] (coef[j:] - coef[j-1]) / (x[j:] - x[j-1])。因为每次只用前一个位置的值来更新后面的值所以可以在原数组上就地覆盖不需要额外开辟二维表。result收集每一步对角位置的值最终构成牛顿插值多项式从常数项到最高次项的系数。实验中如果发现高次插值在区间边缘震荡剧烈不要先怀疑代码写错先用三次样条或分段低次插值作对照两者差距能直观说明 Runge 现象。2.4 数值积分与常微分方程RK4单步格式的代码骨架常微分方程的实验通常从显式欧拉开始但欧拉法的每一步误差较大收敛阶只有一阶用来追求精度的上机实验很难看。直接换成经典 RK4代码量增加不大收敛阶从一阶跳到四阶适合作为通用解法。def rk4(f, t0, y0, h, n_steps): 用 RK4 推进初值问题 y f(t, y)。 f : 右端函数签名 f(t, y) t0, y0 : 初值 h : 步长 n_steps : 推进步数 返回 ts, ys ts, ys [t0], [y0] t, y t0, np.array(y0, dtypefloat) for _ in range(n_steps): k1 f(t, y) k2 f(t h / 2, y h / 2 * k1) k3 f(t h / 2, y h / 2 * k2) k4 f(t h, y h * k3) y y h / 6 * (k1 2 * k2 2 * k3 k4) t t h ts.append(t) ys.append(y) return np.array(ts), np.array(ys)每个k对应一种斜率估计四个斜率加权平均后作为这一区间的平均变化率。h / 6 * (...)的权重是 RK4 方法推导出来的常数不要为了省事改成等权重。不同单步法的属性对比如下方便实验报告里直接引用方法局部截断误差每步调用 f 次数全局收敛阶显式欧拉O(h²)11改进欧拉O(h³)22经典 RK4O(h⁵)443. 数值分析上机实验的代码组织源码、数据与说明书怎么放实验 zip 包拿到手第一步不是看代码而是看结构。一份代码如果所有算法都堆在一个文件里参数写死在函数内部改一个容差要全文搜索那“可自己修改”就是一句空话。我经手的课程实验通常拆成四个目录结构清晰了后面改矩阵、换方法都只在固定位置动刀。3.1 实验zip的标准目录src、data、tests、docs以下是我认为比较稳妥的目录组织方式lab/ ├── src/ │ ├── __init__.py │ ├── roots.py # 非线性方程求根 │ ├── linear.py # 线性方程组直接法 │ ├── interp.py # 插值与拟合 │ └── ode.py # 常微分方程初值问题 ├── data/ │ ├── input_A.csv # 系数矩阵 │ ├── input_b.csv # 右端项 │ └── case_config.ini # 实验参数 ├── tests/ │ └── test_all.py # 随机算例与残差校验 ├── docs/ │ └── 实验说明书.md └── run_lab.py # 主入口脚本src目录只放算法函数每个文件对应一类问题data放输入数据无论是矩阵还是生成的随机数都通过文件读入tests放测试脚本保证改完代码后还能用旧算例回归docs放说明书记录公式来源、变量含义和已知问题。run_lab.py是唯一允许出现print和文件读写的脚本这样算法函数保持纯净换个实验题目时不用动src下的代码。3.2 算法模块与主流程分离一份源码对应多个算例很多实验代码把读文件、调算法、打印结果写在同一个函数里看起来方便实际上每换一组数据都要改主流程。常见的做法是让算法函数只接收 NumPy 数组并返回结果IO 全部放到主入口里import argparse import numpy as np from src.linear import gauss_pivot def load_matrix(path): 从 CSV 读取矩阵自动识别分隔符 return np.loadtxt(path, delimiter,) def main(): ap argparse.ArgumentParser(descriptionnumerical lab runner) ap.add_argument(--a, defaultdata/input_A.csv) ap.add_argument(--b, defaultdata/input_b.csv) args ap.parse_args() A load_matrix(args.a) b load_matrix(args.b).reshape(-1, 1) if A.shape[0] ! b.shape[0]: raise SystemExit(row mismatch between A and b) x gauss_pivot(A, b) # 算法只接收数组 np.savetxt(output_solution.csv, x, delimiter,) print(done) if __name__ __main__: main()--a和--b两个命令行参数允许你指定不同的矩阵文件跑新算例时就无需打开代码改了。reshape(-1, 1)把 b 从一维数组变成列向量保证与增广矩阵拼接时维度对齐。主入口里做维度校验是上机实验最容易漏但最值得加的一步——很多求解器报维度错误根本原因是 b 的形状不对。3.3 参数外置用配置文件代替硬编码tol、max_iter、method这类参数属于运行配置不应该出现在算法函数定义里。把参数写进配置文件后实验者不需要理解代码细节也能调整import configparser config configparser.ConfigParser() config.read(data/case_config.ini) tol float(config[hyper][tol]) max_iter int(config[hyper][max_iter]) method config[hyper][method]对应的case_config.ini内容很简单[hyper] tol 1e-8 max_iter 50 method gauss-pivot实验报告里通常要说明不同容差对结果的影响配置文件的作用就是把这一块隔离出来让你可以在报告里直接列出多组tol的运行结果。另一个好处是多人协同做实验时大家不需要互相改动源码只要传配置文件和输入数据文件即可。3.4 说明书与源码的对应关系拿到包后怎么核对完整性说明书不应该只是把源码贴一遍而是要回答三个问题这个函数实现的是哪个公式、输入的每列数据代表什么、算出的结果在什么条件下可信。我看一份实验说明书重点看它有没有记录“已知失败条件”比如主元为零、迭代发散、矩阵条件数过大。这类信息是教科书里没有的恰恰是上机实验的价值所在。提示如果 hand 内没有tests目录务必先手动构造一个已知解的问题比如x_true np.ones(n)再让b A x_true用回代结果与x_true的差值验证算法正确性再开始改代码。4. 修改数值分析上机实验代码的三个高频场景精度、病态矩阵与新算法标题里写着“可自己修改”但修改不是乱改。课程实验里最常见的改动就三类换精度、换测试矩阵、换求解方法。每类都有固定的改法改错位置会引入隐藏问题。4.1 换精度float换成double以及Python浮点输入的坑C 语言里把float换成double通常会改变主元判断的阈值Python 里则需要显式做类型转换。NumPy 读取 CSV 时如果矩阵元素全是整数默认得到int64数组直接参与除法会丢失小数位。下面的代码是安全的起点A A.astype(np.float64) # 强制提升为双精度 b b.astype(np.float64) # 根据实际精度选择容差 tol 1e-10 if A.dtype np.float64 else 1e-5把数组显式转换成float64后后续所有的除法运算都会按浮点数执行。C/C 环境下还要注意 printf 格式串float传入printf会被自动提升为double所以调试打印时看不出问题只有参与大量累加运算才会暴露精度差异。换精度后第一件事是跑一遍残差确认误差量级随精度的变化符合预期。4.2 换测试矩阵随机矩阵换成Hilbert矩阵很多同学跑完随机矩阵就认为代码没错实际上随机矩阵通常条件数很小掩盖了算法在病态矩阵上的缺陷。Hilbert 矩阵是最容易构造的极端测试用例元素是H[i][j] 1 / (i j 1)阶数越高越病态import numpy as np def hilbert(n): 生成 n×n Hilbert 矩阵 i, j np.indices((n, n)) return 1.0 / (i j 1.0) n 8 H hilbert(n) x_true np.ones(n) b H x_true # 已知真解的右端项 print(cond :, np.linalg.cond(H))用这个脚本生成右侧项后你再把H和b喂给第 2.2 节的gauss_pivot会发现解的前几位可能还正确后面几位已经开始跳动。这是病态矩阵的固有性质不是高斯消去实现错了。实验报告里如果用 Hilbert 矩阵作为算例必须同时给出矩阵条件数说明残差小不等于解误差小。矩阵类型条件数量级特征均匀随机矩阵1~10²求解容易适合快速冒烟测试对角占优矩阵1~10几乎所有直接法都能正确求解Hilbert 矩阵 n8约 10¹⁰双精度下失掉大半有效数字接近奇异矩阵随时间可控调整需要结合残差与条件数一起评估4.3 新增算法在方法注册表里加一条而不是改主流程新增算法时最容易犯的错误是把新方法写进main变成if name custom的分支这样每加一个方法就要动主流程一次。我一般用一个字典做方法注册表# src/solvers.py from .linear import gauss_pivot, solve_lu def solve_cg(A, b, tol1e-8, max_iter1000): 共轭梯度法占位实现返回解向量 ... METHODS { gauss: gauss_pivot, lu: solve_lu, cg: solve_cg, } def solve(A, b, methodgauss): if method not in METHODS: raise KeyError(funknown method: {method}) return METHODS[method](A, b)新增方法时只需要满足一个约定即输入A和b、输出解向量然后在METHODS里注册一行。主入口里的--method cg就能直接切到新算法。这个方法对课程实验的意义在于同一个实验报告可以对比多个方法的残差与耗时而不用复制粘贴多份互不相干的主程序。5. 验证数值分析上机实验代码的三个硬指标残差、收敛阶与打开包后的第一件事5.1 先算残差再下结论不要盯着解的最后几位数字实验报告里最经典的错误是打印出来的解和书中答案前几位一致就直接写“结果正确”。正确做法是计算残差也就是把解代回原方程后的偏差# 解 x 是否可信先看残差的最大分量 residual np.linalg.norm(A x - b, ordnp.inf) assert residual 1e-8, fresidual too large: {residual}ordnp.inf表示取向量各分量的最大值选它是因为它最符合“每个方程都不能偏差太多”的直觉。残差小于容差说明代码实现没有明显错误残差大则说明消去过程、列主元或回代部分需要逐行检查。但如果矩阵本身病态残差小也不能完全证明解可信这时需要结合条件数一起判断。5.2 用收敛阶定位实现错误二阶方法只跑出一阶时的排查点只验证一个步长下的误差不够真正的算法错误要靠收敛阶才能暴露。方法如下步长依次减半记录全局误差再计算相邻误差的比值取以 2 为底的对数import numpy as np # err 是不同步长下的全局误差 err np.array([3.1e-4, 7.8e-5, 1.9e-5, 4.8e-6]) order np.log2(err[:-1] / err[1:]) print(observed order:, order)如果对 RK4 算出的order接近 1 而不是 4优先检查三处时间推进时是否把t和y同时更新了、初值是否被重新赋值、步长是否真的每次都减半而不是固定步长。收敛阶这个指标对实现错误十分敏感小到把h/2写成了h都会反映在阶数上。5.3 打开zip之后的第一件事确认源码、说明书、测试脚本三件套拿到实验 zip 包时不急着解压运行先看清单unzip -l numerical_lab.zip # 列出压缩包内容 unzip -o numerical_lab.zip -d lab # 解压到独立目录 cd lab ls docs src tests # 三部分缺一不可 python -m pytest tests -q # 先跑自带测试unzip -l能在解压之前判断包里是否同时包含源码、数据文件和说明书避免把恶意脚本或损坏文件直接释放到当前目录。解压后第一件事是跑测试测试通过说明环境依赖基本齐全测试失败则先看是算法问题还是路径问题例如 CSV 路径写死就会导致所有用例失败。把这三个检查点固定成习惯下次无论是拿别人的包还是发自己的包都能在五分钟内判断能否继续改。本文还有配套的精品资源点击获取