企业官网建设流程全解析

简介：面向数学、数值分析与计算科学方向的学生和研究者，这份本科毕业论文文档围绕几种常用数值积分方法展开系统比较。内容从数值积分基本思想入手，依次讨论复化求积公式、Newton-Cotes求积公式、Romberg求积公式与高斯型求积公式，并从代数精度、截断误差、绝对与相对误差等角度分析各自优缺点。作者还借助MATLAB上机实验，对不同被积函数下的求积效果进行对照，讨论精度与计算量的权衡，并给出方法选择上的建议。压缩包仅含1个doc文件，约897KB，类型单一但结构完整，包含开题报告、任务书、诚信声明及论文正文等模块。已有117人浏览学习，适合需要理解各类数值积分公式推导思路、误差特性和实验验证方式的读者作为参考资料，也便于课程论文或毕业设计阶段对照借鉴。

1. 从梯形公式说起：为什么工程上绕不开数值积分

大学里学定积分，第一反应是找原函数。真正进了工程计算才发现，能找到初等原函数的被积函数其实是小概率事件。e^{-x²}、sin(x)/x、√(1+x⁴) 这些看起来规整的函数，原函数一个都写不出初等形式；传感器采样得到的温度序列、CFD 里离散的压力场，更谈不上解析表达式。数值积分处理的正是这一类问题：不再追求原函数，而是挑一批节点、算节点处的函数值、做加权求和，把定积分逼近出来。本科阶段常见的四类方法——Newton-Cotes 求积、复化求积、Romberg 外推、Gauss-Legendre——覆盖了教学和工程实现的绝大部分场景，它们之间的差别集中在三件事上：代数精度、误差阶、节点怎么选。这篇把四种方法放在同一套评判标准下摆开，配合 MATLAB 复现脚本，把选型逻辑讲透。

2. 代数精度与余项：数值积分公式的评判标准

2.1 从曲边梯形分割到加权求和的一般形式

任何数值求积公式都能写成同一个形状：

I(f) = ∫_a^b f(x) dx ≈ Σ_{k=0}^{n} A_k f(x_k)

x_k 称为求积节点，A_k 称为求积系数（也叫伴随节点的权）。把定积分拆成小区间再求和，就是数值积分的基本思想——把整块曲边梯形的面积切成若干小曲边梯形的和，去掉取极限那一步，用有限个小曲边梯形的面积和代替整块面积。根据小区间的不同分割方法和分点处 f 值的不同选择，就得到了不同的数值积分公式。

写出一般形式后能立刻注意到一件事：权 A_k 只和节点的选取有关，与被积函数 f 的具体形式完全无关。这是插值型求积公式最核心的性质——一旦节点 x_k 定下来，A_k 就随之确定。这条性质让求积公式可以“一次构造、反复使用”，也是后面比较 Gauss 与 Newton-Cotes 的切入点。

2.2 代数精度：衡量公式“能精确到几阶”的硬指标

代数精度的定义很直白：如果求积公式对任何次数不超过 m 的多项式都精确成立，而对某个 m+1 次多项式不精确，则称该公式的代数精度为 m。

判断方式不需要推导，直接把 f(x) = 1, x, x², x³, … 依次代入，看余项什么时候第一次不为零：

• 梯形公式：代数精度 1（对线性函数精确，对 x² 不准）
• Simpson 公式：代数精度 3（比直觉的 2 高一级）
• Cotes 公式（n=4）：代数精度 5
• 一般 Newton-Cotes：n 为偶数时精度 n+1，n 为奇数时精度 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^(m+1) - a^(m+1)) / (m+1); % 精确值 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²)
• 复化 Simpson：O(h⁴)

这里的 h = (b-a)/n。h 减半，复化梯形误差缩到 1/4，复化 Simpson 缩到 1/16。收敛阶决定了加密步长的“性价比”，n 越大优势越明显。这也是为什么看似“只高了两阶”的 Simpson 在大规模计算里会成为默认选项。

2.4 节点选取自由度决定精度天花板

数一下自由度：n+1 个节点加 n+1 个权，一共 2n+2 个待定参数。Newton-Cotes 强行把节点固定为等距，只剩 n+1 个权能调，代数精度上限是 n+1（n 偶数时）。

如果把节点位置也放开当未知量一起解，理论上能精确到 2n+1 阶。这就是 Gauss 型求积公式的出发点。后面第 4 章会看到，同样用 4 个节点，Newton-Cotes 拿 5 阶精度，Gauss-Legendre 能拿 7 阶。

3. Newton-Cotes 与复化求积：等距节点的两条实现路径

3.1 等距节点上的插值求积：N-C 系数从哪来

构造思路是：在 [a,b] 上取 n+1 个等距节点 x_k = a + (b-a)k/n，用拉格朗日插值多项式 L_n(x) 替代 f(x)，然后对 L_n(x) 逐项积分，就得到求积公式。把 A_k 整理成 (b-a) 乘一个纯数字系数 C_k，这些数字就是 Cotes 系数，与区间无关，只与 n 有关：

MATLAB 实现一个通用的 N-C 求积：

function I = nc_quad(f, a, b, n) % n+1 个等距节点的 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('仅内置 n=1,2,3,4 的低阶 Cotes 系数'); end x = linspace(a, b, n+1); 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/(1+25x²) 在 [-1,1] 上试一下，精确值是 2arctan(5)/5：

f = @(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

实测在 n=8 之后误差不降反升。原因是高次插值多项式在端点附近剧烈振荡，积分时这些振荡被放大。所以 n ≥ 8 的 N-C 公式在工程里基本不用，需要提高精度就转向复化求积。

3.3 复化梯形与复化 Simpson：小步长收敛才是正解

复化梯形把 [a,b] 等分成 n 段，每段用梯形公式再求和：

T_n = h/2 · [ f(a) + 2Σ_{k=1}^{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⁴)。

对应 MATLAB：

function 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) 抓的是所有奇数下标的内部点，系数 4；y(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); end

n 每翻倍，梯形误差减为 1/4，Simpson 减为 1/16。Simpson 在 n=16 左右就已经逼近机器精度，梯形要 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 列的结果相当于一个高阶复化公式，误差阶按 2k+2 跳。

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 指数增长。M=10 单次要 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 把节点固定等距，浪费了 n+1 个自由度。Gauss 型公式让节点和权一起自由，求积公式

∫_{-1}^{1} f(t) dt ≈ Σ_{k=1}^{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))²] 给出。

常用节点和权：

一般区间 [a,b] 通过线性映射还原到 [-1,1]：x = (a+b)/2 + (b-a)/2 · t，dx = (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*(1+x)) / 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 在 x=0 处导数发散，被积函数不在 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 测体会不到方法之间的差别。建议至少覆盖三类：

跑一遍对照：

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 一个实用的选型速查表

实际落到代码时，如果项目里已经用了 SciPy，直接调 scipy.integrate.quad 是最省事的——它内部用自适应 Gauss-Kronrod，本质就是这里讲的 Gauss 节点加误差估计的组合。但要自己控制节点分布、处理离散采样表、或者把求积嵌到别的迭代里，前面这几个手写函数反而更好用：能直接拿 h 和节点位置当参数，不用跟自适应逻辑打交道。

本文还有配套的精品资源，点击获取