简介针对广东工业大学python数值分析课程实验四的配套资料内容聚焦超定方程组最小二乘解与最小二乘拟合并延伸至石油产量预测应用。报告从算法原理讲起结合正规方程组ATAxATb推导给出具体求解步骤和结果分析附有基于numpy、sklearn、matplotlib的完整Python源代码适合正在学习数值分析或需要参考上机实验报告的学生使用。资源包为单一docx文档整体大小112KB无需额外解压配置。已有212人浏览/学习对于课程实验类资料具有一定参考价值。文档中题目三演示了超定方程组a、b的最小二乘解计算题目六给出了观察数据的最小二乘拟合及函数图像绘制实验四则通过石油产量数据展示预测模型建立与参数优化。读者可从中获取实验思路、代码模板和报告排版参考尤其能加深对最小二乘方法从理论到编程实现的理解。1. 数值分析实验四Python 手写数值积分先想清楚要交什么数值分析实验四在很多高校的课程安排里就是数值积分与龙贝格算法这一关给几个被积函数让你用 Python 手写复化梯形、复化辛普森、龙贝格再算精度、分析误差阶最后交一份带图表和结论的实验报告。听起来是套公式的活儿真正动手你会发现公式抄得再对程序跑出来的数字也可能差到离谱。辛普森在光滑函数上收敛得规规矩矩换一个带振荡的被积函数立刻翻车龙贝格表看着收敛遇到不连续点就给出一个假收敛的漂亮数字。这篇笔记把数值积分实现里真正值钱的部分拆开讲方法怎么选、代码怎么写、参数怎么定、哪些坑必须躲以及怎么让报告里的误差分析有说服力。正在赶实验的同学可以直接照抄代码想把手写积分用在工程里的熟手也能看到边界条件。2. 数值积分方法选型梯形、辛普森与龙贝格各自的误差底牌数值积分的核心思路只有一个用容易求积分的函数去逼近那个难搞的被积函数。最常用的逼近对象是插值多项式——在区间上取几个节点构造低次多项式积多项式的解析结果再把节点之差的幂次带进去。这套思路衍生出的 Newton-Cotes 家族就是梯形、辛普森这些公式的源头。选哪个公式本质上是在选多项式次数和误差的衰减速度。2.1 复化梯形与复化辛普森步长减半时误差按什么速度衰减先说梯形公式。用一次多项式连接区间两端点面积就是(b-a)/2 * (f(a) f(b))。误差正比于h²系数是-(b-a)/12 * f(ξ)。注意这里的 h 是整个区间长度。如果区间很大或者函数弯曲得厉害直接套这个公式误差完全不可控所以必须复化把[a,b]均匀切成 n 段每段用梯形再求和。复化梯形公式T_n h/2 * [f(x0) 2*Σf(xi) f(xn)]i 从 1 到 n-1h (b-a)/n关键结论是h 减半n 翻倍误差变成原来的 1/4。因为整体误差是O(h²)这是二阶收敛。复化辛普森在此基础上再进一步。每一段用两个小区间构造二次插值多项式节点取两端点和中点。公式S_n h/3 * [f(x0) 4*Σf(x_odd) 2*Σf(x_even) f(xn)]系数模式就是1,4,2,4,2,...,4,1从 1 开始交替 4 和 2最后收在 1。这里有一个硬条件n 必须是偶数因为辛普森是按两段一组工作的n 为奇数时最后会剩一段找不到配对。整体误差-(b-a)/180 * h⁴ * f⁽⁴⁾(ξ)阶数是O(h⁴)。同样 h 减半误差降到原来的 1/16。同样是二分梯形四次翻倍才消掉一位十进制精度辛普森两次翻倍就能干完差距是实打实的。为什么不用更高次的 Newton-Cotes 公式高次插值在等距节点上会触发 Runge 现象多项式剧烈振荡数值结果完全不可信而且高次的 Newton-Cotes 系数会出现负数求和时引入灾难性抵消。所以工程实践里几乎见不到高次等距 Newton-Cotes大家要么用复化低阶要么用自适应策略。2.2 龙贝格算法理查森外推把梯形误差逐级磨平龙贝格算法的起点是一个很漂亮的数学事实复化梯形公式的误差可以写成 h 的偶次幂展开——T(h) I c1*h² c2*h⁴ c3*h⁶ ...。只要被积函数足够光滑这个展开就成立而且奇数幂全部消失。既然如此我用步长 h 和 h/2 各算一次梯形值两个式子相减就能把c1*h²这一项消掉T(h/2) (T(h/2) - T(h)) / 3 (4*T(h/2) - T(h)) / 3消掉 h² 项之后剩下的首项误差是 h⁴ 量级——这个结果恰好就是复化辛普森公式。这就是理查森外推用两个低精度的结果组合出一个高精度的结果。把同样的操作递归做下去用 h/4 的梯形值和刚才的组合值再消掉 h⁴ 项就得到 h⁶ 精度的结果。每一轮外推误差阶提高 2。实际实现时用一个三角表T(k,m)第 k 行代表步长(b-a)/2^kT(k,0)是该步长下的复化梯形值T(k,m)是对T(k,0)做 m 次外推的结果。外推公式T(k,m) (4^m * T(k,m-1) - T(k-1,m-1)) / (4^m - 1)这个公式里分母的4^m就是每次外推对应的误差阶修正因子。表的第 m 列误差阶为O(h^(2m2))对角线元素T(k,k)是当前能拿到的最高精度估算。龙贝格的工程价值在于它用梯形外推的组合拿到了高精度实现简单而且每一行只新增一半数量的函数求值性价比很高。缺点也同样明显——整套理论依赖光滑性被积函数有不连续点或者导数跳变展开式直接失效右下角的高精度结论全部作废。2.3 选型判断光滑度、精度目标和函数值成本怎么权衡我一般拿到一个积分问题先画函数图像再做选型判断而不是直接抄公式。判断逻辑就三条函数光滑不光滑、目标精度多高、一次函数求值贵不贵。下面这张表是我常用的速查方法误差阶每轮新增函数值适合场景复化梯形O(h²)n粗算、函数本身不光滑、一次求值极贵复化辛普森O(h⁴)n/2光滑函数中等精度实验课主力龙贝格随外推次数递增n/2光滑函数高精度报告加分项自适应辛普森局部 O(h⁴)局部加密区间有陡变但不知道在哪Gauss/Kronrod代数精度 2n-1n高精度且函数值便宜工程首选具体建议是实验课只要求在复化梯形、辛普森、龙贝格里选那就按精度需求来。要求1e-6以下辛普森从 n4 开始翻倍就能到要求1e-10级别直接用龙贝格加容差控制。如果函数在区间里有尖峰、振荡或者断点均匀网格的复化公式全都不合适——不是精度不够而是均匀采样根本打不到关键区域这时候该上自适应或者分段处理。函数求值成本这个因素学生容易忽略如果被积函数本身是一个仿真程序、一个数据库查询结果每算一次要几十毫秒那么梯形和龙贝格每轮只新增一半节点这个特性就非常值钱值得优先考虑。3. 用 Python 实现复化辛普森和龙贝格最小代码加交叉验证理论部分讲清楚后直接上代码。实验课最常见的作业就是手写复化辛普森和龙贝格下面两个实现是我多次用过的最小可靠版本不依赖复杂库只用到 NumPy。代码能跑通是最低要求更重要的是理解每一条语句在算什么、边界条件卡在哪里。3.1 复化辛普森最小实现向量化、奇偶性处理与参数含义import numpy as np def composite_simpson(f, a, b, n): 复化辛普森公式 f: 被积函数要求支持数组输入如 np.sin a, b: 积分区间 [a, b] n: 子区间数偶数传入奇数会自动加 1 返回: 积分近似值 # 辛普森要求两段一组n 必须为偶数 if n % 2 ! 0: n 1 if n 0: raise ValueError(n 必须是正整数) h (b - a) / n x np.linspace(a, b, n 1) y f(x) # 一次向量化求值比循环快一个数量级 # 权重系数: 1, 4, 2, 4, 2, ..., 4, 1 coef np.ones(n 1) coef[1:-1:2] 4 # 奇数下标内部节点权重 4 coef[2:-2:2] 2 # 偶数下标内部节点权重 2 return h / 3 * np.dot(coef, y)这段代码有几个必须理解的细节。第一f必须能直接吃进 numpy 数组。如果你在实验里写了math.sin然后传入x向量第一行f(x)就会报TypeError这是新手最常见的翻车点。要么函数内部用np.sin、np.exp这类向量化版本要么手动把x循环展开。第二权重系数用切片赋值实现coef[1:-1:2]取的是奇数下标内部节点coef[2:-2:2]取的是偶数下标内部节点端点保持 1。第三返回语句用np.dot做加权求和尽量避免 Python 原生循环——n 到几万时性能差距是百倍级的。参数怎么调n 是子区间数不是节点数。以[0,1]区间为例n4 表示切成 4 段、5 个节点权重模式是1,4,2,4,1。h 是每段长度(b-a)/n系数里的h/3来自辛普森公式对每一段的面积贡献。3.2 龙贝格三角表迭代、外推和收敛判据的代码翻译def romberg(f, a, b, tol1e-8, max_iter20): 龙贝格积分 返回: (积分近似值, 迭代层数, 误差估计) R np.zeros((max_iter, max_iter)) h b - a R[0, 0] h / 2 * (f(a) f(b)) # 第 0 层: 最粗的梯形值 for k in range(1, max_iter): # 第 k 层梯形值: 在上一层基础上补算新插入的中点 n_old 2 ** (k - 1) # 上一层的子区间数 h_old (b - a) / n_old # 上一层的步长 x_new a h_old / 2 h_old * np.arange(n_old) # 所有新区间中点 R[k, 0] 0.5 * R[k - 1, 0] h_old / 2 * np.sum(f(x_new)) # 理查森外推: 逐列消掉低阶误差项 for m in range(1, k 1): R[k, m] (4 ** m * R[k, m - 1] - R[k - 1, m - 1]) / (4 ** m - 1) # 收敛判断: 对角线相邻两层的差值 err abs(R[k, k] - R[k - 1, k - 1]) if err tol * max(1.0, abs(R[k, k])): # 绝对相对组合判据 return R[k, k], k, err return R[max_iter - 1, max_iter - 1], max_iter, \ abs(R[max_iter - 1, max_iter - 1] - R[max_iter - 2, max_iter - 2])核心逻辑在第 5 到第 9 行。R[k, 0]不是重新算整个梯形而是复用上一层的梯形值只补算新插入的区间中点。理解x_new的计算上一层有n_old个子区间每个新区间有一个中点这些中点的位置是a h_old/2 h_old * i其中i从 0 到n_old-1。这样每层只新增2^(k-1)次函数求值总代价比每次重算低得多。外推循环里的4 ** m是误差阶修正因子公式分母的4^m - 1保证外推组合的系数和为 1。收敛判据用的是对角线相邻两层的相对组合误差——单独用绝对误差在积分值本身很大时会过度苛刻单独用相对误差在积分值接近零时又会失效。这里有个工程细节R表开到max_iter * max_iter实际每层只用一行和一列。我在生产环境里会改成只保留上一行和当前行的两个一维数组省内存还能减少缓存 miss。实验课无所谓但养成这个意识不亏。3.3 用 scipy.integrate.quad 交叉验证手写实现到底准不准手写算法最大的风险是公式或者下标写错但结果看起来挺像那么回事。我的习惯是用成熟库做交叉验证一次性暴露问题。from scipy.integrate import quad # 以 sin(x) 在 [0,1] 的积分为例真值 1 - cos(1) f lambda x: np.sin(x) exact 1 - np.cos(1) val_simp composite_simpson(f, 0, 1, 64) val_romb, k, err_romb romberg(f, 0, 1, tol1e-10) val_quad, err_quad quad(f, 0, 1) # quad 返回 (积分值, 绝对误差估计) print(f辛普森: {val_simp:.12f}, 误差 {abs(val_simp-exact):.3e}) print(f龙贝格: {val_romb:.12f}, 迭代 {k} 层, 误差 {abs(val_romb-exact):.3e}) print(fquad : {val_quad:.12f}, 估计误差 {err_quad:.3e})用真值1 - cos(1)做标尺三种方法的结果应该在小数点后 8 位以内一致。如果对不上优先怀疑手写代码而不是 quad——quad 是 QUADPACK 的自适应 Gauss-Kronrod 实现对光滑函数非常可靠。需要注意两点一是quad返回的是一个(value, abserr)元组别只取第一个值然后拿整个元组去做算术二是 quad 也不是万能基准遇到高振荡函数它同样会给出错误结果所以交叉验证的测试函数要选sin、exp这类良态函数。4. 实验参数怎么定剖分数 n、龙贝格容差与误差阶验证实验报告能不能拿高分就看参数选得讲不讲道理。很多同学直接n1000一跑了事结果是算对了但答不出来为什么选 1000。这一章把三个最关键的参数决策讲清楚。4.1 剖分数 n 的选取翻倍策略和误差上界公式选 n 不是拍脑袋两个抓手误差上界公式和相邻结果对照。以复化辛普森为例误差上界是(b-a)/180 * h⁴ * M₄其中M₄是被积函数四阶导数绝对值的最大值。理论上你可以解这个不等式算出需要的 n但实际上算M₄常常比算积分还难所以我更推荐翻倍策略从 n4 开始每次都翻倍直到相邻两次结果差小于目标容差。以∫₀¹ sin(x) dx为例真值0.45969769413两种方法的误差随 n 变化n梯形误差辛普森误差29.6e-31.6e-442.4e-31.0e-586.0e-46.3e-7161.5e-43.9e-8看这个表就知道梯形每翻倍误差除以 4辛普森每翻倍误差除以 16。同样要 1e-8 精度辛普森算到 n32 左右梯形要 n 破千。所以辛普森是均匀网格下的默认选择梯形只适合粗算。翻倍策略还有一个好处它能顺便给你收敛速度的证据写报告时把相邻差值列出来老师一眼就能看出你理解了误差阶的概念。4.2 龙贝格容差设置绝对误差与相对误差的组合判据龙贝格代码里的tol怎么给直接决定程序是收敛后干净退出还是把 max_iter 跑穿后带着警告退出。我见过有人把tol1e-15然后抱怨程序怎么还不收敛——float64 的机器精度就在2.2e-16附近你要求一个1e-15的绝对误差在积分值本身是 1 到 10 的场景下等于是让算法去抠浮点表示的极限它抠不动。正确的做法是组合判据err tol * max(1.0, |integral|)。这个公式的含义是——积分值小于 1 时按绝对误差控制大于 1 时按相对误差控制max(1.0, ...)防止积分值恰好为零时判据失效。实验课我一般给tol1e-10既能压到 float64 的合理精度又不会触发极限。另一个必设参数是max_iter。龙贝格一层比一层贵如果被积函数不光滑导致怎么都不收敛max_iter就是保底的安全阀。我通常设 20对应2^20个节点对多数场景足够。4.3 误差阶验证log-log 拟合确认算法实现正确实验报告里最有说服力的一张图是误差随步长 h 变化的 log-log 曲线斜率就是误差阶。写一段代码把误差阶量化出来def verify_order(f, a, b, exact, scheme, ns(4, 8, 16, 32, 64, 128)): 验证某积分方案的误差阶 scheme: 函数句柄形如 scheme(f, a, b, n) 返回拟合出的误差阶斜率 hs, errs [], [] for n in ns: val scheme(f, a, b, n) err abs(val - exact) if err 1e-16: # 误差为 0 时 log 会崩给个下限 err 1e-16 hs.append((b - a) / n) errs.append(err) # log(err) log(C) p * log(h)斜率 p 就是误差阶 slope np.polyfit(np.log(hs), np.log(errs), 1)[0] return slope, hs, errs slope_simp, _, _ verify_order(f, 0, 1, exact, composite_simpson) slope_trap, _, _ verify_order(f, 0, 1, exact, composite_trapezoid) print(f辛普森误差阶 {slope_simp:.2f}预期 4) print(f梯形误差阶 {slope_trap:.2f}预期 2)np.polyfit在 log-log 空间做一次线性拟合返回的斜率就是误差阶。辛普森应该拟合出 4.0 左右梯形 2.0 左右。如果斜率明显低于理论值说明实现里有 bug或者 h 已经小到误差进入浮点饱和区——这解释了为什么前面要给误差一个下限。实践里我还会把饱和段的数据点在拟合前去掉否则斜率会被压弯。这张图和斜率值放实验报告里比贴十行运行输出都有用它直接证明你理解误差阶这个概念而不只是会调库。5. 数值积分避坑指南五个翻车现场和排查办法这一章是我自己踩坑攒出来的清单。每一条都是先给现象再给原因最后给解决路径对照排查比重新读理论快得多。5.1 高振荡函数辛普森的结果看起来收敛实际差得离谱现象对f(x) sin(100x)在[0, π]上做积分n64 时结果和 n128 只差零点几看起来收敛了和真值一对账差出一整个量级。原因辛普森假设每个子区间内的函数可以用二次多项式近似但 100Hz 的正弦在几个子区间内已经振荡了好几圈三次精度假设被暴力破坏而且均匀网格对高于网格分辨率的频率会直接失采样——采样点根本抓不到振荡的细节。解决这类问题不要用等距复化公式改用自适应积分或者先估一下函数的主频率保证每个振荡周期内至少 8 到 10 个采样点。用之前的romberg函数去算振荡积分也是同样的结局别把龙贝格当万能药。5.2 端点奇异直接算就是 NaN换元才能救现象算∫₀¹ 1/√x dxcomposite_simpson直接输出inf或者nan控制台还飘着 RuntimeWarning。原因被积函数在 x0 处趋于无穷零点是可积的瑕点但辛普森公式要求在被积函数在闭区间上有定义端点处的一个inf就把整个加权和污染了。解决换元把奇点消掉。令x t²则dx 2t dt积分变成∫₀¹ 2 dt 2完完全全的良性函数。这类问题更通用的处理是分段积分或者用含端点的自适应算法但换元永远是最省事的手段。注意scipy.integrate.quad对这类瑕积分一般能算出来但也会给出警告所以看到警告先别慌对照真值确认。5.3 大区间与远离原点浮点灾难性抵消让误差不降反升现象同一个函数在[0, 2]上积分精度 1e-14换到[1e8, 1e82]上积分误差直接掉到 1e-8 级别n 再加大也压不回去。原因float64 的相对精度固定约 1e-16x 到 1e8 时绝对精度只有 1e-8 量级。被积函数值在 1e16 附近积分结果约 2e16辛普森求和过程中 1e16 数量级的大数相减把低位的信息全部吞掉。这不是算法问题是浮点表示本身的边界。解决把积分区间平移或归一化让采样点回到原点附近再算如果必须在大数区间积分用math.fsum做求和它对灾难性抵消有一定抵抗能力或者把区间拆成小段分别积分。判断依据很简单看被积函数值和积分结果的数量级差差得越多越危险。5.4 龙贝格假收敛不连续函数不能硬套外推现象对被积函数为分段常数比如 x0.5 取 1否则取 2的积分跑romberg对角线误差一路降到 1e-7你以为收敛了但和真值 1.5 一比差 1e-2 量级。原因龙贝格的全部理论建立在梯形误差的 h 偶次幂展开上这个展开要求被积函数足够光滑。遇到间断点展开式失效理查森外推不仅消不掉低阶误差反而可能放大它表的右下角看起来很高精度实际全是废数据。解决先画函数图像检查连续性有间断就手动分段每段独立积分再求和或者改用自适应积分让算法自己去处理不连续点。写报告时如果实验给了这类函数主动分段比硬算然后拟合误差阶要体面得多因为你拟合出的斜率肯定会异常恰好证明你懂原理。5.5 lambda 闭包陷阱批量跑实验时结果全变同一个值现象循环里生成多个积分函数批量计算后所有结果一模一样像是函数没变。原因这是 Python 闭包的经典坑——lambda 内部引用的循环变量 n 是延迟绑定的循环结束后 n 绑定到了最后一个值所有 lambda 实际都是同一个函数。解决用默认参数把当前值钉死lambda x, nn: f(x, n)。或者干脆写一个普通的函数工厂返回内部函数。这个陷阱在实验课里频繁出现因为你经常要批量测多组参数、多个被积函数一旦中招错误是静默的——结果格式都正常就是数值不对排查起来很痛苦。我现在写批量脚本之前都会先跑一个最小用例验证闭包绑定这是省时间的习惯。6. 把实验四做厚误差阶可视化和自适应辛普森扩展实验做完只是起点。把结果可视化并且加一个自适应版本这份实验的价值会完全不同甚至能直接沉淀成一个可复用的小工具。6.1 log-log 误差阶图一页纸证明实现与算法一致把上一章的verify_order结果画出来代码很短import matplotlib.pyplot as plt slope_simp, hs_simp, errs_simp verify_order(f, 0, 1, exact, composite_simpson, ns(2, 4, 8, 16, 32, 64)) slope_trap, hs_trap, errs_trap verify_order(f, 0, 1, exact, composite_trapezoid, ns(2, 4, 8, 16, 32, 64)) plt.loglog(hs_trap, errs_trap, o-, labelf梯形, 斜率 {slope_trap:.2f}) plt.loglog(hs_simp, errs_simp, s-, labelf辛普森, 斜率 {slope_simp:.2f}) plt.loglog(hs_simp, errs_simp[0] * (np.array(hs_simp)/hs_simp[0])**4, --, labelh^4 参考线) plt.xlabel(步长 h) plt.ylabel(绝对误差) plt.legend() plt.grid(True, whichboth, ls--, alpha0.4) plt.show()读图的方法参考线和辛普森的实线如果平行说明误差确实按 h⁴ 衰减斜率数值也印证这一点。数据点如果在 h 特别小的地方开始变平那是浮点饱和区不是 bug。这张图放进实验报告配合一段三行的文字说明斜率含义比大段的运行截图有说服力得多。6.2 自适应辛普森跳过均匀网格的暴力采样自适应版本的思路是先用辛普森算整个区间再分成两半各算一次如果两半之和与整区间的差在容差范围内就接受并做一次误差补偿否则递归细分。真正的高效之处在于它只在函数变化剧烈的地方加密节点平缓区域保持粗网格。def adaptive_simpson(f, a, b, tol1e-8, max_depth50): def recurse(a, b, fa, fm, fb, segment_val, tol, depth): m (a b) / 2.0 lm, rm (a m) / 2.0, (m b) / 2.0 flm, frm f(lm), f(rm) left (m - a) / 6.0 * (fa 4.0 * flm fm) right (b - m) / 6.0 * (fm 4.0 * frm fb) delta left right - segment_val if depth max_depth or abs(delta) 15.0 * tol: return left right delta / 15.0 # 误差补偿项 return recurse(a, m, fa, flm, fm, left, tol/2, depth1) \ recurse(m, b, fm, frm, fb, right, tol/2, depth1) m (a b) / 2.0 return recurse(a, b, f(a), f(m), f(b), (b - a) / 6.0 * (f(a) 4.0 * f(m) f(b)), tol, 0)这里15*tol和delta/15来自辛普森误差估计的理论推导作用是让容差判断更准确并在返回时把已知误差补回去。max_depth是递归安全阀防止病态函数把栈压爆。用它对sin(100x)做积分结果明显优于均匀网格的复化辛普森因为递归自动扎堆到了振荡密集的地方。我现在的习惯是凡是工程里要用数值积分先花三十秒画函数草图和跨几个数量级的采样检查再决定用自适应还是分段还是龙贝格凡是给别人的积分代码一定带一个交叉验证的用例和明确的容差参数。这个流程帮我少踩了数不清的坑——很多问题不是算法不对而是参数没想清楚就往下跑了。希望帮到你。本文还有配套的精品资源点击获取