简介面向数学、数值分析与计算科学方向的学生和研究者这份本科毕业论文文档围绕几种常用数值积分方法展开系统比较。内容从数值积分基本思想入手依次讨论复化求积公式、Newton-Cotes求积公式、Romberg求积公式与高斯型求积公式并从代数精度、截断误差、绝对与相对误差等角度分析各自优缺点。作者还借助MATLAB上机实验对不同被积函数下的求积效果进行对照讨论精度与计算量的权衡并给出方法选择上的建议。压缩包仅含1个doc文件约897KB类型单一但结构完整包含开题报告、任务书、诚信声明及论文正文等模块。已有117人浏览学习适合需要理解各类数值积分公式推导思路、误差特性和实验验证方式的读者作为参考资料也便于课程论文或毕业设计阶段对照借鉴。1. 从梯形公式说起为什么工程上绕不开数值积分大学里学定积分第一反应是找原函数。真正进了工程计算才发现能找到初等原函数的被积函数其实是小概率事件。e^{-x²}、sin(x)/x、√(1x⁴) 这些看起来规整的函数原函数一个都写不出初等形式传感器采样得到的温度序列、CFD 里离散的压力场更谈不上解析表达式。数值积分处理的正是这一类问题不再追求原函数而是挑一批节点、算节点处的函数值、做加权求和把定积分逼近出来。本科阶段常见的四类方法——Newton-Cotes 求积、复化求积、Romberg 外推、Gauss-Legendre——覆盖了教学和工程实现的绝大部分场景它们之间的差别集中在三件事上代数精度、误差阶、节点怎么选。这篇把四种方法放在同一套评判标准下摆开配合 MATLAB 复现脚本把选型逻辑讲透。2. 代数精度与余项数值积分公式的评判标准2.1 从曲边梯形分割到加权求和的一般形式任何数值求积公式都能写成同一个形状I(f) ∫_a^b f(x) dx ≈ Σ_{k0}^{n} A_k f(x_k)x_k 称为求积节点A_k 称为求积系数也叫伴随节点的权。把定积分拆成小区间再求和就是数值积分的基本思想——把整块曲边梯形的面积切成若干小曲边梯形的和去掉取极限那一步用有限个小曲边梯形的面积和代替整块面积。根据小区间的不同分割方法和分点处 f 值的不同选择就得到了不同的数值积分公式。写出一般形式后能立刻注意到一件事权 A_k 只和节点的选取有关与被积函数 f 的具体形式完全无关。这是插值型求积公式最核心的性质——一旦节点 x_k 定下来A_k 就随之确定。这条性质让求积公式可以“一次构造、反复使用”也是后面比较 Gauss 与 Newton-Cotes 的切入点。2.2 代数精度衡量公式“能精确到几阶”的硬指标代数精度的定义很直白如果求积公式对任何次数不超过 m 的多项式都精确成立而对某个 m1 次多项式不精确则称该公式的代数精度为 m。判断方式不需要推导直接把 f(x) 1, x, x², x³, … 依次代入看余项什么时候第一次不为零梯形公式代数精度 1对线性函数精确对 x² 不准Simpson 公式代数精度 3比直觉的 2 高一级Cotes 公式n4代数精度 5一般 Newton-Cotesn 为偶数时精度 n1n 为奇数时精度 n对称性会白送一阶。区间 [a,b] 和节点都关于中点对称时奇次幂项的误差相互抵消所以 Simpson 对三次多项式也精确。用下面这段代码可以直接验证function deg alg_precision(quadfun, a, b) % 判断求积公式 quadfun 在 [a,b] 上的代数精度 % quadfun 接受 (f, a, b) 三个参数返回积分近似值 for m 0:10 f (x) x.^m; exact (b^(m1) - a^(m1)) / (m1); % 精确值 approx quadfun(f, a, b); if abs(approx - exact) 1e-12 * max(1, abs(exact)) deg m - 1; return; end end deg Inf; end公差 1e-12 用来防止浮点噪声误判。m 增大时 x^m 数值会很大实际跑的时候建议把 a、b 限制在 [0,1] 内。调用方式trap (f,a,b) (b-a)/2 * (f(a) f(b)); fprintf(梯形公式代数精度%d\n, alg_precision(trap, 0, 1));2.3 余项、误差阶与收敛速度代数精度评判的是“对多项式”误差阶评判的是“对光滑函数”。两者是互补的。梯形公式的余项R_T -(b-a)³/12 · f(ξ)ξ ∈ (a,b)。Simpson 公式的余项R_S -(b-a)⁵/2880 · f⁽⁴⁾(ξ)。当把区间等分并复化后误差随步长 h 缩复化梯形O(h²)复化 SimpsonO(h⁴)这里的 h (b-a)/n。h 减半复化梯形误差缩到 1/4复化 Simpson 缩到 1/16。收敛阶决定了加密步长的“性价比”n 越大优势越明显。这也是为什么看似“只高了两阶”的 Simpson 在大规模计算里会成为默认选项。2.4 节点选取自由度决定精度天花板数一下自由度n1 个节点加 n1 个权一共 2n2 个待定参数。Newton-Cotes 强行把节点固定为等距只剩 n1 个权能调代数精度上限是 n1n 偶数时。如果把节点位置也放开当未知量一起解理论上能精确到 2n1 阶。这就是 Gauss 型求积公式的出发点。后面第 4 章会看到同样用 4 个节点Newton-Cotes 拿 5 阶精度Gauss-Legendre 能拿 7 阶。3. Newton-Cotes 与复化求积等距节点的两条实现路径3.1 等距节点上的插值求积N-C 系数从哪来构造思路是在 [a,b] 上取 n1 个等距节点 x_k a (b-a)k/n用拉格朗日插值多项式 L_n(x) 替代 f(x)然后对 L_n(x) 逐项积分就得到求积公式。把 A_k 整理成 (b-a) 乘一个纯数字系数 C_k这些数字就是 Cotes 系数与区间无关只与 n 有关nCotes 系数代数精度常见名字11/2, 1/21梯形公式21/6, 4/6, 1/63Simpson 公式31/8, 3/8, 3/8, 1/833/8 公式47/90, 32/90, 12/90, 32/90, 7/905Cotes 公式MATLAB 实现一个通用的 N-C 求积function I nc_quad(f, a, b, n) % n1 个等距节点的 Newton-Cotes 求积 % n 建议取 1,2,4再大不要用 switch n case 1, C [1 1] / 2; case 2, C [1 4 1] / 6; case 3, C [1 3 3 1] / 8; case 4, C [7 32 12 32 7] / 90; otherwise error(仅内置 n1,2,3,4 的低阶 Cotes 系数); end x linspace(a, b, n1); I (b - a) * sum(C .* f(x)); end参数说明n 是节点数减 1直接对应插值多项式的次数。xz 用 linspace 生成等距节点(b-a) 乘 C 得到真实的权 A_k。系数之和恒为 1 可以用作一个快速自检。3.2 高次 N-C 为什么不能用高次插值的 Runge 现象出现在积分近似上一点不含糊。用 f(x) 1/(125x²) 在 [-1,1] 上试一下精确值是 2arctan(5)/5f (x) 1 ./ (1 25*x.^2); exact 2*atan(5)/5; for n [1 2 4 6 8 10 12] I nc_quad(f, -1, 1, n); fprintf(n%2d I%.6f err%.2e\n, n, I, abs(I-exact)); end实测在 n8 之后误差不降反升。原因是高次插值多项式在端点附近剧烈振荡积分时这些振荡被放大。所以 n ≥ 8 的 N-C 公式在工程里基本不用需要提高精度就转向复化求积。3.3 复化梯形与复化 Simpson小步长收敛才是正解复化梯形把 [a,b] 等分成 n 段每段用梯形公式再求和T_n h/2 · [ f(a) 2Σ_{k1}^{n-1} f(x_k) f(b) ]h (b-a)/n误差 O(h²)。复化 Simpson 每两个小区间合成一段用 Simpson 公式S_n h/3 · [ f(a) 4Σ_{奇数下标} f(x_k) 2Σ_{偶数下标非端点} f(x_k) f(b) ]误差 O(h⁴)。对应 MATLABfunction T comp_trap(f, a, b, n) % 复化梯形公式n 为分段数 h (b - a) / n; x a (0:n) * h; y f(x); % 向量化调用 T h * (y(1) y(end)) / 2 h * sum(y(2:end-1)); end function S comp_simpson(f, a, b, n) % 复化 Simpson 公式要求 n 为偶数 if mod(n, 2) ~ 0 error(n must be even for composite Simpson); end h (b - a) / n; x a (0:n) * h; y f(x); S h/3 * (y(1) y(end) ... 4*sum(y(2:2:end-1)) ... % 奇数下标系数 4 2*sum(y(3:2:end-2))); % 偶数内点系数 2 end参数说明n 是分段数复化梯形对奇偶无要求复化 Simpson 必须偶数。f 要支持向量输入内部用 .^ 和 .* 而不是 ^ 和 *。y(2:2:end-1) 抓的是所有奇数下标的内部点系数 4y(3:2:end-2) 抓偶数下标的内部点系数 2。3.4 用 e^x 实测收敛阶∫₀¹ e^x dx e - 1直接对比两种方法的误差随 n 的下降速度f (x) exp(x); exact exp(1) - 1; for n [4 8 16 32 64 128] eT abs(comp_trap(f, 0, 1, n) - exact); eS abs(comp_simpson(f, 0, 1, n) - exact); fprintf(n%3d 梯形 err%.3e Simpson err%.3e\n, n, eT, eS); endn 每翻倍梯形误差减为 1/4Simpson 减为 1/16。Simpson 在 n16 左右就已经逼近机器精度梯形要 n 上百次才追得上。这就是误差阶的直接体现O(h⁴) 相对于 O(h²) 的优势在大规模计算里非常可观。4. Romberg 外推与 Gauss-Legendre少节点高精度的两条路线4.1 Richardson 外推从复化梯形到 Romberg 表复化梯形的误差展开有个重要性质——只含 h 的偶次幂T(h) I c₁h² c₂h⁴ c₃h⁶ …既然只含偶次项那么拿两个不同步长的结果做线性组合就能消掉 h² 项T₁(h) ( 4T(h/2) - T(h) ) / 3这个 T₁(h) 正好就是复化 Simpson 公式。继续消T₂(h) ( 16T₁(h/2) - T₁(h) ) / 15得到 Cotes T₃(h) ( 64T₂(h/2) - T₂(h) ) / 63一般形式 T_m(h) ( 4^m · T_{m-1}(h/2) - T_{m-1}(h) ) / (4^m - 1)。按这个递推往上搭就得到 Romberg 三角表。第 k 列的结果相当于一个高阶复化公式误差阶按 2k2 跳。MATLAB 实现function R romberg(f, a, b, M) % Romberg 积分返回 M x M 三角表 % M 为外推层数工程上 5~8 足够 R zeros(M, M); h b - a; R(1,1) h/2 * (f(a) f(b)); % 复化梯形1 段 for k 2:M h h / 2; % 第 k 层复化梯形新增的节点之和 x_new a h * (1:2:(2^k - 1)); R(k,1) 0.5 * R(k-1,1) h * sum(f(x_new)); % Richardson 外推 for j 2:k R(k,j) (4^(j-1) * R(k,j-1) - R(k-1,j-1)) / (4^(j-1) - 1); end end end参数说明M外推层数。M 越大精度越高但每层新节点数按 2^k 指数增长。M10 单次要 1000 多次函数求值。R(k,1)第 k 层复化梯形结果为避免重复计算只加新增节点。R(k,j)第 j-1 阶外推结果j-1 阶对应误差阶 2j。对角线相邻两项之差 |R(k,k) - R(k-1,k-1)| 可以直接当误差估计用不需要额外计算。4.2 Gauss-Legendre把节点也变成自由度Newton-Cotes 把节点固定等距浪费了 n1 个自由度。Gauss 型公式让节点和权一起自由求积公式∫_{-1}^{1} f(t) dt ≈ Σ_{k1}^{n} A_k f(t_k)能够达到 2n-1 阶代数精度。节点 t_k 恰好是 n 次 Legendre 多项式 P_n(t) 的零点权由 A_k 2 / [(1-t_k²)(P_n(t_k))²] 给出。常用节点和权n节点 t_k权 A_k代数精度10212±1/√3 ≈ ±0.57735021, 1330, ±√(3/5) ≈ ±0.77459678/9, 5/9, 5/954±0.3399810, ±0.86113630.6521452, 0.34785487一般区间 [a,b] 通过线性映射还原到 [-1,1]x (ab)/2 (b-a)/2 · tdx (b-a)/2 dt。MATLAB 里不手抄节点表用对称三对角矩阵的特征值分解生成任意 n 的节点和权function [x, w] lgwt(n, a, b) % 生成 n 点 Gauss-Legendre 节点和权 % 基于 Jacobi 矩阵特征值分解n 到 20 也不掉精度 i (1:n-1); beta i ./ sqrt(4*i.^2 - 1); T diag(beta, 1) diag(beta, -1); % 对称三对角 [V, D] eig(T); x diag(D); % 特征值即节点 [x, idx] sort(x); w 2 * V(1, idx).^2; % 第一行分量平方乘 2 即权 x (a*(1-x) b*(1x)) / 2; % 映射到 [a,b] w w * (b - a) / 2; % 权乘以 Jacobian end function I gauss_legendre(f, a, b, n) [t, w] lgwt(n, -1, 1); x (a b)/2 (b - a)/2 * t; I (b - a)/2 * sum(w .* f(x)); end逻辑说明Legendre 多项式的零点就是那个 Jacobi 矩阵的特征值第一行分量的平方正比于对应的权。这种算法避免了直接牛顿迭代求零点时的初值问题n 到 20 依然稳定。参数说明n点数。Gauss 用 n 个点拿到 2n-1 阶精度。同样精度下复化梯形要 O(n²) 量级的点数才追上。映射线性变换会自动把权里的 Jacobian (b-a)/2 吸收进去。节点不包含端点这一点在奇异性场景下是有利的下一节讲。4.3 光滑性前提与常见误用Gauss-Legendre 的高精度有个硬前提被积函数在 [a,b] 上充分光滑。f 有端点奇异、跳变、导数不连续时2n-1 阶的理论精度会直接掉阶。例如 ∫₀¹ x^{0.5} dx用 4 点 Gauss-Legendre 出来的误差达不到理论量级。原因是 √x 在 x0 处导数发散被积函数不在 C⁴ 内。碰到这种场景有两个改法一是用 Gauss-Jacobi 型公式节点按 x^{α} 加权的正交多项式选α 取 0.5二是先做变量代换 x t²把奇异吸进去再套标准 Gauss-Legendre。另一个容易误用的点Gauss 节点不算端点如果被积函数在端点值恰好是 NaN 或 Inf但积分收敛反而不会被采样到。这和 Newton-Cotes 会踩到端点形成对比。真正需要端点信息的时候可以考虑 Gauss-Radau 或 Gauss-Lobatto那两种会把端点纳入节点。5. MATLAB 对比实验与选型判定5.1 三个测试用例别只测光滑函数只用 e^x 测体会不到方法之间的差别。建议至少覆盖三类编号被积函数区间精确值特性f₁exp(x)[0,1]e-1光滑f₂sin(20x)[0,2π]0振荡f₃sqrt(x)[0,1]2/3弱奇异跑一遍对照exact [exp(1)-1, 0, 2/3]; fs {(x) exp(x), (x) sin(20*x), (x) sqrt(x)}; ivals [0 1; 0 2*pi; 0 1]; names {光滑, 振荡, 弱奇异}; for k 1:3 f fs{k}; a ivals(k,1); b ivals(k,2); fprintf(--- %s ---\n, names{k}); for n [2 4 8 16] T comp_trap(f, a, b, n*4); S comp_simpson(f, a, b, n*4); G gauss_legendre(f, a, b, n); fprintf(n%2d 梯形%.2e Simpson%.2e Gauss%.2e\n, ... n, abs(T-exact(k)), abs(S-exact(k)), abs(G-exact(k))); end end弱奇异那一行会看到 Gauss 的误差下降到某个量级后卡住不动这正是 4.3 节讲的掉阶现象。5.2 loglog 图看收敛阶比看绝对误差靠谱绝对误差受比例因子影响收敛阶才是硬指标。横轴取 h纵轴取误差画双对数图N 2.^(2:8); h 2*pi ./ N; f (x) sin(20*x); errT arrayfun((n) abs(comp_trap(f, 0, 2*pi, n)), N); errS arrayfun((n) abs(comp_simpson(f, 0, 2*pi, n)), N); loglog(h, errT, o-, h, errS, s-, ... h, 1e-1*h.^2, k--, h, 1e0*h.^4, k:); legend(复化梯形,复化Simpson,参考 O(h^2),参考 O(h^4)); grid on;复化梯形的点会落在 O(h²) 参考线上复化 Simpson 落在 O(h⁴) 上。振荡函数 sin(20x) 有个前提每个振荡周期至少放 8~10 个点才会进入渐近区。n 太小的时候 Simpson 不一定比梯形好——这不是方法的问题是还没到渐近区就被误判了。5.3 一个实用的选型速查表场景首选方法关键理由被积函数光滑精度中等复化 Simpson实现简单O(h⁴) 够用函数值来自采样表等距复化梯形 / Simpson只能接受等距节点选不了 Gauss函数可任意求值要极高精度Gauss-Legendre少量节点逼近机器精度端点弱奇异Gauss-Jacobi 或先做变量代换Legendre 在高阶处丧失精度想要误差估计怕局部跳变Romberg 或自适应 Simpson三角表对角线之差可作误差上界实际落到代码时如果项目里已经用了 SciPy直接调 scipy.integrate.quad 是最省事的——它内部用自适应 Gauss-Kronrod本质就是这里讲的 Gauss 节点加误差估计的组合。但要自己控制节点分布、处理离散采样表、或者把求积嵌到别的迭代里前面这几个手写函数反而更好用能直接拿 h 和节点位置当参数不用跟自适应逻辑打交道。本文还有配套的精品资源点击获取
数值积分方法对比：复化求积、Romberg与Gauss-Legendre选型
简介面向数学、数值分析与计算科学方向的学生和研究者这份本科毕业论文文档围绕几种常用数值积分方法展开系统比较。内容从数值积分基本思想入手依次讨论复化求积公式、Newton-Cotes求积公式、Romberg求积公式与高斯型求积公式并从代数精度、截断误差、绝对与相对误差等角度分析各自优缺点。作者还借助MATLAB上机实验对不同被积函数下的求积效果进行对照讨论精度与计算量的权衡并给出方法选择上的建议。压缩包仅含1个doc文件约897KB类型单一但结构完整包含开题报告、任务书、诚信声明及论文正文等模块。已有117人浏览学习适合需要理解各类数值积分公式推导思路、误差特性和实验验证方式的读者作为参考资料也便于课程论文或毕业设计阶段对照借鉴。1. 从梯形公式说起为什么工程上绕不开数值积分大学里学定积分第一反应是找原函数。真正进了工程计算才发现能找到初等原函数的被积函数其实是小概率事件。e^{-x²}、sin(x)/x、√(1x⁴) 这些看起来规整的函数原函数一个都写不出初等形式传感器采样得到的温度序列、CFD 里离散的压力场更谈不上解析表达式。数值积分处理的正是这一类问题不再追求原函数而是挑一批节点、算节点处的函数值、做加权求和把定积分逼近出来。本科阶段常见的四类方法——Newton-Cotes 求积、复化求积、Romberg 外推、Gauss-Legendre——覆盖了教学和工程实现的绝大部分场景它们之间的差别集中在三件事上代数精度、误差阶、节点怎么选。这篇把四种方法放在同一套评判标准下摆开配合 MATLAB 复现脚本把选型逻辑讲透。2. 代数精度与余项数值积分公式的评判标准2.1 从曲边梯形分割到加权求和的一般形式任何数值求积公式都能写成同一个形状I(f) ∫_a^b f(x) dx ≈ Σ_{k0}^{n} A_k f(x_k)x_k 称为求积节点A_k 称为求积系数也叫伴随节点的权。把定积分拆成小区间再求和就是数值积分的基本思想——把整块曲边梯形的面积切成若干小曲边梯形的和去掉取极限那一步用有限个小曲边梯形的面积和代替整块面积。根据小区间的不同分割方法和分点处 f 值的不同选择就得到了不同的数值积分公式。写出一般形式后能立刻注意到一件事权 A_k 只和节点的选取有关与被积函数 f 的具体形式完全无关。这是插值型求积公式最核心的性质——一旦节点 x_k 定下来A_k 就随之确定。这条性质让求积公式可以“一次构造、反复使用”也是后面比较 Gauss 与 Newton-Cotes 的切入点。2.2 代数精度衡量公式“能精确到几阶”的硬指标代数精度的定义很直白如果求积公式对任何次数不超过 m 的多项式都精确成立而对某个 m1 次多项式不精确则称该公式的代数精度为 m。判断方式不需要推导直接把 f(x) 1, x, x², x³, … 依次代入看余项什么时候第一次不为零梯形公式代数精度 1对线性函数精确对 x² 不准Simpson 公式代数精度 3比直觉的 2 高一级Cotes 公式n4代数精度 5一般 Newton-Cotesn 为偶数时精度 n1n 为奇数时精度 n对称性会白送一阶。区间 [a,b] 和节点都关于中点对称时奇次幂项的误差相互抵消所以 Simpson 对三次多项式也精确。用下面这段代码可以直接验证function deg alg_precision(quadfun, a, b) % 判断求积公式 quadfun 在 [a,b] 上的代数精度 % quadfun 接受 (f, a, b) 三个参数返回积分近似值 for m 0:10 f (x) x.^m; exact (b^(m1) - a^(m1)) / (m1); % 精确值 approx quadfun(f, a, b); if abs(approx - exact) 1e-12 * max(1, abs(exact)) deg m - 1; return; end end deg Inf; end公差 1e-12 用来防止浮点噪声误判。m 增大时 x^m 数值会很大实际跑的时候建议把 a、b 限制在 [0,1] 内。调用方式trap (f,a,b) (b-a)/2 * (f(a) f(b)); fprintf(梯形公式代数精度%d\n, alg_precision(trap, 0, 1));2.3 余项、误差阶与收敛速度代数精度评判的是“对多项式”误差阶评判的是“对光滑函数”。两者是互补的。梯形公式的余项R_T -(b-a)³/12 · f(ξ)ξ ∈ (a,b)。Simpson 公式的余项R_S -(b-a)⁵/2880 · f⁽⁴⁾(ξ)。当把区间等分并复化后误差随步长 h 缩复化梯形O(h²)复化 SimpsonO(h⁴)这里的 h (b-a)/n。h 减半复化梯形误差缩到 1/4复化 Simpson 缩到 1/16。收敛阶决定了加密步长的“性价比”n 越大优势越明显。这也是为什么看似“只高了两阶”的 Simpson 在大规模计算里会成为默认选项。2.4 节点选取自由度决定精度天花板数一下自由度n1 个节点加 n1 个权一共 2n2 个待定参数。Newton-Cotes 强行把节点固定为等距只剩 n1 个权能调代数精度上限是 n1n 偶数时。如果把节点位置也放开当未知量一起解理论上能精确到 2n1 阶。这就是 Gauss 型求积公式的出发点。后面第 4 章会看到同样用 4 个节点Newton-Cotes 拿 5 阶精度Gauss-Legendre 能拿 7 阶。3. Newton-Cotes 与复化求积等距节点的两条实现路径3.1 等距节点上的插值求积N-C 系数从哪来构造思路是在 [a,b] 上取 n1 个等距节点 x_k a (b-a)k/n用拉格朗日插值多项式 L_n(x) 替代 f(x)然后对 L_n(x) 逐项积分就得到求积公式。把 A_k 整理成 (b-a) 乘一个纯数字系数 C_k这些数字就是 Cotes 系数与区间无关只与 n 有关nCotes 系数代数精度常见名字11/2, 1/21梯形公式21/6, 4/6, 1/63Simpson 公式31/8, 3/8, 3/8, 1/833/8 公式47/90, 32/90, 12/90, 32/90, 7/905Cotes 公式MATLAB 实现一个通用的 N-C 求积function I nc_quad(f, a, b, n) % n1 个等距节点的 Newton-Cotes 求积 % n 建议取 1,2,4再大不要用 switch n case 1, C [1 1] / 2; case 2, C [1 4 1] / 6; case 3, C [1 3 3 1] / 8; case 4, C [7 32 12 32 7] / 90; otherwise error(仅内置 n1,2,3,4 的低阶 Cotes 系数); end x linspace(a, b, n1); I (b - a) * sum(C .* f(x)); end参数说明n 是节点数减 1直接对应插值多项式的次数。xz 用 linspace 生成等距节点(b-a) 乘 C 得到真实的权 A_k。系数之和恒为 1 可以用作一个快速自检。3.2 高次 N-C 为什么不能用高次插值的 Runge 现象出现在积分近似上一点不含糊。用 f(x) 1/(125x²) 在 [-1,1] 上试一下精确值是 2arctan(5)/5f (x) 1 ./ (1 25*x.^2); exact 2*atan(5)/5; for n [1 2 4 6 8 10 12] I nc_quad(f, -1, 1, n); fprintf(n%2d I%.6f err%.2e\n, n, I, abs(I-exact)); end实测在 n8 之后误差不降反升。原因是高次插值多项式在端点附近剧烈振荡积分时这些振荡被放大。所以 n ≥ 8 的 N-C 公式在工程里基本不用需要提高精度就转向复化求积。3.3 复化梯形与复化 Simpson小步长收敛才是正解复化梯形把 [a,b] 等分成 n 段每段用梯形公式再求和T_n h/2 · [ f(a) 2Σ_{k1}^{n-1} f(x_k) f(b) ]h (b-a)/n误差 O(h²)。复化 Simpson 每两个小区间合成一段用 Simpson 公式S_n h/3 · [ f(a) 4Σ_{奇数下标} f(x_k) 2Σ_{偶数下标非端点} f(x_k) f(b) ]误差 O(h⁴)。对应 MATLABfunction T comp_trap(f, a, b, n) % 复化梯形公式n 为分段数 h (b - a) / n; x a (0:n) * h; y f(x); % 向量化调用 T h * (y(1) y(end)) / 2 h * sum(y(2:end-1)); end function S comp_simpson(f, a, b, n) % 复化 Simpson 公式要求 n 为偶数 if mod(n, 2) ~ 0 error(n must be even for composite Simpson); end h (b - a) / n; x a (0:n) * h; y f(x); S h/3 * (y(1) y(end) ... 4*sum(y(2:2:end-1)) ... % 奇数下标系数 4 2*sum(y(3:2:end-2))); % 偶数内点系数 2 end参数说明n 是分段数复化梯形对奇偶无要求复化 Simpson 必须偶数。f 要支持向量输入内部用 .^ 和 .* 而不是 ^ 和 *。y(2:2:end-1) 抓的是所有奇数下标的内部点系数 4y(3:2:end-2) 抓偶数下标的内部点系数 2。3.4 用 e^x 实测收敛阶∫₀¹ e^x dx e - 1直接对比两种方法的误差随 n 的下降速度f (x) exp(x); exact exp(1) - 1; for n [4 8 16 32 64 128] eT abs(comp_trap(f, 0, 1, n) - exact); eS abs(comp_simpson(f, 0, 1, n) - exact); fprintf(n%3d 梯形 err%.3e Simpson err%.3e\n, n, eT, eS); endn 每翻倍梯形误差减为 1/4Simpson 减为 1/16。Simpson 在 n16 左右就已经逼近机器精度梯形要 n 上百次才追得上。这就是误差阶的直接体现O(h⁴) 相对于 O(h²) 的优势在大规模计算里非常可观。4. Romberg 外推与 Gauss-Legendre少节点高精度的两条路线4.1 Richardson 外推从复化梯形到 Romberg 表复化梯形的误差展开有个重要性质——只含 h 的偶次幂T(h) I c₁h² c₂h⁴ c₃h⁶ …既然只含偶次项那么拿两个不同步长的结果做线性组合就能消掉 h² 项T₁(h) ( 4T(h/2) - T(h) ) / 3这个 T₁(h) 正好就是复化 Simpson 公式。继续消T₂(h) ( 16T₁(h/2) - T₁(h) ) / 15得到 Cotes T₃(h) ( 64T₂(h/2) - T₂(h) ) / 63一般形式 T_m(h) ( 4^m · T_{m-1}(h/2) - T_{m-1}(h) ) / (4^m - 1)。按这个递推往上搭就得到 Romberg 三角表。第 k 列的结果相当于一个高阶复化公式误差阶按 2k2 跳。MATLAB 实现function R romberg(f, a, b, M) % Romberg 积分返回 M x M 三角表 % M 为外推层数工程上 5~8 足够 R zeros(M, M); h b - a; R(1,1) h/2 * (f(a) f(b)); % 复化梯形1 段 for k 2:M h h / 2; % 第 k 层复化梯形新增的节点之和 x_new a h * (1:2:(2^k - 1)); R(k,1) 0.5 * R(k-1,1) h * sum(f(x_new)); % Richardson 外推 for j 2:k R(k,j) (4^(j-1) * R(k,j-1) - R(k-1,j-1)) / (4^(j-1) - 1); end end end参数说明M外推层数。M 越大精度越高但每层新节点数按 2^k 指数增长。M10 单次要 1000 多次函数求值。R(k,1)第 k 层复化梯形结果为避免重复计算只加新增节点。R(k,j)第 j-1 阶外推结果j-1 阶对应误差阶 2j。对角线相邻两项之差 |R(k,k) - R(k-1,k-1)| 可以直接当误差估计用不需要额外计算。4.2 Gauss-Legendre把节点也变成自由度Newton-Cotes 把节点固定等距浪费了 n1 个自由度。Gauss 型公式让节点和权一起自由求积公式∫_{-1}^{1} f(t) dt ≈ Σ_{k1}^{n} A_k f(t_k)能够达到 2n-1 阶代数精度。节点 t_k 恰好是 n 次 Legendre 多项式 P_n(t) 的零点权由 A_k 2 / [(1-t_k²)(P_n(t_k))²] 给出。常用节点和权n节点 t_k权 A_k代数精度10212±1/√3 ≈ ±0.57735021, 1330, ±√(3/5) ≈ ±0.77459678/9, 5/9, 5/954±0.3399810, ±0.86113630.6521452, 0.34785487一般区间 [a,b] 通过线性映射还原到 [-1,1]x (ab)/2 (b-a)/2 · tdx (b-a)/2 dt。MATLAB 里不手抄节点表用对称三对角矩阵的特征值分解生成任意 n 的节点和权function [x, w] lgwt(n, a, b) % 生成 n 点 Gauss-Legendre 节点和权 % 基于 Jacobi 矩阵特征值分解n 到 20 也不掉精度 i (1:n-1); beta i ./ sqrt(4*i.^2 - 1); T diag(beta, 1) diag(beta, -1); % 对称三对角 [V, D] eig(T); x diag(D); % 特征值即节点 [x, idx] sort(x); w 2 * V(1, idx).^2; % 第一行分量平方乘 2 即权 x (a*(1-x) b*(1x)) / 2; % 映射到 [a,b] w w * (b - a) / 2; % 权乘以 Jacobian end function I gauss_legendre(f, a, b, n) [t, w] lgwt(n, -1, 1); x (a b)/2 (b - a)/2 * t; I (b - a)/2 * sum(w .* f(x)); end逻辑说明Legendre 多项式的零点就是那个 Jacobi 矩阵的特征值第一行分量的平方正比于对应的权。这种算法避免了直接牛顿迭代求零点时的初值问题n 到 20 依然稳定。参数说明n点数。Gauss 用 n 个点拿到 2n-1 阶精度。同样精度下复化梯形要 O(n²) 量级的点数才追上。映射线性变换会自动把权里的 Jacobian (b-a)/2 吸收进去。节点不包含端点这一点在奇异性场景下是有利的下一节讲。4.3 光滑性前提与常见误用Gauss-Legendre 的高精度有个硬前提被积函数在 [a,b] 上充分光滑。f 有端点奇异、跳变、导数不连续时2n-1 阶的理论精度会直接掉阶。例如 ∫₀¹ x^{0.5} dx用 4 点 Gauss-Legendre 出来的误差达不到理论量级。原因是 √x 在 x0 处导数发散被积函数不在 C⁴ 内。碰到这种场景有两个改法一是用 Gauss-Jacobi 型公式节点按 x^{α} 加权的正交多项式选α 取 0.5二是先做变量代换 x t²把奇异吸进去再套标准 Gauss-Legendre。另一个容易误用的点Gauss 节点不算端点如果被积函数在端点值恰好是 NaN 或 Inf但积分收敛反而不会被采样到。这和 Newton-Cotes 会踩到端点形成对比。真正需要端点信息的时候可以考虑 Gauss-Radau 或 Gauss-Lobatto那两种会把端点纳入节点。5. MATLAB 对比实验与选型判定5.1 三个测试用例别只测光滑函数只用 e^x 测体会不到方法之间的差别。建议至少覆盖三类编号被积函数区间精确值特性f₁exp(x)[0,1]e-1光滑f₂sin(20x)[0,2π]0振荡f₃sqrt(x)[0,1]2/3弱奇异跑一遍对照exact [exp(1)-1, 0, 2/3]; fs {(x) exp(x), (x) sin(20*x), (x) sqrt(x)}; ivals [0 1; 0 2*pi; 0 1]; names {光滑, 振荡, 弱奇异}; for k 1:3 f fs{k}; a ivals(k,1); b ivals(k,2); fprintf(--- %s ---\n, names{k}); for n [2 4 8 16] T comp_trap(f, a, b, n*4); S comp_simpson(f, a, b, n*4); G gauss_legendre(f, a, b, n); fprintf(n%2d 梯形%.2e Simpson%.2e Gauss%.2e\n, ... n, abs(T-exact(k)), abs(S-exact(k)), abs(G-exact(k))); end end弱奇异那一行会看到 Gauss 的误差下降到某个量级后卡住不动这正是 4.3 节讲的掉阶现象。5.2 loglog 图看收敛阶比看绝对误差靠谱绝对误差受比例因子影响收敛阶才是硬指标。横轴取 h纵轴取误差画双对数图N 2.^(2:8); h 2*pi ./ N; f (x) sin(20*x); errT arrayfun((n) abs(comp_trap(f, 0, 2*pi, n)), N); errS arrayfun((n) abs(comp_simpson(f, 0, 2*pi, n)), N); loglog(h, errT, o-, h, errS, s-, ... h, 1e-1*h.^2, k--, h, 1e0*h.^4, k:); legend(复化梯形,复化Simpson,参考 O(h^2),参考 O(h^4)); grid on;复化梯形的点会落在 O(h²) 参考线上复化 Simpson 落在 O(h⁴) 上。振荡函数 sin(20x) 有个前提每个振荡周期至少放 8~10 个点才会进入渐近区。n 太小的时候 Simpson 不一定比梯形好——这不是方法的问题是还没到渐近区就被误判了。5.3 一个实用的选型速查表场景首选方法关键理由被积函数光滑精度中等复化 Simpson实现简单O(h⁴) 够用函数值来自采样表等距复化梯形 / Simpson只能接受等距节点选不了 Gauss函数可任意求值要极高精度Gauss-Legendre少量节点逼近机器精度端点弱奇异Gauss-Jacobi 或先做变量代换Legendre 在高阶处丧失精度想要误差估计怕局部跳变Romberg 或自适应 Simpson三角表对角线之差可作误差上界实际落到代码时如果项目里已经用了 SciPy直接调 scipy.integrate.quad 是最省事的——它内部用自适应 Gauss-Kronrod本质就是这里讲的 Gauss 节点加误差估计的组合。但要自己控制节点分布、处理离散采样表、或者把求积嵌到别的迭代里前面这几个手写函数反而更好用能直接拿 h 和节点位置当参数不用跟自适应逻辑打交道。本文还有配套的精品资源点击获取