简介广东工业大学数值分析课程实验四的Python实验报告适合正在学习数值分析、科学计算或Python编程的高校学生使用。资源围绕三个典型计算问题展开超定方程组的最小二乘解、给定观察数据的最小二乘拟合、以及石油产量预测模型的建立。其中题目三分别求解两组具体超定方程组的最小二乘解完整推导了正规方程组ATAxATb并调用numpy库完成计算题目六基于给定观察数据采用形如yc0c1/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 为什么“看起来收敛了”不等于“算对了”这是实验四最容易让人翻车的认知盲区。一个算法收敛指的是随着步长趋于零数值结果趋于某个极限值但这个极限值不一定是真值——如果公式本身推导有误或者实现里漏了系数你会观察到结果逐渐稳定在一个错误的数上还误以为“节点多了就是准”。我见过一份实验报告复合梯形公式从 n10 到 n10^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(fSimpson result: {result:.10f})这个实现的关键在于weights[1:-1:2] 4和weights[2:-2:2] 2两行——奇数下标对应每个小区间中点的函数值系数为 4偶数下标对应相邻小区间的共享端点除了全局端点外系数为 2。n必须为偶数因为辛普森公式要求每个子区间由两个小区间组成。h/3是公式的公共系数挪到sum外面统一乘避免在循环里重复计算。参数修改上n100对大多数光滑函数已经足够实验报告如果想画“误差随 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 (ab)/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(fn{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(fGauss-{k} points err{abs(gauss - true_val):.2e})代码里值得注意的一行是w (b - a) / 2 * w_ref——这是区间变换的雅可比因子。很多同学抄高斯求积代码时会漏掉这个因子导致结果永远差(b-a)/2倍。另外leggauss返回的节点是从小到大排列的但求和时节点顺序不影响结果所以不需要额外排序。跑完这段代码你会发现一个很有趣的现象复合梯形 n80 的误差大概在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 # 注意: 这个函数在 x0 处需要特殊处理 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(fn{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]上积分复合辛普森 n20 的误差就能到10^-9再往上加密网格纯属浪费机器时间。反过来如果被积函数在区间内有尖锐变化比如sqrt(x)在接近 0 的地方、或者1/(1100x^2)这类有“山脊”的函数n 要加到几百甚至上千误差才肯降下来。这时候观察误差随 n 的变化曲线如果出现“平台期”——误差不再随 n 增加而显著下降——那不是算法问题而是浮点舍入误差开始主导了。复合辛普森的舍入误差大约正比于n * eps * |f|eps是双精度浮点的机器精度约2.2e-16所以 n 加到10^6级别时舍入误差会盖过截断误差再多算也没用。一个务实的做法先取 n10、20、40、80、160 跑一组看误差下降趋势。如果误差下降规律符合理论阶数Simpson 是 4那 n 取到误差小于10^-8就够了如果下降趋势不符合多半是代码写错了而不是 n 不够大。这也是上面收敛阶分析最大的价值——它给你一个“算法有没有正常工作的判据”。4.2 积分区间 [a, b]端点奇异和“看起来正常”的陷阱实验四给的积分区间往往很干净比如[0, 1]、[0, π]但换一个类似的形式就可能有坑。最典型的坑是端点为奇异点对sin(x)/x在[0, π]积分数学上真值存在但计算机在 x0 处会得到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(2t1) 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, a2.0: np.exp(a * x)或者定义闭包。我见过同学在循环里改了全局a但f捕获的是旧值导致误差分析表里一堆莫名其妙的数字。这种问题有个快速排查法——打印f(0)和f(1)手动验算一下一分钟就能定位。5. 实验四避坑指南四个高发翻车点与排查方法5.1 现象n 加大后误差不降反升甚至出现负的收敛阶这是实验四报告里最让人抓狂的情况理论上 n 越大越准实际跑出来 n160 的误差比 n80 还大。原因通常有两类——第一是被积函数本身有奇异点你没有处理第二是浮点舍入误差开始主导前面 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 或返回 NaNsin(x)/x这类函数在 0 处要单独处理。常见做法是写一个分段定义的函数在奇点附近用极限值替代。比如sin(x)/x在 0 处的极限是 1可以写成np.where(x 0, 1.0, np.sin(x)/x)。但要注意——如果x刚好在数值上不为 0 但非常接近 0比如1e-16np.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% 的“看着对但实际错”的结果。实验四这个作业说实话过了就过了但数值计算的基本功不过关后面做优化、做仿真、做数据拟合都会还债。至少我出去做项目时遇到过不止一次同事拿数值积分结果当真理最后发现是端点奇异没处理——这类血泪教训完全可以早点避开。希望这篇笔记能帮你把实验四做得又快又稳也希望你养成那个“先验证再相信”的习惯它比任何一次实验分数都值钱。本文还有配套的精品资源点击获取