简介：广东工业大学数值分析课程实验四的Python实验报告，适合正在学习数值分析、科学计算或Python编程的高校学生使用。资源围绕三个典型计算问题展开：超定方程组的最小二乘解、给定观察数据的最小二乘拟合、以及石油产量预测模型的建立。其中题目三分别求解两组具体超定方程组的最小二乘解，完整推导了正规方程组ATAx=ATb并调用numpy库完成计算；题目六基于给定观察数据，采用形如y=c0+c1/x的线性组合模型进行拟合，并利用线性回归工具输出参数；实验四则进一步将最小二乘方法用于石油产量预测，体现数值分析在实际数据场景中的应用。压缩包内为1个docx文档，大小仅112KB，但内容覆盖算法原理、源代码与结果分析。已有212人学习下载，尤其适合需要完成同类上机实验、撰写数值分析报告或复习最小二乘相关知识的读者，有助于快速掌握Python在数值计算中的实际应用。
1. 数值分析实验四：算法选对了，结果才不是玄学
数值分析实验四在很多高校的 Python 课程里都是分水岭——前三周还在写循环、调列表，到这一周突然要面对误差阶、收敛性、截断误差这些真正属于“数值方法”的概念。我见过不少同学把实验四做成“调库大会”：quad一行出结果，图表一画就交差，但问一句“为什么用复合辛普森而不是梯形”就答不上来。这篇笔记想说的就是：实验四表面上考的是“把定积分算出来”，实际上考的是“你知不知道自己在算什么”。
这套实验真正练的是三件事：第一，理解不同求积公式的代数精度，知道梯形、辛普森、高斯求积各自能吃到多少精度；第二，学会用误差阶判断一个算法“有没有在正常工作”，而不是只看结果像不像；第三，把数学公式翻译成可维护的 Python 代码，并养成用解析解或高精度参考值验证的习惯。适合正在被实验四折磨的同学，也适合想补一补“数值计算为什么不能盲信机器结果”的从业者。接下来直接进入正题。
2. 数值积分的算法选型：三种求积公式，各自能吃多少精度
2.1 从矩形到辛普森：代数精度是怎么一步步涨上去的
数值积分的基本思路很简单——用有限个点上的函数值去逼近区间上的积分。最粗糙的矩形公式把被积函数看成常数，误差正比于步长的一次方，记作O(h)；梯形公式把相邻两点连成直线，误差降到O(h^2)；辛普森公式用三点拟合一条抛物线，误差直接跳到O(h^4)。这个“误差阶每升一级，步长减半时误差缩小的倍数就翻一档”的规律，是整个实验四的核心。
代数精度这个概念值得单独拎出来说。一个求积公式的代数精度是“它能被精确满足的最高次多项式次数”：矩形公式对常数精确，精度是 0；梯形公式对一次多项式精确，精度是 1；辛普森公式因为有 3 个自由参数，对二次和三次多项式都精确，精度是 3。别小看这个数字——它决定了当你把一个光滑函数拆成泰勒展开后，公式到底能吃住前面几项。吃住的项越多，截断误差越小。
实验四里最常见的一个误区是“节点越多结果越准”。这在理论上成立，但前提是你用的公式本身精度够高。用梯形公式硬算 10 万个子区间，误差可能还不如辛普森公式算 100 个子区间。原因就是梯形公式的误差阶是O(h^2)，而辛普森是O(h^4)——步长缩小 10 倍，前者误差缩小 100 倍，后者缩小 10000 倍。这就是为什么算法选型比盲目加密网格更重要。
2.2 高斯求积：为什么它用很少的节点就能达到很高的精度
高斯求积是另一条路线：与其均分节点，不如把节点位置本身也当作未知数来优化。n 点高斯-勒让德公式有 2n 个自由参数（n 个位置加 n 个权重），所以代数精度能做到 2n-1。两点高斯求积的代数精度就是 3，和三个节点的辛普森公式持平，但只用了两个函数值——这就是“效率”二字的含义。
代价也很明显：高斯求积的节点不是等距的，而且节点位置依赖积分区间。实际做法是先通过线性变换把任意区间[a, b]映射到[-1, 1]，然后查表或用numpy的leggauss拿到节点和权重。实验四如果要求手写高斯求积，难点往往不在公式本身，而在区间变换时忘了乘雅可比因子(b-a)/2——这个系数漏了，结果差一倍起步。
选型建议我给三条，都是我实际带实验时总结的：被积函数光滑且不知道任何先验信息，优先复合辛普森，代码短、行为稳、误差阶够写报告；被积函数在区间端点有奇异性，改用高斯求积或分段处理，别硬用辛普森；被积函数是震荡函数（比如sin(100x)），等距节点公式全部失效，这时候该上自适应积分或变量替换，这个在第 6 章会展开。
2.3 为什么“看起来收敛了”不等于“算对了”
这是实验四最容易让人翻车的认知盲区。一个算法收敛，指的是随着步长趋于零，数值结果趋于某个极限值；但这个极限值不一定是真值——如果公式本身推导有误，或者实现里漏了系数，你会观察到结果逐渐稳定在一个错误的数上，还误以为“节点多了就是准”。
我见过一份实验报告，复合梯形公式从 n=10 到 n=10^6 都稳定在 1.5707 附近，学生以为这就是 π/2。实际检查发现他写的被积函数把x^2写成了x，程序一直在精确地计算另一个积分。这件事给实验四提了个醒：收敛性只能说明算法自洽，不能说明结果正确。正确性要靠外部参照来验证——解析解、高精度库函数、或者不同算法之间的交叉验证，三者至少占一样。这也是为什么第 3 章我坚持在代码里加一个scipy.integrate.quad的高精度结果做对照。
3. 用 Python 把实验四跑通：从被积函数到误差分析表的完整代码
3.1 最小可复现版本：手写复合辛普森求积
先给一个可以直接跑通的最小实现。这个版本刻意不依赖scipy，只用numpy，方便你在实验报告里展示“手写算法”的过程。
import numpy as np def composite_simpson(f, a, b, n): """ 复合辛普森求积公式 参数: f: 被积函数，接受标量或数组输入 a, b: 积分区间 [a, b] n: 子区间数量，必须是偶数 返回: 积分近似值 """ if n % 2 != 0: raise ValueError("n 必须为偶数") h = (b - a) / n x = np.linspace(a, b, n + 1) y = f(x) # 辛普森权重: 端点系数为1, 奇数下标系数为4, 偶数下标(非端点)系数为2 weights = np.ones(n + 1) weights[1:-1:2] = 4 weights[2:-2:2] = 2 return h / 3 * np.sum(weights * y) # 测试: 计算 ∫_0^π sin(x) dx，真值应为 2 result = composite_simpson(np.sin, 0, np.pi, 100) print(f"Simpson result: {result:.10f}")这个实现的关键在于weights[1:-1:2] = 4和weights[2:-2:2] = 2两行——奇数下标对应每个小区间中点的函数值，系数为 4；偶数下标对应相邻小区间的共享端点，除了全局端点外系数为 2。n必须为偶数，因为辛普森公式要求每个子区间由两个小区间组成。h/3是公式的公共系数，挪到sum外面统一乘，避免在循环里重复计算。
参数修改上，n=100对大多数光滑函数已经足够，实验报告如果想画“误差随 n 的变化”，把n取[10, 20, 40, 80, 160]就够看趋势。如果被积函数在区间内变化剧烈，n提到 1000 也不过分，但超过 10000 后计算时间的增长会明显变快，这时应该考虑自适应积分而不是继续加密网格。
3.2 梯形与高斯求积：同一份接口，三种算法横向对比
实验四如果只做一种算法，报告会显得单薄。我一般建议把三种算法放进同一个函数签名格式里，横向对比才方便。下面的代码实现了复合梯形和两点/三点高斯求积，并打印误差对比表。
import numpy as np from numpy.polynomial.legendre import leggauss def composite_trapezoid(f, a, b, n): """复合梯形公式: 误差阶 O(h^2)""" h = (b - a) / n x = np.linspace(a, b, n + 1) y = f(x) # 梯形权重: 端点为1, 内部点为2, 总系数为 h/2 weights = np.ones(n + 1) weights[1:-1] = 2 return h / 2 * np.sum(weights * y) def gauss_legendre(f, a, b, n_points): """ n_points 点高斯-勒让德求积 节点和权重由 leggauss 给出，定义在 [-1, 1] """ x_ref, w_ref = leggauss(n_points) # 区间变换: [a, b] -> [-1, 1], 对应公式 x = (b-a)/2 * t + (a+b)/2 x = (b - a) / 2 * x_ref + (a + b) / 2 w = (b - a) / 2 * w_ref return np.sum(w * f(x)) # 测试函数: ∫_0^1 e^x dx 真值 = e - 1 f = lambda x: np.exp(x) true_val = np.e - 1 for n in [10, 20, 40, 80]: simpson = composite_simpson(f, 0, 1, n) trapezoid = composite_trapezoid(f, 0, 1, n) print(f"n={n:4d} Simpson err={abs(simpson - true_val):.2e} Trapezoid err={abs(trapezoid - true_val):.2e}") # 两点/三点高斯求积 for k in [2, 3]: gauss = gauss_legendre(f, 0, 1, k) print(f"Gauss-{k} points err={abs(gauss - true_val):.2e}")代码里值得注意的一行是w = (b - a) / 2 * w_ref——这是区间变换的雅可比因子。很多同学抄高斯求积代码时会漏掉这个因子，导致结果永远差(b-a)/2倍。另外leggauss返回的节点是从小到大排列的，但求和时节点顺序不影响结果，所以不需要额外排序。
跑完这段代码你会发现一个很有趣的现象：复合梯形 n=80 的误差大概在10^-4量级，三点高斯求积误差直接到10^-8量级。这就是代数精度差异带来的实际效果——高斯求积只用了 3 个点，精度吊打等距节点 81 个点。报告里如果能把这张对比表放进去，实验四的深度立马上一个台阶。
3.3 误差分析表：打印收敛阶才是实验报告的点睛之笔
实验报告最容易被扣分的部分，是“只有结果数字，没有收敛阶分析”。收敛阶的数值验证很简单：如果你有两组步长h1和h2对应的误差e1和e2，那么收敛阶 p 近似等于log(e1/e2) / log(h1/h2)。代码实现如下。
def convergence_order(errors, ns): """ 根据相邻两组误差估计收敛阶 errors: 误差列表, 从粗网格到细网格排列 ns: 对应的子区间数量列表 """ orders = [] for i in range(len(ns) - 1): # 步长与 n 成反比, 因此 log(h1/h2) = log(n2/n1) p = np.log(errors[i] / errors[i + 1]) / np.log(ns[i + 1] / ns[i]) orders.append(p) return orders ns = [10, 20, 40, 80, 160] f = lambda x: np.sin(x) / x # 注意: 这个函数在 x=0 处需要特殊处理 a, b = 1e-8, np.pi # 端点避开了 0 errs_simpson = [] for n in ns: val = composite_simpson(f, a, b, n) true_val, _ = integrate.quad(f, a, b) # 需要 import scipy.integrate as integrate errs_simpson.append(abs(val - true_val)) orders = convergence_order(errs_simpson, ns) for n, e, p in zip(ns, errs_simpson, orders): print(f"n={n:4d} err={e:.2e} order={p:.2f}")理论上复合辛普森的收敛阶应该逼近 4，复合梯形逼近 2。如果你跑出来的收敛阶在 3.5 到 4.2 之间，说明算法实现和误差分析都正确；如果收敛阶在 2 附近，说明你实际用的是梯形公式；如果收敛阶忽高忽低或者甚至为负，那基本可以断定被积函数或区间设置有问题，先去排查边界。这段代码是我个人强烈建议写进实验报告的——它让老师一眼看出你不是在“跑程序”，而是在“验证方法”。
4. 步长、节点与误差：三个必调参数，以及实验结果的解读门道
4.1 子区间数量 n：选多少才能让误差低于“报告门槛”
实验四里n怎么选，是所有参数里最影响结果的一个。我常用的判断标准不是“n 越大越好”，而是“让误差降到你实验方法本身的精度瓶颈之下”。如果你的被积函数和区间都相对普通，比如exp(x)在[0,1]上积分，复合辛普森 n=20 的误差就能到10^-9，再往上加密网格纯属浪费机器时间。
反过来，如果被积函数在区间内有尖锐变化，比如sqrt(x)在接近 0 的地方、或者1/(1+100x^2)这类有“山脊”的函数，n 要加到几百甚至上千误差才肯降下来。这时候观察误差随 n 的变化曲线，如果出现“平台期”——误差不再随 n 增加而显著下降——那不是算法问题，而是浮点舍入误差开始主导了。复合辛普森的舍入误差大约正比于n * eps * |f|，eps是双精度浮点的机器精度约2.2e-16，所以 n 加到10^6级别时，舍入误差会盖过截断误差，再多算也没用。
一个务实的做法：先取 n=10、20、40、80、160 跑一组，看误差下降趋势。如果误差下降规律符合理论阶数（Simpson 是 4），那 n 取到误差小于10^-8就够了；如果下降趋势不符合，多半是代码写错了，而不是 n 不够大。这也是上面收敛阶分析最大的价值——它给你一个“算法有没有正常工作的判据”。
4.2 积分区间 [a, b]：端点奇异和“看起来正常”的陷阱
实验四给的积分区间往往很干净，比如[0, 1]、[0, π]，但换一个类似的形式就可能有坑。最典型的坑是端点为奇异点：对sin(x)/x在[0, π]积分，数学上真值存在，但计算机在 x=0 处会得到0/0的 NaN。这种函数在物理和工程里很常见，处理办法通常是避开端点——积分从1e-8开始，而不是从 0 开始。前面 3.3 节代码已经演示了这个写法。
另一个隐蔽陷阱是区间太大导致的高频震荡。比如∫_0^100 sin(x^2) dx，这类函数即使光滑且无奇异点，等距节点求积也得用非常多的节点才能跟上震荡节奏。更合理的方案是用高斯求积搭配大型 n（比如 50 点），或者做变量替换把震荡“压平”。实验四如果遇到这种函数，第一步不是调 n，而是停下来想想：这个函数是不是换一种表示方法更好积分。比如x = sqrt(t)的替换，可能把震荡频率降低一个量级。
区间参数还有一个容易被忽略的细节——积分区间的线性变换会影响精度。对同一个积分∫_0^2 f(x) dx，写成∫_{-1}^1 2*f(2t+1) dt后，公式的舍入误差表现会不一样。高精度计算时，我一般会尽量把区间平移和缩放做到“函数值不过大也不过小”，避免中间结果溢出或下溢。实验四阶段不用太追求这个，但报告中如果能提到“将区间映射到 [-1,1] 以避免函数值尺度问题”，会显得很有思考深度。
4.3 被积函数的写法：传入方式决定你能不能横向对比
被积函数用def还是lambda、能不能接收数组输入，这些看起来是 Python 语法的小事，实际上直接影响算法函数的通用性。上面的代码里我用f(x)直接传入数组，这样numpy的向量化运算可以一次性算出全部函数值，比循环快两个量级。如果被积函数里用了math.sin而不是np.sin，传入数组时会直接报错——这是实验四最常见的 TypeError 来源之一。
另一个容易踩的坑是函数内部的全局变量。比如你要对不同的参数a计算∫_0^1 e^(a*x) dx，如果a定义在函数外面，在循环里改a时必须保证函数能读到新的值。更稳的做法是把参数写进函数签名，比如f = lambda x, a=2.0: np.exp(a * x)，或者定义闭包。我见过同学在循环里改了全局a，但f捕获的是旧值，导致误差分析表里一堆莫名其妙的数字。这种问题有个快速排查法——打印f(0)和f(1)手动验算一下，一分钟就能定位。
5. 实验四避坑指南：四个高发翻车点与排查方法
5.1 现象：n 加大后误差不降反升，甚至出现负的收敛阶
这是实验四报告里最让人抓狂的情况：理论上 n 越大越准，实际跑出来 n=160 的误差比 n=80 还大。原因通常有两类——第一是被积函数本身有奇异点，你没有处理；第二是浮点舍入误差开始主导，前面 4.1 节提过，复合公式的舍入误差随 n 增大而累积。解决方法是先画出被积函数在区间内的图像，确认函数值范围和变化形态；如果函数值跨越多个数量级，考虑换区间变换或改用高斯求积。另一个排查点是你有没有把常数系数写对——h/3写成了h*3，误差会差 9 倍。
5.2 现象：高斯求积和梯形、辛普森的结果对不上
三种算法的结果差了10^-4以上时，不要急着怀疑高斯求积“更高级所以是对的”——先检查高斯求积的权重和节点从哪拿的。常见问题包括：leggauss返回的权重没有乘以(b-a)/2；或者你用的节点是[-1, 1]上的，但没有做线性变换就把被积函数原封不动地代进去。第二个问题是很多人会犯的：节点定义域不同，函数输入必须跟着变换。排查方法很简单——对常数函数f(x)=1测试积分，结果应该是b-a。如果这个测试都不通过，说明区间变换写错了，与算法精度无关。
5.3 现象：被积函数在端点处能算，但程序直接报 RuntimeWarning 或返回 NaN
sin(x)/x这类函数在 0 处要单独处理。常见做法是写一个分段定义的函数：在奇点附近用极限值替代。比如sin(x)/x在 0 处的极限是 1，可以写成np.where(x == 0, 1.0, np.sin(x)/x)。但要注意——如果x刚好在数值上不为 0 但非常接近 0，比如1e-16，np.sin(x)/x仍然能算出 1.0，不会出问题。真正的风险在解析解和数值解对比时：如果你的“真值”也用了同样的奇异点处理方式，两边会共享同一个系统性偏差，误差看起来很小但实际是错的。所以奇异函数最好用级数展开或者查高精度参考值做对照。
5.4 现象：三种算法结果很接近，但和解析解差了一个常数
这是最隐蔽的一类 bug——某个系数写错了，但被三种算法一致地继承下来。比如复合辛普森公式的h/3，如果你写成h/6，算法依然是收敛的，误差阶依然是 4，但收敛到的极限值是真实积分的一半。所有横向对比都会显示“三种算法互相吻合”，唯独看起来无比正确的解析解对不上。排查思路只有一个：拿一个你知道解析解的简单函数做基准测试，比如∫_0^1 x^2 dx = 1/3、∫_0^π sin(x) dx = 2，逐项验证你的实现。这个习惯在实验四做完之后依然重要——工作里用任何数值方法，第一步永远是拿“已知答案的问题”去验证实现，而不是直接信任第一次跑出来的结果。
6. 从“过实验”到“真会用”：自适应积分与解析解对照的验证习惯
实验四如果只满足于手写三种公式并对比误差，那只是把书本走了一遍。真正把数值积分用出价值的人，会进一步接触自适应积分——让程序自动判断哪里该加密网格，哪里可以稀疏。这个思想在工程里的意义是：同样达到10^-8的误差，自适应算法可能只要几百次函数求值，而均匀网格需要几万次。简化版的思路是这样：先对整个区间做一个辛普森估计，再对半分两个子区间分别估计，如果两半的估计值之和与整体估计值的差足够小，就接受；否则对两半各自递归。
def adaptive_simpson(f, a, b, tol): """简化版自适应辛普森积分""" c = (a + b) / 2 whole = composite_simpson(f, a, b, 2) # 整体用 2 个子区间 left = composite_simpson(f, a, c, 2) right = composite_simpson(f, c, b, 2) if abs(left + right - whole) < 15 * tol: # 误差足够小, 接受 — 15 是辛普森公式的经验安全系数 return left + right + (left + right - whole) / 15 h = (b - a) / 2 return adaptive_simpson(f, a, c, tol/2) + adaptive_simpson(f, c, b, tol/2)这个版本没有处理递归深度问题，但用在实验四的氛围函数上已经够用。跑一下adaptive_simpson(np.sin, 0, np.pi, 1e-10)，你会发现函数求值次数远少于均匀网格，精度却反超。这就是数值方法“用思考换计算量”的典型例子。
我自己的习惯是：所有数值计算结果，先问自己三个问题。第一，有没有解析解可以对照？有就用解析解验证；没有就找一个高精度的参考值，比如scipy.integrate.quad容差调到1e-14的结果。第二，不同方法的结果是否一致？一致不代表对，但不一致一定哪里有问题。第三，收敛阶是否符合理论预期？符合预期说明算法实现正确，不符合就回到第 5 章的排查流程。这三个问题花不了两分钟，但能拦住 90% 的“看着对但实际错”的结果。
实验四这个作业，说实话，过了就过了，但数值计算的基本功不过关，后面做优化、做仿真、做数据拟合都会还债。至少我出去做项目时，遇到过不止一次同事拿数值积分结果当真理，最后发现是端点奇异没处理——这类血泪教训完全可以早点避开。希望这篇笔记能帮你把实验四做得又快又稳，也希望你养成那个“先验证再相信”的习惯，它比任何一次实验分数都值钱。
本文还有配套的精品资源，点击获取
