简介：面向理工科数值分析课程的上机作业文档，围绕三次样条插值、重积分计算与常微分方程初值问题数值解三类实验，梳理公式推导、算法步骤和MATLAB编程实现要点，适合本科生、研究生对照课程题目完成上机练习与报告整理。压缩包内仅有1个docx文件，大小约566KB，篇幅虽小但内容集中，便于按上机题号快速定位样条插值、重积分和初值问题三个模块。已有904人学习，文档结合具体算例给出车门曲线三次样条插值、复化梯形公式配合逐次二分与外推求二重积分、RK4与AB4及AB4-AM4预测校正法的计算过程和结果比较。读者可从中获得各题算法框架、MATLAB实现思路、误差对比与精度排序结论，也能借鉴作者对边界条件、步长选择和程序调试的反思，用于查漏补缺和上机报告写作。
1. 数值分析上机题的三条主线：插值、积分与微分方程
数值分析上机题给人的错觉是“就三道编程题”，实际它把三门课的核心工具压进了一条线：插值负责把一堆离散点造成连续函数，数值积分负责把算不出的定积分变成能累加的项，常微分方程数值解负责把连续的导数关系变成一步步往前推的递推式。三块共用同一套底子——误差阶、步长、收敛判据、稳定性，所以按这条线复习，比按教材章节顺序背公式省力得多。
这份笔记按三次样条插值、重积分计算、初值问题数值解三块展开，中间插入一圈最容易翻车的位置和补法，末尾给一份可以直接照着过的自查清单。适合正在赶课程上机、准备期末或者考研复试题的人对照着查漏，也适合写完程序跑不出结果、想找人对着公式核一遍的人。
2. 三次样条插值：从三弯矩方程到 S(i+0.5) 的完整链路
三次样条是这三道题里“公式最长、但程序最短”的一道。它的难点不在写代码，而在于把 4n 个待定系数消成 n+1 个节点二阶导 M，再把边界条件塞进这个三对角系统里。理解这条消元路径，比死记任何一种写法都管用。
2.1 三弯矩方程：把 4n 个系数压成 n+1 个 M
先把未知量数清楚。整个区间被分成 n 个子区间，每个子区间上是一个三次多项式，4 个系数，一共 4n 个未知量。约束来自哪里：每个节点函数值相等给 2(n-1) 个，内节点一阶导连续给 n-1 个，二阶导连续给 n-1 个，再加两个边界条件，正好 4n 个。直接解 4n 元方程组在数值上很蠢，工程上的标准做法是改用节点处的二阶导 M_j = S''(x_j) 作为未知量，因为三次样条的二阶导是分段线性的，只要 M 定了，整条曲线就定了。
得到的三弯矩方程是这一章的核心：
μ_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}] （二阶差商乘 6）这里 f[x_{j-1}, x_j, x_{j+1}] 是二阶差商，展开成步长形式就是 ((y_{j+1}-y_j)/h_j - (y_j-y_{j-1})/h_{j-1}) / (h_{j-1}+h_j)。均匀网格 h 全相等时，它退化成 (y_{j-1} - 2y_j + y_{j+1}) / h²，乘 6 之后正好是 6(y_{j-1} - 2y_j + y_{j+1})/h²，这个退化形式最好记，也最容易在考场上临时推。
边界条件决定了系统怎么闭合。第一型给两端的一阶导，方程是：
2·M_0 + M_1 = (6/h_0) · ( f[x_0,x_1] - y_0' ) M_{n-1} + 2·M_n = (6/h_{n-1}) · ( y_n' - f[x_{n-1},x_n] )第二型给两端的二阶导，直接 M_0 = y_0''、M_n = y_n''，系统降成 n-1 阶。第三型是周期条件，首尾方程要联立。三种边界条件的差别值得单独列一张表记：
| 边界类型 | 已知量 | 方程处理 | 系统阶数 |
|---|---|---|---|
| 第一型（固支） | y_0'、y_n' | 补两个含 M_0、M_n 的方程，系数为 2 和 1 | n+1 |
| 第二型（自然） | y_0''、y_n'' | M_0、M_n 直接赋值，不参与求解 | n-1 |
| 第三型（周期） | y_0 = y_n、导数周期 | 首尾方程合并，矩阵变成循环三对角 | n |
重点在哪：均匀网格下 μ_j = λ_j = 0.5，对角元是 2，非对角元之和是 1，矩阵严格对角占优。这条性质直接保证追赶法不用选主元、误差不被放大，程序里可以放心地写两趟递推而不用做列主元消去。理解了这一点，才算真懂为什么教材非要绕到 M 上去。
2.2 车门曲线型值点：h、d 与右端项的逐项对照
拿题干给的那组车门曲线数据走一遍。x_i = 0, 1, …, 10，一共 11 个点，n = 10 个子区间，h 全部等于 1，所以 μ_j = λ_j = 0.5。由于 h = 1，d_j = 6(y_{j-1} - 2y_j + y_{j+1})，可以直接心算核对：
| j | y_{j-1} | y_j | y_{j+1} | y_{j-1}-2y_j+y_{j+1} | d_j |
|---|---|---|---|---|---|
| 1 | 2.51 | 3.30 | 4.04 | -0.05 | -0.30 |
| 2 | 3.30 | 4.04 | 4.70 | -0.08 | -0.48 |
| 3 | 4.04 | 4.70 | 5.22 | -0.14 | -0.84 |
| 4 | 4.70 | 5.22 | 5.54 | -0.20 | -1.20 |
| 5 | 5.22 | 5.54 | 5.78 | -0.08 | -0.48 |
| 6 | 5.54 | 5.78 | 5.40 | -0.62 | -3.72 |
| 7 | 5.78 | 5.40 | 5.57 | 0.55 | 3.30 |
| 8 | 5.40 | 5.57 | 5.70 | -0.04 | -0.24 |
| 9 | 5.57 | 5.70 | 5.80 | -0.03 | -0.18 |
第 6 个点附近 d_6 = -3.72、d_7 = +3.30 这两项绝对值突然变大，对应的是 y 在 i=6 到 i=7 之间从 5.78 掉到 5.40 那个折点——数据本身不光滑，插值函数在这一段必然的曲率变化就大。反过来看，如果程序跑出来 M_6、M_7 和别处差不多，那多半是 d 数组下标错位了，这是一个非常好的自查信号。
边界条件用第一型，端点导数是 y_0' = 0.8、y_10' = 0.2。h_0 = h_9 = 1，于是：
左端：2M_0 + M_1 = 6 × ((3.30 - 2.51)/1 - 0.8) = 6 × (-0.01) = -0.06 右端：M_9 + 2M_10 = 6 × (0.2 - (5.80 - 5.70)/1) = 6 × 0.10 = 0.60合起来是 11 阶三对角系统，第一行和最后一行系数是 (2, 1) 和 (1, 2)，中间九行都是 (0.5, 2, 0.5)。解出 M 之后，第 j 个区间上的插值函数写成：
S(x) = M_j · (x_{j+1} - x)^3 / (6h_j) + M_{j+1} · (x - x_j)^3 / (6h_j) + (y_j - M_j · h_j²/6) · (x_{j+1} - x) / h_j + (y_{j+1} - M_{j+1}·h_j²/6) · (x - x_j) / h_j题目要打印 S(i+0.5)，i = 0…9。注意 x = 0.5 落在第 0 个区间，x = 9.5 落在第 9 个区间，每个半点都要先判断落在哪个子区间，再代对应那一段的四项式。这是最容易写错的地方之一：直接拿一个统一公式套全区间，结果只会在节点上对，别处全歪。
2.3 程序骨架：组装—求解—回代三段式
三段式结构最适合这道题，也最好查错。第一段根据 x、y 和两个端点导数组装三对角矩阵，第二段解方程得 M，第三段逐点回代求 S。三段各自可以单独打印中间量核对。
function M = splineM3(x, y, dy0, dyn) % 第一型三次样条：求节点二阶导 M（三弯矩方程） n = length(x) - 1; % 子区间个数 h = diff(x); % h(j) = x(j+1) - x(j) A = zeros(n+1); rhs = zeros(n+1, 1); % 左边界：2*M(1) + M(2) = 6/h(1) * (f[x0,x1] - dy0) A(1,1) = 2; A(1,2) = 1; rhs(1) = 6 / h(1) * ((y(2) - y(1)) / h(1) - dy0); % 中间节点：mu*M(j-1) + 2*M(j) + lam*M(j+1) = d(j) for j = 2:n mu = h(j-1) / (h(j-1) + h(j)); lam = 1 - mu; d = 6 * ((y(j+1)-y(j))/h(j) - (y(j)-y(j-1))/h(j-1)) / (h(j-1)+h(j)); A(j, j-1) = mu; A(j, j) = 2; A(j, j+1) = lam; rhs(j) = d; end % 右边界：M(n) + 2*M(n+1) = 6/h(n) * (dyn - f[x_{n-1},x_n]) A(n+1, n) = 1; A(n+1, n+1) = 2; rhs(n+1) = 6 / h(n) * (dyn - (y(n+1) - y(n)) / h(n)); M = A \ rhs; % 规模大时换成追赶法，可省一个量级的内存 end这段代码有两个索引细节必须说清。MATLAB 下标从 1 开始，而公式里的 j 从 1 到 n-1 对应的是 MATLAB 的 2 到 n，所以循环写成for j = 2:n，每一行填的是 A(j, j-1)、A(j, j)、A(j, j+1)；h由diff得到，长度是 n，h(j-1) 和 h(j) 都合法。参数 dy0、dyn 是端点一阶导，传进来之前先确认单位与 x、y 一致，别拿角度值往弧度制公式里塞。
求值部分单独写一个函数，把“落在哪个子区间”这件事显式处理掉：
function S = splineEval(x, y, M, xi) % 在 xi 各点上回代求 S(xi)，xi 可以是向量 S = zeros(size(xi)); for k = 1:numel(xi) j = find(x <= xi(k), 1, 'last'); % 定位所属子区间 j = min(j, length(x) - 1); % 右端点归到最后一个区间 hj = x(j+1) - x(j); S(k) = M(j) * (x(j+1) - xi(k))^3 / (6*hj) ... + M(j+1) * (xi(k) - x(j))^3 / (6*hj) ... + (y(j) - M(j) *hj^2/6) * (x(j+1) - xi(k)) / hj ... + (y(j+1) - M(j+1)*hj^2/6) * (xi(k) - x(j)) / hj; end endfind(x <= xi(k), 1, 'last')拿到的是不大于 xi 的最后一个节点下标，xi 恰好等于某个节点时会落到该节点右侧的区间上，所以后面要用min(j, length(x)-1)兜一下右端点，防止 x = x_n 时 j 越界。写完这两个函数，正常流程是：先用splineEval在全部节点上验一遍 S(x_i) 是否等于 y_i（误差应在 1e-12 量级），再用中心差分估一下两端的 S'，看是否贴近 0.8 和 0.2。这两个检查做完，程序基本就不会有隐藏错误了。
3. 重积分数值计算：复化梯形 + 逐次二分 + 外推
重积分这道题的题面很长，拆开只有三件事：把二重积分化成两层定积分、每层用复化梯形算、用逐次二分加 Richardson 外推把精度顶上去。真正的坑都在“什么时候停”和“怎么少算点”上。
3.1 两层嵌套：内层对 x、外层对 y 的复化梯形
题目给的积分是 I(f) = ∫_c^d ( ∫_a^b f(x,y) dx ) dy。括号里那层先算：固定 y，对 x 在 [a,b] 上做复化梯形，得到一个只跟 y 有关的函数 g(y)；再对 g(y) 在 [c,d] 上做一次复化梯形。设 x 方向 n 等分、步长 h = (b-a)/n，y 方向 m 等分、步长 k = (d-c)/m，二重复化梯形的完整形式是：
T(n,m) = (k/2) · Σ_{i=0}^{m} w_i · g(y_i) g(y_i) = (h/2) · [ f(a,y_i) + f(b,y_i) + 2·Σ_{j=1}^{n-1} f(x_j, y_i) ] w_i = 1 (i = 0 或 m)，w_i = 2 (其余)把内层展开写全，就是四个部分：四个角点各算一次、四条边上的内部点各算两次、内部网格点各算四次。写程序时不需要真的按“角、边、内”分类，直接用权重因子做外积最省事，也最不容易漏项。
重点在于理解“误差是两层叠加”的。如果内层固定用很小的 n、外层死命加密 m，总误差会被内层的 O(h²) 卡住，外层加得再细也没用。所以逐次二分必须是两个方向同时二分，n 和 m 一起翻倍，让 h² 和 k² 同步下降；或者强制内层精度远高于外层，比如内层容差取外层容差的百分之一。这一条是这道题最容易被忽略的设计原则。
3.2 Richardson 外推：一个公式，两种身份
外推的出发点是复化梯形的误差展开：
T_n = I + c·h² + O(h⁴) T_2n = I + c·(h/2)² + O(h⁴) = I + c·h²/4 + O(h⁴)两式一联立，把 c·h² 这一项消掉：
S_n = (4/3)·T_2n - (1/3)·T_n = I + O(h⁴)这个组合式有两个身份：从外推的角度看，它是 Richardson 外推一次的结果；从求积公式的角度看，它恰好就是复化 Simpson 公式。也就是说，“用两次梯形值线性组合提高一阶精度”和“直接上 Simpson 公式”在数学上是同一件事，只是一个从误差展开发推导、一个从插值多项式推。把这两条线对上，就不需要死记 Simpson 的系数了。
需要注意的方向：外推要求步长严格按 h、h/2、h/4 的序列二分，两次值必须对应同一积分区间，且 h 的比例正好是 2。如果中途改了区间或者跳着二分，系数 4/3 和 -1/3 就不成立，这是手工验算时最常见的对不上的原因。
还有一个副产品：由 T_n = I + c·h² 可推出 I - T_n ≈ (T_2n - T_n)/3，也就是说两次梯形值之差的三分之一就是梯形值本身的误差估计。这个式子在做停机判据时非常好用，可以顺便打印出来看收敛速度是不是二阶。
3.3 逐次二分的复用技巧与被积函数拆分
朴素写法是每次二分都重新算全部网格点，工作量直接翻四倍（二维情况下 n、m 各翻倍）。实际上二分之后老网格点全部保留，只需要算新增的中点，一次二分的新增点数约等于原来的点数，而不是四倍。这个技巧写出来只有一行：
function T2 = refineTrapz(f, a, b, n, Tn) % 把 [a,b] 上 n 等分的复化梯形值 Tn 加密到 2n 等分 % f 是函数句柄，n 是当前等分数，Tn 是当前的梯形近似值 h = (b - a) / n; xnew = a + h/2 : h : b - h/2; % 只枚举新增的中点，共 n 个 T2 = Tn/2 + h/2 * sum(f(xnew)); % 老点权重减半 + 新点全权重 end参数含义：a、b是积分上下限，n是加密前的等分数（必须是 2 的幂次才好一路二分），Tn是加密前的梯形值，f是函数句柄。逻辑上就是把原来的 T_n = (h/2)(f_0 + 2Σf_old + f_n) 改写成 T_2n = (h/4)(f_0 + 2Σf_old + f_n) + (h/2)Σf_new，前一块正好是 T_n/2，后一块是新增中点的一次求和。这段代码在二重积分里要套两层：内层对每个 y_i 调用一次，外层再对 g(y) 做一次。
被积函数和积分器一定要拆开。题目里的被积函数是 tan(x² + y²)，把它单独写成一个myfun.m，积分器只接收函数句柄，这样换一道题只需要换这一个文件，不用动主程序。同理，计算梯形值的部分写成getT.m，主程序exp5.m只负责组织循环、判断收敛、打印结果。这种拆法在上机时不显眼，但调试时价值很大：怀疑积分器有问题就换成 f = 1 去积，解析结果是区间面积，一测就知道；怀疑被积函数有问题就单独画图看有没有奇点。
3.4 停机判据与运行速度：慢在哪、怎么省
题干要求的停机条件是“相邻两次近似值之差的绝对值不超过 ε”，ε 取 0.5×10⁻⁵。这里有个口径问题必须先定死：判据是加在外推后的 S 上，还是加在二分前后的 T 上？两种写法收敛速度不同，最终迭代次数也不同。稳妥的做法是两者都监控，主判据用 |S_new - S_old| < ε，同时打印 |T_2n - T_n|/3 作为误差估计，两个量如果数量级差太多，说明还没进入渐近收敛区，继续二分。
第二个口径问题是内层。二重积分每算一个外层点都要做一次完整的内层积分，如果内层只固定一个 n，总误差就被锁死了。合理的做法是内外层同步二分，让 n 和 m 一起翻倍，外推也在两个方向上分别做或一次性做二维外推。
运行速度慢是这道题的常见反馈，原因基本能归成三条。第一，tan 的调用次数太多，二维逐次二分下总调用量按 4 倍增长，几十万次调用在解释型语言里开销很直观。第二，循环写成了标量循环，MATLAB 里对 x 方向的求和完全可以向量化成矩阵运算，把for j换成对矩阵列求和，速度差一个量级。第三，函数句柄在循环里被反复解析，把f(x,y)提前算成矩阵缓存起来能省不少。真要在报告里写“程序慢”，最好给出调用次数和耗时两组数，而不是一句“运行速度慢是很大的不足”——后者是感受，前者是可供改进的依据。
自测时可以用一个已知解析解的积分验证：把被积函数换成 1，积分值应等于 (b-a)(d-c)；换成 x·y，积分值应等于 (b²-a²)(d²-c²)/4。这两个测试能在十行代码内定位“是积分器错了还是被积函数写错了”。
4. 常微分方程初值问题：RK4 与 Adams 系列的四个程序
这道题要求写四个程序、跑同一道初值问题、和精确解比误差，最后回答“能得到哪些结论”。四个方法的公式都能查到，真正考的是：为什么同是四阶，误差能差两个数量级；以及预测校正到底在补什么。
题目给的初值问题是 y' = -x²y²，0 ≤ x ≤ 1.5，y(0) = 3，步长 h = 0.1，精确解 y = 3/(1 + x³)。先验一下：对 y = 3(1+x³)⁻¹ 求导得 y' = -9x²/(1+x³)²，而 -x²y² = -x²·9/(1+x³)²，两边一致，说明精确解没抄错。这一步别省，精确解写错的话后面四个方法的误差全是错的，反而看不出问题。
4.1 RK4：四个斜率与 h/2 的取点方式
四阶 Runge-Kutta 的标准形式是：
y_{i+1} = y_i + (h/6)·(k1 + 2k2 + 2k3 + k4) 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)四个斜率的取点顺序是：起点、中点（用 k1 推）、中点（用 k2 推）、终点（用 k3 推）。权重是 1:2:2:1，除 6。这里的记忆锚点是“两个中点，权重各 2”，因为 Simpson 公式的权重也是 1:4:1 的变体，两者同源，都是从积分节点权重来的。
function [x, y] = rk4(f, a, b, y0, h) % 定步长四阶 Runge-Kutta x = (a:h:b)'; n = length(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参数说明：f是右端函数句柄，a、b是区间端点，y0是初值，h是步长。x = (a:h:b)'生成列向量，长度 n 由步长决定，如果 (b-a)/h 不是整数，最后一个点会被截掉，所以上机前先确认 1.5/0.1 = 15 能整除。RK4 是单步法，只需要当前点的信息，天然自启动，这也是它能给后面几个多步法供初值的原因。
RK4 的局部截断误差是 O(h⁵)，全局是 O(h⁴)。这意味着步长减半，误差大约降到 1/16。这个比例可以作为验算工具：用 h = 0.1 和 h = 0.05 各跑一遍，看误差比是不是接近 16，偏离太远说明程序里某处系数写错了。
4.2 AB4 的启动问题：为什么必须用 RK4 供初值
四阶 Adams-Bashforth（AB4）是显式多步法：
y_{i+1} = y_i + (h/24)·[55·f_i - 59·f_{i-1} + 37·f_{i-2} - 9·f_{i-3}]注意这里的系数是 55、-59、37、-9，配的是 h/24，不是 h/12。这四个系数来自对 f 在最近四个点上做三次插值再积分，交错的正负号是插值多项式系数的直接结果，符号写反一个，结果会明显偏离。
AB4 的麻烦在于它需要前四个点的 f 值才能启动，也就是 f_i、f_{i-1}、f_{i-2}、f_{i-3} 都得已知。第一个点 i = 0 时根本没有前三个点，所以必须用单步法（题面指定 RK4）先把 y_1、y_2、y_3 算出来，再从 i = 3 开始跑 AB4 循环。程序里循环起点写成for i = 4:n-1（MATLAB 下标），一旦写成从 1 开始，f(i-3) 立刻越界，报错还算好的，某些语言里会静默读到脏数据。
这里有个容易被忽略的连带影响：启动段的误差会一路传下去。AB4 的误差常数本身就比 RK4 大，再加上启动用的是 RK4 的近似值而不是精确值，整体误差自然比 RK4 高出一截。等会儿看误差表就明白了，AB4 的 2-范数误差是 4.8×10⁻³，比 RK4 的 3.8×10⁻⁵ 大了两个数量级。
4.3 AB4-AM4 预测校正与 251/270 加权
四阶 Adams-Moulton（AM4）是隐式公式：
y_{i+1} = y_i + (h/24)·[9·f(x_{i+1}, y_{i+1}) + 19·f_i - 5·f_{i-1} + f_{i-2}]右端含未知的 y_{i+1}，不能直接算。解决办法是拿 AB4 先“预测”一个 y_{i+1}^{(p)}，代入 AM4 右端算出“校正”值 y_{i+1}^{(c)}，这就是预测校正系统。因为只迭代一次、不反复收敛，所以它严格来说不是隐式解法，而是显式化的近似。
带改进的版本再做一步加权平均：
y_{i+1} = (251/270)·y_{i+1}^{(c)} + (19/270)·y_{i+1}^{(p)}两个权重加起来正好等于 1，这是检验有没有抄错系数的最快办法。改进的动机是：预测值和校正值各自带有不同符号、不同大小的误差，按固定比例混合可以把主误差项再抵消掉一部分，代价几乎为零。
% 前提：f 已按 x 网格预先算好存入向量 fv，y(1:4) 由 RK4 提供 for i = 4:n-1 % AB4 预测 yp = y(i) + h/24*(55*fv(i) - 59*fv(i-1) + 37*fv(i-2) - 9*fv(i-3)); % AM4 校正：右端第一项要用预测值处的斜率 fc = f(x(i+1), yp); yc = y(i) + h/24*(9*fc + 19*fv(i) - 5*fv(i-1) + fv(i-2)); % 改进加权 y(i+1) = 251/270*yc + 19/270*yp; fv(i+1) = f(x(i+1), y(i+1)); % 为下一步准备斜率 end关键在fc = f(x(i+1), yp)这一行：AM4 用的是校正前的预测值处的斜率，不是 y(i) 处的，也不是校正后的。写成f(x(i+1), y(i+1))会变成隐式迭代的味道，但只迭代一次并不收敛，结果反而更差。另外fv(i+1)必须用最终的加权值重算，不能沿用 fc，否则下一步的 AB4 就用了错的斜率。
4.4 误差对比：四种方法在同一个步长下的排序与原因
把四种方法在同一网格上跑完，与精确解逐点比较，取 2-范数误差，题面给出的结果是：
| 方法 | 误差（2-范数） | 精度排序 | 每步右端函数调用次数 |
|---|---|---|---|
| RK4 | 3.8×10⁻⁵ | 1 | 4 |
| AB4 | 4.809×10⁻³ | 4 | 1 |
| AB4-AM4 | 9.20×10⁻⁴ | 3 | 2（含预测） |
| 改进的 AB4-AM4 | 5.31×10⁻⁴ | 2 | 2（含预测） |
排序是 RK4 > 改进 AB4-AM4 > AB4-AM4 > AB4，这个结论本身好记，但报告里只写排序是不够的，得解释为什么。
第一，四个方法的阶数其实都是四阶，理论收敛速度一样，差距全在误差常数上。RK4 的误差常数在这类问题上最小，所以同样的 h 下误差最小。第二，AB4 是显式外推，稳定性区域比 RK4 窄，误差常数也大，这两个因素叠加，导致它在 h = 0.1 这种“不算小”的步长下误差明显偏大。第三，AM4 的隐式结构本身误差常数更小，但只做一次预测校正并没有完全拿到隐式的好处，所以落在中间；改进版把预测值和校正值按 251:19 混合，等于又消掉了一部分主误差项，于是比 AB4-AM4 再好一档。第四， RK4 每步要算四次 f，AB4 只要一次，所以如果按“达到同样精度所需的函数调用次数”来评价，多步法的优势才会显出来——这也是这类题最后那道“能得到哪些结论”该答的方向：精度和计算量要一起看，不能只看误差表的排序。
还有一个容易被忽略的观察点：把 h 减半再跑一遍，四种方法的误差应该都按 1/16 左右下降，如果某一种下降得明显更慢，说明它的启动段或者索引处理有问题。这个二次验证比反复读代码有效得多。
5. 易错点排查：五个最容易翻车的位置与补法
公式都能查到，程序也都能跑通，最后分数差在细节上。下面这五条是这类上机题里反复出现的翻车点，按“公式层—索引层—判据层”三层排下来，每条都写清漏在哪、为什么漏、怎么补。
5.1 公式层：符号镜像与端点符号
易漏点一：三弯矩方程里 μ_j 和 λ_j 配错了 M。标准形式是 μ_j 乘 M_{j-1}、λ_j 乘 M_{j+1}，而 μ_j 只跟左边的步长 h_{j-1} 有关，λ_j 只跟右边的 h_j 有关。
为什么漏：这两个希腊字母在教材里经常被反过来定义，有的书把 μ 定义为右步长占比，有的定义为左步长占比，抄公式时只看了一眼字母位置就下笔，很容易镜像。均匀网格下 μ = λ = 0.5，写反了完全看不出来，只有非均匀网格才会暴露。
怎么补：组装完矩阵后加一行断言，检查每行非对角元之和是否等于 1（第一型边界行除外）。再检查对称性之外的一个特征：把 x 数组整体反转、y 同步反转，解出的 M 序列应该也反转，用这个不变量测一次，比逐行核公式快。
易漏点二：第一型边界右端项的符号。左端是 (6/h_0)(f[x_0,x_1] - y_0')，右端是 (6/h_{n-1})(y_n' - f[x_{n-1},x_n])，两边的减数位置是反的。
为什么漏：两端的推导方向相反，左端是“用差商去逼近 y_0'”，右端是“用 y_n' 去逼近差商”，直觉上容易写成同一个顺序。
怎么补：用一个最简单的算例验——取 y = x²，端点导数 0 和 2，插值应该精确复现二次函数，此时 M 全是常数 2。如果符号写反，M 会在两端出现明显偏离，一眼就能看出。这个测试三十秒能跑
本文还有配套的精品资源，点击获取
