企业官网建设流程全解析

简介：「数值分析上机题」配套上机实践资料，面向修读数值分析课程、需要用 MATLAB 完成课后编程作业的高校学生与自学者，内容围绕三次样条插值函数、重积分计算与常微分方程初值问题数值解三个上机专题展开。其中，三次样条插值部分给出第一型三次样条插值函数的通用程序编制思路与车门曲线型值点算例；重积分部分采用复化梯形公式，配合逐次二分步长与外推思想实现通用计算程序；常微分方程部分则对 RK4、AB4、AB4-AM4 及改进预测校正四种方法进行编程实现，并以二范数误差比较各方法精度。压缩包共 1 个 docx 文件，约 566KB，是一份完整的上机作业文档。已有 905 人学习，可对照公式推导、程序组织与误差结论查漏补缺，也适合作为 MATLAB 数值算法入门的参考范例。

1. 三道上机题，最后拼的是排错能力

数值分析上机题这套东西，代码量加起来不到三百行，但真正让人熬夜的不是公式推不出来，是推出来的公式对不上程序跑出来的数。我印象最深的一次翻车：三次样条的线性方程组解出来了，代回去一打印，第一段曲线在端点附近明显翘起来一块，看着像插值点选错了，其实是首行边界条件的符号写反，一个正负号折腾了半宿。还有一次 RK4 单步看着很准，把四种方法摆在一张表里比误差，排序却和课本结论对不上，最后发现是 AB4 的历史项从第二个点就开始取，取到了未初始化的零。用 MATLAB 写这些通用程序，语言本身不难，难的是每个下标对应公式里的哪一个、边界条件怎么进方程、二分的终止判据挂在哪一层。下面是三道题的完整骨架，参数含义、失败时看哪里都写在代码旁边，新手能照着跑通，写过一遍的人能对上自己当时没想明白的那一处。

2. 三次样条插值：从 11 个型值点到 S(x) 的完整链路

2.1 为什么拿节点二阶导当未知量

三次样条的做法是分段三次、整体二阶连续，直接拿每段的四个系数当未知量，会得到 4n 个方程，写起来整齐、解起来笨重。课本走的是另一条路：把每个节点上的二阶导数值 M_j 当作未知量，因为三次多项式在区间上是二阶导线性函数，只要两端的 M 定下来，整段的表达式就唯一确定，剩下的约束自动满足。这一换元把问题压成 n+1 个未知量的三对角方程组，也就是常见的弯矩法。

给定 n+1 个型值点，最后一条公式长这样：

μ_j M_{j-1} + 2 M_j + α_j M_{j+1} = d_j， j = 1, 2, …, n-1

其中 μ_j = h_{j-1}/(h_{j-1}+h_j)，α_j = h_j/(h_{j-1}+h_j) = 1 - μ_j，d_j = 6·f[x_{j-1}, x_j, x_{j+1}]。注意 μ 用的是左边那一段的步长，α 用的是右边那一段，两者交换位置是这一步最容易写错的地方。手里这道车门轮廓的题，x 从 0 到 10 均匀分布，h 全是 1，于是 μ_j = α_j = 0.5，方程退化成 (M_{j-1} + 4M_j + M_{j+1})/2 = d_j，形式上简单很多，但也正因为简单，很多人会直接把 0.5 硬编码进去，换成非均匀节点就整个崩掉。

提示：先把 h、一阶差商、二阶差商三样算成独立向量，再组装矩阵。不要在循环里边算边填，出错了根本分不清是差商错了还是系数错了。

2.2 第一型边界条件怎么塞进首尾两行

中间 n-1 个方程只够定 n-1 个未知量，剩下的两个靠端点条件补。这道题给的是第一型，即两端的一阶导数值 y_0' = 0.8、y_10' = 0.2，也就是固定斜率。把 S'(x_0) 用 M_0、M_1 表示出来，整理后首行是：

2M_0 + M_1 = (6/h_0)·(f[x_0, x_1] - y_0')

尾行对称：

M_{n-1} + 2M_n = (6/h_{n-1})·(y_n' - f[x_{n-1}, x_n])

这两行的符号方向不一样，首行是「差商减导数」，尾行是「导数减差商」，原因是一阶导在左端点由区间内部向左取、在右端点由左向右取，方向反了。我见过的符号性翻车，八成出在这里：把尾行也写成差商减导数，结果整条曲线在右端点附近往下压，插值点还是准的，但两端斜率全错，只看中间几个数根本发现不了。

如果换成第二型自然边界，首尾两行变成 M_0 = 0 和 M_n = 0，代码只改这两行的赋值，其余不动；第三型周期边界则要把首尾两行都换成周期关系，并把第一个和最后一个未知量绑在一起。写通用程序的时候，我一般会留一个 switch 参数控制边界类型，而不是三份代码各写一遍，因为三份代码里必有一份是过时的。

2.3 组装与求解：一份能直接跑的骨架

下面这段是我当时那版的核心，参数全是显式传入，方便换数据验证。

function [M, S] = spline3_clamped(x, y, dy0, dyn) % x, y 为型值点（等长行向量或列向量） % dy0, dyn 为两端一阶导数值（第一型边界条件） x = x(:); y = y(:); n = length(x) - 1; % 区间个数 h = diff(x); % h(j) = x(j+1) - x(j) f1 = diff(y) ./ h; % 一阶差商 f1(j) = f[x(j), x(j+1)] % 二阶差商 f2(k) = f[x(k), x(k+1), x(k+2)]，共 n-1 个 f2 = (f1(2:end) - f1(1:end-1)) ./ (x(3:end) - x(1:end-2)); A = zeros(n+1); rhs = zeros(n+1, 1); % 首行：2*M0 + M1 = 6/h0 * (f[x0,x1] - y0') A(1,1) = 2; A(1,2) = 1; rhs(1) = 6 / h(1) * (f1(1) - dy0); for j = 2:n % 对应节点 x(j)，j 从 2 到 n mu = h(j-1) / (h(j-1) + h(j)); alp = 1 - mu; A(j, j-1) = mu; A(j, j) = 2; A(j, j+1) = alp; rhs(j) = 6 * f2(j-1); % 注意下标错位一位 end % 尾行：M(n-1) + 2*Mn = 6/hn * (yn' - f[x(n-1),x(n)]) A(n+1, n) = 1; A(n+1, n+1) = 2; rhs(n+1) = 6 / h(n) * (dyn - f1(n)); M = A \ rhs; % 三对角也可手写追赶法，结果一致 end

几个地方值得单独说。f2(j-1)这一处下标错位是整段代码里最隐蔽的：j 是节点编号，而 f2 的第 k 个元素对应的是以 x(k) 开头的三点差商，两者天然差一位，写成f2(j)在 n=10 时会直接越界报错，运气好报错，运气不好节点数刚好够，结果整体偏移一个区间，插值点全部对不上。A\rhs用的是全矩阵左除，教学上要求写追赶法的话，把这三行换成前代回代两段循环即可，三对角矩阵不必担心数值稳定性，主对角占优是成立的。

求出 M 之后，区间 [x_j, x_{j+1}] 上的表达式按弯矩形式写：

S(x) = a·y_j + b·y_{j+1} + [(a³ - a)·M_j + (b³ - b)·M_{j+1}]·h_j²/6

其中 a = (x_{j+1} - x)/h_j，b = (x - x_j)/h_j。a + b 恒等于 1，这一点可以直接拿来当自检：随便取一个 x，算出 a、b，加起来不为 1 就说明区间定位错了。

function v = evalS(x, y, M, t) % 在 t（可为向量）处求三次样条值 x = x(:); y = y(:); M = M(:); v = zeros(size(t)); for k = 1:numel(t) j = find(x <= t(k), 1, 'last'); j = min(j, numel(x) - 1); % 右端点保护，避免 j 超出区间数 hj = x(j+1) - x(j); a = (x(j+1) - t(k)) / hj; b = (t(k) - x(j)) / hj; v(k) = a*y(j) + b*y(j+1) + ((a^3-a)*M(j) + (b^3-b)*M(j+1))*hj^2/6; end end

题目要的是 S(i+0.5)，i 从 0 到 9，直接evalS(x, y, M, 0.5:1:9.5)一行出结果。

2.4 跑完先做三件事再交

第一件，把原始型值点代回去，S(x_i) 必须精确等于 y_i，差超过 1e-10 就是矩阵或下标有问题，不是精度问题。第二件，用数值差分估一下两端斜率，(S(0.001) - S(0))/0.001 应该贴着 0.8，(S(10) - S(9.999))/0.001 应该贴着 0.2，对不上就是边界条件那两行符号反了。第三件，把 x 加密画图，看二阶导 M 的曲线，正常情况下 M 应该是连续变化的，如果出现锯齿状跳变，多半是某一段的 h 用错了邻段。

这三件事加起来不超过十行代码，但能挡掉九成的返工。我当时就是跳过了第二件，直接对着打印出来的 S(i+0.5) 一列数字看，看了半小时也看不出哪里不对，因为中间那些数看着都「合理」。

3. 重积分与外推：把复化梯形到 Simpson 的链路拆成三个函数

3.1 动手之前先看被积函数在这个区域里有没有奇点

这道题要算的是外层对 y、内层对 x 的二重积分，被积函数 tan(x² + y²)，x 在 [0, π/3]，y 在 [0, π/6]，精度要求 ε = 12×10⁻⁵。tan 在 π/2 附近会炸，所以在写任何循环之前，先算一下区域内 x² + y² 的最大值：x 取 π/3 时 x² ≈ 1.0966，y 取 π/6 时 y² ≈ 0.2742，加起来 1.3708，离 π/2 ≈ 1.5708 还有一段距离，被积函数在整个矩形区域上光滑有界，可以放心二分。

这一步看着像废话，但它决定了后面所有判断的合法性。如果区域边界刚好贴到 π/2，复化梯形在奇点附近根本不收敛，二分再多次也是白搭，得先做变量替换或者分段处理。很多人拿到题直接敲循环，跑了十几层发现两次结果差得越来越远，才回头怀疑公式，其实问题在更前面。

注意：判断收敛性之前，先判断被积函数在区域上有没有奇点、有没有剧烈振荡。这两条不过关，任何步长策略都是空谈。

3.2 复化梯形写成独立子函数，权重向量只算一次

题目要求用复化梯形公式，逐次二分步长，再对外推。最直接的写法是三层循环：外层 y 网格、内层 x 网格、最外层二分。三层套在一起，改一个参数要重看一遍缩进，非常难受。我后来的做法是把它拆成三块：被积函数一个文件、复化梯形一个函数、主控脚本一个文件。函数签名把矩形区域和两个方向的等分数都显式传进去，谁调用谁负责。

function T = getT(f, a, b, c, d, nx, ny) % 复化梯形求二重积分 % 内层 x 在 [a,b] 上 nx 等分，外层 y 在 [c,d] 上 ny 等分 hx = (b - a) / nx; hy = (d - c) / ny; wx = ones(nx+1, 1); wx([1 end]) = 0.5; % 梯形权：端点 1/2，内部 1 wy = ones(ny+1, 1); wy([1 end]) = 0.5; T = 0; for i = 1:ny+1 y = c + (i-1)*hy; s = 0; for j = 1:nx+1 x = a + (j-1)*hx; s = s + wx(j) * f(x, y); end T = T + wy(i) * s; end T = T * hx * hy; % 两个方向的步长都要乘，漏一个结果差 h 倍 end

权重向量wx、wy是复化梯形的本质：内部节点权重 1，两个端点权重 1/2。把权重提出来单独构造，比在循环里写if j==1 || j==nx+1干净，而且不会在二分之后忘记更新端点判断。最后那一行hx * hy是最容易漏的乘子，漏掉的话结果会系统性偏大，而且偏大的倍数恰好等于步长乘积的倒数，二分几层之后数值会越来越离谱，反而不像「忘了乘」的样子。

被积函数单独一个文件，写成function z = myfun(x, y); z = tan(x.^2 + y.^2); end，这样换一道题只改这一个文件，主控脚本一个字不用动。

3.3 4/3 和 -1/3 这对系数怎么用、什么时候停

复化梯形的误差按 h² 展开，把步长减半前后的两个值 T_n 和 T_2n 做一次线性组合，就能把 h² 项消掉，剩下的主项是 h⁴ 阶，组合出来的正是复化 Simpson：

S_n = (4/3)·T_2n - (1/3)·T_n

系数怎么来的可以自己推一遍：设 T(h) = I + C·h² + …，那么 T(h/2) = I + C·h²/4 + …，要消掉 C·h²，需要 4·T(h/2) - T(h)，再除以 3 归一，就是上面的 4/3 和 -1/3。亲手推一次，比背公式牢靠，因为换成四点外推（消 h⁴ 项）的时候，你会知道系数该怎么改，而不是干瞪眼。

主控脚本的关键是把「上一次的外推值」和「这一次的外推值」存成两个变量，判据落在两者之差上：

f = @(x, y) tan(x.^2 + y.^2); a = 0; b = pi/3; c = 0; d = pi/6; eps = 1.2e-4; nx = 4; ny = 4; kmax = 12; % 层数上限，防死循环 Tprev = getT(f, a, b, c, d, nx, ny); nx = 2*nx; ny = 2*ny; Tcur = getT(f, a, b, c, d, nx, ny); Sprev = 4/3*Tcur - 1/3*Tprev; % 第一次外推 for k = 1:kmax Tprev = Tcur; nx = 2*nx; ny = 2*ny; Tcur = getT(f, a, b, c, d, nx, ny); Scur = 4/3*Tcur - 1/3*Tprev; if abs(Scur - Sprev) <= eps break; end Sprev = Scur; end fprintf('I = %.10f, 层数 = %d, nx = %d, ny = %d\n', Scur, k, nx, ny);

两个方向必须同时翻倍，只翻 x 不翻 y，收敛阶会被 y 方向拖住，看着也在变小，但达不到 h⁴ 的下降速度，往往多跑四五层还是差一点点。kmax这个保险很必要：浮点数的有效位就那么多，步长小到一定程度，舍入误差开始主导，差值反而不再下降，没有上限就是死循环。

3.4 双层循环慢在哪，向量化怎么改

原始写法跑起来确实慢，因为每次二分，节点总数变成四倍，而双层循环是纯标量运算。在 MATLAB 里这类循环的开销主要不是算术，是解释器逐次调度。改成矩阵一次算完，能快一到两个数量级：

xs = linspace(a, b, nx+1); ys = linspace(c, d, ny+1); [X, Y] = meshgrid(xs, ys); % 生成整张网格 Z = f(X, Y); % 一次调用算完所有节点 Wx = ones(1, nx+1); Wx([1 end]) = 0.5; Wy = ones(ny+1, 1); Wy([1 end]) = 0.5; T = Wy' * Z * Wx' * hx * hy; % 两次矩阵乘完成加权求和

这里的前提是被积函数写得支持数组输入，tan(x.^2 + y.^2)天然支持，只要注意点乘点幂。向量化之后单次求值的时间基本可以忽略，二分十几层也是秒级。当时作业里写的是循环版，交上去能跑对，但把 nx 提到几千就会明显卡顿，如果老师要求现场演示大网格，这一改就是生死差别。

4. 常微分方程四法对齐：RK4 起步、AB4 接棒、预测校正收口

4.1 先把初值问题摆正，再谈四种方法

题目给的初值问题是 y' = -x²y²，0 ≤ x ≤ 1.5，y(0) = 3，步长 h = 0.1，精确解 y(x) = 3/(1 + x³)。拿到手第一件事不是写 RK4，是验一下精确解对不对：求导得 y' = -9x²/(1+x³)²，而 y² = 9/(1+x³)²，两者乘上 -x² 完全吻合，y(0) = 3 也对。这步验算只需要两分钟，但能挡掉后面所有「误差表整体偏大」的怀疑。我吃过一次亏，精确解抄错一个系数，四种方法的误差全是同一个方向偏大，还以为是步长不够，把 h 改了三次都没用。

四个方法共用同一个右端函数和同一组节点，所以写代码的时候接口必须统一：都是[x, y] = method(f, x0, xN, y0, h)，返回完整节点和完整解向量。只有接口一致，后面才能用一行norm(y - yex(x))把四种方法摆进一张表。很多人每个方法各写一个脚本、各自打印，最后做对比时手工抄数，抄错一位就白忙。

4.2 RK4 的四段斜率和 AB4 的历史项缓存

RK4 是四阶单步法，每一步要在区间内取四个斜率加权：

function [x, y] = rk4(f, x0, xN, y0, h) x = (x0:h:xN)'; n = numel(x); y = zeros(n, 1); y(1) = y0; for i = 1:n-1 k1 = f(x(i), y(i)); k2 = f(x(i)+h/2, y(i) + h/2*k1); k3 = f(x(i)+h/2, y(i) + h/2*k2); k4 = f(x(i)+h, y(i) + h*k3); y(i+1) = y(i) + h/6 * (k1 + 2*k2 + 2*k3 + k4); end end

课本上四段斜率的排版经常挤在一起，k2 和 k3 的第二个参数都是 y(i) 加上半步长乘以斜率，k4 是整步长，而且 k4 的自变量是 x(i)+h 而不是 x(i)+h/2。这四处错任何一处，方法就掉阶，但在 h = 0.1 这种步长下，前几步的误差肉眼看不出来，只有把整条曲线和精确解画在一起才露馅。

AB4 是四阶显式多步法，用前面四个点的导数值外推下一步：

y_{i+1} = y_i + h/24 · [55f_i - 59f_{i-1} + 37f_{i-2} - 9f_{i-3}]

注意这里系数是 h/24，不是 h/12。有些资料写成 h/12 · [55, -59, 37, -9] 的等价位，那括号里就得整体除以 2，混着用就差了整整一倍。多步法需要四个历史点，第一个点以外的前三个点没有历史，只能由 RK4 提供同精度的起步值，这就是题目里「由 RK4 提供初值」的含义。

function [x, y] = ab4(f, x0, xN, y0, h) x = (x0:h:xN)'; n = numel(x); y = zeros(n, 1); [~, yr] = rk4(f, x0, x0+3*h, y0, h); % 只用 RK4 算前 4 个点 y(1:4) = yr(1:4); F = zeros(n, 1); F(1:4) = f(x(1:4), y(1:4)); % 先把已知点的导数值缓存下来 for i = 4:n-1 % 从第 5 个点开始由 AB4 推进 y(i+1) = y(i) + h/24 * (55*F(i) - 59*F(i-1) + 37*F(i-2) - 9*F(i-3)); F(i+1) = f(x(i+1), y(i+1)); % 边算边存，下一步还要用 end end

F这个导数缓存数组是多步法的核心数据结构，显式方法和校正方法都要靠它取历史。当年我第一次写的时候没缓存，每步重新调三次 f，结果数值上没错但速度慢了一半，而且因为重复计算浮点舍入不同，跟别人的结果在小数点后八位对不上，一度以为方法写错了。

4.3 预测校正与加权修正：251/270 和 19/270 到底在修正什么

AB4 和 AM4 各有短板：显式的 AB4 快但稳定性差，隐式的 AM4 精度高但每一步都要解方程。把两者拼起来就是预测校正：先用 AB4 算一个预测值，把它代进 AM4 的右端当作新点导数的估计，再算一次校正值。校正那一步的公式是：

y_{i+1}^{(c)} = y_i + h/24 · [9f(x_{i+1}, y_{i+1}^{(p)}) + 19f_i - 5f_{i-1} + f_{i-2}]

这一步把 AB4 的四阶精度继承了，同时用上了隐式公式的结构，误差比纯 AB4 小一个量级，代价只是每步多算一次 f。

再往上一层是带改进的版本，在预测值和校正值之间做一次加权：

y_{i+1} = (251/270)·y_{i+1}^{(c)} + (19/270)·y_{i+1}^{(p)}

两个权重加起来正好是 1，校正值占 93%，预测值只占 7%，所以这个加权不是简单平均，而是「以校正值为主、用预测值做微调」。它修的是局部截断误差里的低阶项，本质上是把两者的误差主项对齐后互相抵消。

function [x, y] = ab4am4(f, x0, xN, y0, h, mode) % mode = 1 预测校正；mode = 2 带改进的加权版本 x = (x0:h:xN)'; n = numel(x); y = zeros(n, 1); [~, yr] = rk4(f, x0, x0+3*h, y0, h); y(1:4) = yr(1:4); F = zeros(n, 1); F(1:4) = f(x(1:4), y(1:4)); for i = 4:n-1 yp = y(i) + h/24 * (55*F(i) - 59*F(i-1) + 37*F(i-2) - 9*F(i-3)); Fp = f(x(i+1), yp); % 预测点上的导数估计 yc = y(i) + h/24 * (9*Fp + 19*F(i) - 5*F(i-1) + F(i-2)); if mode == 2 y(i+1) = 251/270 * yc + 19/270 * yp; % 加权修正低阶误差项 else y(i+1) = yc; end F(i+1) = f(x(i+1), y(i+1)); end end

一个容易忽略的点：加权之后 y(i+1) 变了，那么 F(i+1) 必须用加权后的值重新算，不能用 Fp 代替。用 Fp 的话，下一步的历史导数就是预测值导数，误差会一代一代累积，跑十几步之后反而比不加权的版本还差。这类「改了结果忘了改缓存」的问题在多步法里特别常见，因为缓存是隐形的，代码看起来完全正常。

4.4 误差 2 范数怎么算，为什么排序比数值更重要

四个方法跑完，把解向量和精确解相减，取 2 范数：

f = @(x, y) -x.^2 .* y.^2; yex = @(x) 3 ./ (1 + x.^3); h = 0.1; x0 = 0; xN = 1.5; y0 = 3; [xx, y1] = rk4(f, x0, xN, y0, h); [~, y2] = ab4(f, x0, xN, y0, h); [~, y3] = ab4am4(f, x0, xN, y0, h, 1); [~, y4] = ab4am4(f, x0, xN, y0, h, 2); e = @(yv) sqrt(sum((yv - yex(xx)).^2)); % 未加权的 2 范数 fprintf('RK4=%g AB4=%g AB4-AM4=%g 改进=%g\n', e(y1), e(y2), e(y3), e(y4));

范数取法会改变数值大小，但不会改变四个方法的相对排序，因为它们是同一组节点上同一把尺子量出来的。这一点很关键：作业里要的结论是「哪种方法精度高」，不是某个具体小数，所以范数的加权方式只要在报告里写清楚，就不影响判断。反过来，如果排序和课本结论不一致，那一定是代码有问题，而不是范数选得不好看。

5. 踩坑与排查：五个让我重跑整晚的具体问题

5.1 样条曲线在端点附近翘起一块

现象：中间所有插值点都能对上，S(i+0.5) 打印出来也是平滑的一列数，但把曲线画出来，x 靠近 0 的那一段明显向上鼓，靠近 10 的那一段向下压。

原因：首行和尾行的公式方向写反了。首行应该是 6/h·(差商 - 端点导数)，尾行是 6/h·(端点导数 - 差商)，我把尾行也写成了差商在前，等于把 y_10' = 0.2 当成了 -0.2 在用。中间节点完全不参与首尾行，所以中间结果看着毫无异常。

解决：写完之后用数值差分反推两端斜率，(S(1e-3) - S(0))/1e-3应该等于 0.8，误差量级 1e-3 以内。对不上就回去检查首尾两行的减数被减数顺序，别去动中间的三对角。

5.2 AB4 从第 5 个点开始就偏，而且越跑越偏

现象：RK4 那条曲线和精确解基本重合，AB4 从第 5、6 个点开始肉眼可见地往上抬，跑到 x = 1.5 差了快一个百分点。

原因：起步点只用了 RK4 算前 3 个点，而 AB4 需要 i、i-1、i-2、i-3 四个历史导数，循环从 i = 3 开始，F(i-3)取到的是索引 0，MATLAB 里索引 0 直接报错；但如果改成从 i = 4 开始而只初始化了 y(1:3)，F(4) 就是初始化的零，f 在零处的值被当成真实历史导数用了一次，误差立刻被放大四倍进入主循环。

解决：起步必须算满四个点，y(1:4)用 RK4，F(1:4)立刻用 y(1:4) 填满，循环从i = 4开始。检查方法很简单，在循环第一轮打印 i 和 F(i-3:i)，看四个数是不是都非零、是不是对应前四个节点的导数值。

5.3 重积分二分停得太早，结果差一个量级

现象：判据用的是abs(Tcur - Tprev) <= eps，程序三五层就退出，结果和手算的参考值差挺多，把 eps 调小十倍还是差不多。

原因：复化梯形的收敛是 O(h²)，相邻两次梯形值之差本身就比较小，拿它当判据相当于放宽了标准；而要求的外推值 S 是 O(h⁴)，它的实际误差比梯形差值还小，用梯形的差值判断 S 是否收敛，等于提前收工。

解决：判据落在相邻两次外推值上，abs(Scur - Sprev) <= eps。同时保留一个层数上限，因为浮点精度有限，差值小到一定程度会开始抖动，没有上限就是死循环。

5.4 只翻倍了内层网格，收敛速度上不去

现象：每次都把 nx 翻倍，ny 保持不变，跑了七八层差值还在 1e-3 上下晃，怎么都进不了 1e-4。

原因：二重积分的误差是两个方向误差之和，外层 y 的步长没变，它的 O(hy²) 项一直挂在那里，内层再怎么细化也压不下去。表现出来就是收敛「卡住」，而不是「变慢」。

解决：两个方向的等分数同时翻倍，nx 和 ny 一起乘 2。判据里可以顺带打印 nx、ny，看到只有一个在变，就知道是这里的问题。

5.5 误差表整体偏大，但排序是对的

现象：四种方法的相对排序和课本结论一致，RK4 最好，AB4 最差，但每个数值都比同学的大了一截，差得还挺均匀。

原因：精确解写错了，或者 f 里的系数抄错了。这道题的精确解是 3/(1+x³)，如果误写成 3/(1+x²) 或者 3/(1+x)³，误差会整体抬升但不会改变方法之间的相对关系，所以排序看起来完全正常，特别有欺骗性。

解决：动手写代码之前，先把精确解代回微分方程验一遍，两边是否恒等、初值是否满足。这一步两分钟，能省掉后面几个小时的自查。另一种可能是 f 里的负号漏了，那误差会大到离谱，反而容易发现。

6. 误差排序读法：四种方法的结论怎么用在下一次选题上

把四种方法在 h = 0.1 下的 2 范数误差摆在一起，结论很干脆：

RK4 比 AB4 小两个数量级，这个差距不是靠调步长能追回来的，因为两者阶数相同但误差常数差得远，RK4 的误差常数在四阶方法里属于特别小的那一档。改进的 AB4-AM4 比不加权的版本好将近一倍，说明 251/270 和 19/270 这对权重确实在抵消低阶项，而不是摆设。AB4 垫底也在意料之中，显式多步法没有校正环节，误差只能靠缩小步长硬压。

真正值得记住的是这个排序背后的取舍。RK4 精度最高，代价是每步四次函数求值；AB4 每步只要一次，但精度差两个数量级；AB4-AM4 每步两次，精度介于两者之间。如果在实际计算里 f 的求值很贵（比如每次都要解一个子问题），那 AB4-AM4 的性价比就出来了，用两倍的计算量换二十倍的精度提升。如果 f 很便宜，直接用 RK4，省心。

回过头看，这道题真正教给我的不是四个公式，是「误差表要怎么读」。第一次做完，我盯着那四个数字只想着哪个最小，第二次重做才发现，改进版比不加权版好一倍这件事，恰好说明加权那一步不能省；而 AB4 单独用误差大这件事，恰好说明多步法必须配校正。每一行数字都对应一个设计决策，读不懂数字就等于没做这道题。

从那以后我每次交数值作业之前，都强制走一遍三件事：先把精确解代回原方程验一遍，再用最简单的特例（比如 h 取大一点、节点取少一点）手算一步对一下程序，最后把误差表按数量级排一遍，看排序是否符合方法阶数的预期。这三步加起来不到二十分钟，但救回来的通宵不止一次。希望帮到你。

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