MATLAB数值分析上机实战：三次样条、二重积分外推与ODE方法对比

简介：面向学习数值分析与科学计算的高校学生，这份东南大学出版配套上机题文档围绕MATLAB编程实践展开，聚焦三次样条插值函数、重积分计算和常微分方程初值问题数值解三类典型题目。文档不仅给出第一型三次样条插值、复化梯形公式与外推、RK4、AB4及AB4-AM4预测校正等方法的算法思路，还结合车门曲线型值点、二重积分和给定初值问题展示计算过程与结果对比，并附有精度排序和编程反思，便于读者理解公式如何落地为程序、误差从何而来。文档同时强调主程序与子函数分离的组织方式，提到逐次二分步长、外推提高精度以及公式输入错误引发的排错经历，对课程作业自查和MATLAB实现复盘较有参考价值。资源包仅含1个docx文件，约566KB，篇幅紧凑，适合上机前预习与课后对照。目前已有904人学习，可作为数值分析上机训练的实用样例。

1. 三道上机题，卡人的从来不是公式

数值分析上机题最容易给人错觉：三次样条的三弯矩方程、复化梯形加外推、RK4 和 AB4-AM4 的系数，公式都在书上印着，抄下来似乎就完事。真坐到 MATLAB 前面，翻车的往往不是推导能力，而是把 μ 和 λ 的位置写反、把h/24抄成h/12、把内层外层积分变量搞混——这些地方不报错，只给你一条看着挺光滑、其实整体偏了的曲线。

这份数值分析上机题资料里连着做了三道题：三次样条插值（车门曲线）、二重积分（复化梯形逐次二分加外推）、常微分方程初值问题（RK4、AB4、AB4-AM4 和改进版对比）。下面按这三道题的顺序拆，把每道题的算法落点、可抄的 MATLAB 骨架、以及我实测过会踩的坑都摆出来。新手可以照着把程序跑通，熟手能对一下判断依据；做课程作业的、临时要补 MATLAB 数值程序的，都能直接拿走。

2. 三次样条插值：三弯矩方程怎么补边界、怎么解三对角

三次样条这道题的价值不在结果，在于它把“公式推导 → 线性方程组 → 程序实现 → 结果回代”四件事完整走了一遍。第一型（给定两端一阶导数）的三次样条，本质就是解一个三对角线性方程组求二阶导数 M，再用 M 把每个小区间上的三次多项式写出来。思路清楚了，代码不到五十行；思路没理清，改一天也调不出来。

2.1 型值点先核对一遍，再谈公式

车门曲线的型值点一共 11 个，x 从 0 到 10 等距，步长 h=1。资料里那行数据在复制、截图、转 Word 的过程中很容易粘连：2.51和3.30挨在一起，一不留神就变成13.30，而后面 4.04、4.70、5.22 一路平滑上去，唯独第二个点跳到 13，插值出来必然是一段剧烈振荡。我一般会先把数据摆成表，肉眼扫一遍量级，再往程序里灌。

边界条件是两端一阶导数已知：左端 y'(0)=0.8，右端 y'(10)=0.2。这两个数是第一型（Clamped）条件的全部输入，少了任何一个，方程组就欠两行，M 解不出来。等距情况下 h_j 全等于 1，这件事后面能用来手算校验，别浪费。

2.2 三弯矩方程的三个系数，缺一个都解不出

对内部节点 j=1,…,n-1，三弯矩方程是：

μ_j·M_{j-1} + 2·M_j + λ_j·M_{j+1} = d_j

其中 h_j = x_{j+1} - x_j，μ_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 倍二阶差商。

资料里把 λ 写成了 α，这是符号习惯问题，不影响计算，但要清楚它对应的是 M_{j+1} 的系数，不是 M_{j-1} 的。μ 和 λ 一个挂在左邻、一个挂在右邻，写反了方程组依然能解、程序依然不报错，只是解出来的 M 序列会整体错位，画出来的曲线在节点处一阶导不连续——肉眼不细看根本发现不了，这是这道题最阴的地方。

二阶差商我习惯显式写出来，避免嵌套调用出错：

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=1 的时候，这个式子退化成 (y_{j+1} - 2y_j + y_{j-1})/2，可以直接口算几个 d_j 来对程序输出。这是我最常用的一招：先手算 j=1 那一行的 μ、λ、d 三个数，跟程序打印的第一行对一遍，对不上就别往下跑了，越跑越乱。

2.3 第一型边界条件怎么补两行方程

内部的 n-1 个方程不够，一共 n+1 个未知量 M_0 到 M_n，还差两行，由端点导数条件补。第一型的补法是：

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] )

注意两行里差商和导数的先后顺序是反的：左端是“差商减导数”，右端是“导数减差商”。这个顺序反了，误差不会暴露成崩溃，只会让整条曲线两端微微翘起。第二型（给定二阶导数）的补法不一样，是 2M_0 + M_1 = 6·f[x_0,x_1,x_2] 这种形式；第三型（周期性）更不一样。三类型别在同一个程序里混用，是最常见的翻车点之一，做题前先确认题目给的是哪一型。

拼起来就是一个 (n+1)×(n+1) 的三对角方程组 A·M = d，直接用 MATLAB 的反斜杠A\d能解，但既然题目要求“编制通用程序”，把追赶法（Thomas 算法）写出来更符合要求，而且 O(n) 的复杂度在节点多的时候优势明显。

2.4 通用程序的骨架：构造—追赶—回代

先写主函数，负责组装三对角矩阵、调用追赶法、返回 M：

function [M, h] = spline3_clamped(x, y, dy0, dyn) % 第一型三次样条：输入型值点 x,y 与两端一阶导数 dy0,dyn % 输出各节点二阶导数 M（长度与 x 相同）和区间步长 h n = length(x) - 1; % 子区间个数 h = diff(x); % h(j) = x(j+1) - x(j) A = zeros(n+1); d = zeros(n+1,1); % 左边界：2*M0 + M1 = 6/h0 * (f[x0,x1] - y0') A(1,1) = 2; A(1,2) = 1; d(1) = 6/h(1) * ((y(2)-y(1))/h(1) - dy0); % 内部节点：mu*M(j-1) + 2*M(j) + lam*M(j+1) = 6*f[x(j-1),x(j),x(j+1)] for j = 2:n mu = h(j-1) / (h(j-1) + h(j)); lam = 1 - mu; A(j,j-1) = mu; A(j,j) = 2; A(j,j+1) = lam; d(j) = 6 * ((y(j+1)-y(j))/h(j) - (y(j)-y(j-1))/h(j-1)) ... / (h(j-1) + h(j)); end % 右边界：M(n-1) + 2*Mn = 6/h(n) * (yn' - f[x(n-1),x(n)]) A(n+1,n) = 1; A(n+1,n+1) = 2; d(n+1) = 6/h(n) * (dyn - (y(n+1)-y(n))/h(n)); M = thomas(A, d); % 追赶法求解，见下 end

追赶法把三对角矩阵压成三个向量，前代消元、回代求解，注意第一行和最后一行的消元边界：

function x = thomas(A, b) % 三对角方程组追赶法：A 为方阵，只取三条对角线 n = length(b); a = diag(A,-1); bb = diag(A,0); c = diag(A,1); for i = 2:n m = a(i-1) / bb(i-1); % 消元因子 bb(i) = bb(i) - m * c(i-1); b(i) = b(i) - m * b(i-1); end x = zeros(n,1); x(n) = b(n) / bb(n); % 回代 for i = n-1:-1:1 x(i) = (b(i) - c(i)*x(i+1)) / bb(i); end end

diag(A,-1)取的是下标 (i, i-1)，对应 a(i-1)；diag(A,1)取的是 (i, i+1)，对应 c(i)。这两个索引偏移错一位，追赶法解出来的结果会差得很离谱，但矩阵本身仍然可解，不报错。所以写完先拿max(abs(A*x - b))验残差，残差在 1e-12 量级才算通过。

2.5 打印 S(i+0.5) 的定位技巧

题目要打印 S(0.5)、S(1.5)……S(9.5)。插值函数是分段定义的，得先判断查询点落在哪个子区间。我用find(x <= xq, 1, 'last')定位左端点下标，再做上下界钳位，防止 xq 恰好等于最后一个节点时越界：

function v = eval_spline(x, y, M, xq) % 用二阶导数 M 回代计算三次样条在 xq 处的值 h = diff(x); v = zeros(size(xq)); for k = 1:numel(xq) j = find(x <= xq(k), 1, 'last'); j = min(max(j, 1), length(x)-1); % 钳位，防止越界 hj = h(j); t1 = (x(j+1) - xq(k)) / hj; % 左基函数 t2 = (xq(k) - x(j)) / hj; % 右基函数 v(k) = t1*y(j) + t2*y(j+1) ... + ((t1^3 - t1)*M(j) + (t2^3 - t2)*M(j+1)) * hj^2 / 6; end end

这个式子是三次样条的标准回代公式，(t^3 - t)那一项乘的是 M 和 h²/6，不是 h/6，量纲才对得上。写完之后做两个自检：一是eval_spline在节点处的值必须等于 y（误差 1e-12 以内），二是两端的一阶导数应该逼近 0.8 和 0.2（差分数值近似即可）。这两条过了，S(i+0.5) 的结果基本不会错。

3. 重积分外推：从 Tn 到 Sn 的那一步才是精度来源

二重积分这题看着比样条简单，其实考点全在“逐次二分 + 外推”上。题目给的是 I(f) = ∫_c^d ( ∫_a^b f(x,y) dx ) dy，要求用复化梯形公式、逐次二分步长、外推，算到相邻两次近似值之差的绝对值不超过 ε。如果只是把两次梯形值一减就交，那外推等于没做，精度停在 O(h²)，二分十几次也追不上要求。

3.1 为什么是 4/3 和 -1/3：一次外推把误差阶抬上去

复化梯形公式的误差展开是 T(h) = I + C·h² + O(h⁴)。步长减半后 T(h/2) = I + C·h²/4 + O(h⁴)。两式一消，把 h² 项干掉：

S = (4/3)·T(h/2) − (1/3)·T(h) = I + O(h⁴)

这就是资料里写的 S_n = (4/3)T_{2n} − (1/3)T_n。系数 4/3 和 -1/3 不是拍脑袋来的，是消 h² 项消出来的：4/3 − 1/3 = 1 保证常数项对，4/3·(1/4) − 1/3 = 0 保证 h² 项对。写代码时这两个系数写反、或者写成 3/4 和 1/4，结果会稳稳地偏离真值，而且二分越多次偏离越“稳定”，看起来特别像收敛了。

理解这一点之后，判断依据就很清楚：外推后的 S 序列收敛速度应该明显快于 T 序列。如果打印出来发现 S 和 T 收敛得一样慢，八成是系数或者 T_{2n}/T_n 的顺序搞反了。

3.2 内层 x、外层 y：二重循环的顺序不能乱

题目写的是 ∫_c^d (∫_a^b f(x,y) dx) dy，内层对 x 积、外层对 y 积。程序里两层循环必须严格对应：外层循环变量是 y，从 c 到 d；内层是 x，从 a 到 b。我见过不少把两层换过来的写法，程序照样出数、照样收敛，只是算出来的是另一个积分的值——因为积分区域是矩形的时候，交换次序数学上等价，但如果被积函数里的 x、y 在调用时传反了（比如f(x,y)写成了f(y,x)），结果就悄悄错了。

被积函数 f(x,y) = tan(x² + y²) 里 x、y 是平方和，交换了看不出差别，这类题最容易掩盖参数顺序错误。我的习惯是给临时测试函数f(x,y) = x*10 + y，跑一遍看外层积分变量是不是真的在外面，两分钟就能确认，比事后调半天强。

积分限这一题是 a=0、b=π/3、c=0、d=π/6，上界小、下界大这种写法在题目里也常见，(b-a)一减如果是负的，整个积分符号就反了。写完先把积分限代进hx = (b-a)/nx打一行看看正负。

3.3 三个文件的拆分：myfun / getT / 主程序

资料里把被积函数、复化梯形子函数、主程序分成三个文件，这个拆法是对的——被积函数一改就能换题，梯形子函数可以单独测。被积函数：

function z = myfun(x, y) % 被积函数，注意 x 是内层积分变量、y 是外层 z = tan(x.^2 + y.^2); end

复化梯形二重积分，注意权重：角点乘 1/4，边上乘 1/2，内部乘 1。写代码时把权重拆成wx * wy最不容易错：

function T = getT(f, a, b, c, d, nx, ny) % 矩形区域上的复化梯形二重积分，nx/ny 为 x/y 方向等分数 hx = (b-a)/nx; hy = (d-c)/ny; T = 0; for i = 0:ny y = c + i*hy; wy = 1; if i==0 || i==ny, wy = 0.5; end % y 方向端点权重 inner = 0; for j = 0:nx x = a + j*hx; wx = 1; if j==0 || j==nx, wx = 0.5; end % x 方向端点权重 inner = inner + wx * f(x, y); end T = T + wy * inner; end T = T * hx * hy; % 最后统一乘步长 end

hx*hy放在循环外面统一乘，比每次内层都乘一遍快，也避免累加过程中浮点误差被放大。这个写法在 nx、ny 都到几千的时候，差别能看出来。

3.4 二分停止条件写在哪儿，决定了会不会死循环

题目要求“相邻两次近似值之差的绝对值不超过 ε”，ε 取 0.5e-5。关键是这个“近似值”指的是外推后的 S，不是原始的 T。判据写成abs(T2 - T1) < eps就偏了题意，写成abs(S2 - S1) < eps才对。

eps0 = 0.5e-5; nx = 2; ny = 2; Tn = getT(@myfun, 0, pi/3, 0, pi/6, nx, ny); % 粗网格 nx = 2*nx; ny = 2*ny; T2n = getT(@myfun, 0, pi/3, 0, pi/6, nx, ny); Sn = 4/3*T2n - 1/3*Tn; % 第一次外推 it = 0; while it < 20 % 上限兜底，防死循环 Tn = T2n; nx = 2*nx; ny = 2*ny; T2n = getT(@myfun, 0, pi/3, 0, pi/6, nx, ny); S2n = 4/3*T2n - 1/3*Tn; if abs(S2n - Sn) < eps0 break; end Sn = S2n; it = it + 1; end fprintf('积分值 = %.10f，二分次数 = %d\n', S2n, it);

it < 20这个上限不是可选项，是必需品。浮点数在 1e-5 附近会进入震荡区，两个相邻的 S 值可能来回跳，永远不满足小于 ε，没上限就是死循环。而且二分到 nx、ny 上千以后，单次 getT 的耗时是平方级增长的，卡死一次能等十几分钟。

3.5 运行慢的真因：每次二分都在重算旧节点

资料里反思“运行速度慢是很大的不足”，这个判断是准确的，但原因不是“MATLAB 慢”，而是每次getT都从头把所有节点重算了一遍。实际上复化梯形公式有递推关系：

T_{2n} = T_n / 2 + (h/2) · Σ f(新增节点)

新增节点的横坐标是 a + (2k-1)h/2，数量正好是 n 个。用这个递推，每次二分的计算量从 O(4n²) 降到 O(n²) 的一小部分，尤其两个方向同时二分的时候，节省非常明显。改法是把getT拆成“首算”和“加细”两个函数，主循环里只算新增点。做作业不一定要求，但把这段写进反思里，比写“还需要更深入学习 MATLAB”有说服力得多。

4. 常微分方程初值：RK4、AB4、AB4-AM4 和改进版的误差差在哪

第三道题最有意思，因为它把四种方法放在同一个初值问题上横向比。初值问题是 y' = −x²y²（0 ≤ x ≤ 1.5），y(0) = 3，步长 h = 0.1，精确解 y = 3/(1+x³)。四种方法、同一组数据、同一个误差范数，排序出来一目了然。这道题真正考的不是“会不会写 RK4”，而是“知不知道为什么 AB4 会比 RK4 差两个数量级”。

4.1 先把精确解和右端函数准备好

右端函数和精确解各写一个文件，后面所有方法共用：

function f = fxy(x, y) % 一阶导数右端项 f = -x.^2 .* y.^2; end function y = fx(x) % 精确解，用于误差对比 y = 3 ./ (1 + x.^3); end

精确解这个文件很容易被当成“可有可无”，实际它是整道题的裁判。没有它，四种方法的误差没法排序，题目第（6）问的“能得到哪些结论”就只能凭感觉写。我一般会把精确解和数值解画在同一张图上，再加一张误差随 x 变化的图，误差曲线的形状比单个数字能说明问题——RK4 的误差曲线几乎是平的，AB4 的误差曲线一路往上爬，这个差别一眼就能看出来。

验证一下精确解和右端项一致：对 y = 3/(1+x³) 求导得 −9x²/(1+x³)²，而 −x²y² = −9x²/(1+x³)²，两边相等。这一步花三十秒，能避免后面所有对比都建立在一个写错的 fx 上。

4.2 RK4 是单步自启动，AB4 不是

RK4 是单步法，只依赖 y_i 就能算出 y_{i+1}，所以第一个点起步毫无障碍：

function [x, y] = rk4(f, x0, xend, y0, h) % 经典四阶 Runge-Kutta x = (x0:h:xend)'; N = length(x); y = zeros(N,1); y(1) = y0; for i = 1:N-1 xi = x(i); yi = y(i); k1 = f(xi, yi); k2 = f(xi + h/2, yi + h/2*k1); k3 = f(xi + h/2, yi + h/2*k2); k4 = f(xi + h, yi + h*k3); y(i+1) = yi + h/6 * (k1 + 2*k2 + 2*k3 + k4); end end

四个 k 的权重是 1:2:2:1，除以 6，这是四阶精度的来源。k2、k3 都是在中点取值，区别只在于用 k1 还是 k2 去估斜率，这一点写反了阶数会掉到二阶，误差大约变 h² 倍——h=0.1 时差 100 倍左右，但程序不报错，只能靠和精确解比才看得出来。

AB4 是四阶显式 Adams 方法，用前面四个点的 f 值外推：

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

这里有个极易抄错的地方：系数是 55、−59、37、−9，前面的因子是 h/24，不是 h/12。资料正文里这一行写的是 h/12，标准四阶 AB4 的权系数分母是 24（对 55、−59、37、−9 来说，24 才能让它们加起来等于 24，保证常数函数精确成立：55−59+37−9 = 24）。这是个非常实用的自检：四个系数加起来必须等于分母，不等就说明抄错了。

AB4 自身不能启动，因为 i=1 时需要 f_0、f_{-1}、f_{-2}，前三个点得由别的单步法提供。题目要求用 RK4 提供初值，那就先跑 RK4 得到 y_1、y_2、y_3，从 i=4 开始用 AB4 递推。这里有个坑：RK4 提供的初值必须用同一个步长 h=0.1，不能图省事用更小的步长先算前三个点，否则后面 AB4 用的 f 序列步长不一致，整个递推全废。

4.3 AB4-AM4 的预测校正：先猜后修

AB4 是显式方法，用已知点外推，稳定性差；AM4 是隐式方法，公式右边含 y_{i+1}，得迭代解。把两者接起来就是预测校正：先用 AB4 猜一个 y_{i+1}^{(p)}，代进 AM4 的右端算 f_{i+1}，再修正一次得到 y_{i+1}^{(c)}。

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}]

注意这里只做一次校正，不迭代到收敛。迭代到收敛当然更准，但计算量上去了，而且四阶方法的截断误差摆在那儿，多迭代几次收益有限。只做一次校正是 PECE 模式（Predict-Evaluate-Correct-Evaluate），工程上最常用的折中。

校正这一步把 AB4 的误差压下去大约一个量级，因为隐式公式的稳定域比显式大得多。这也是为什么资料里 AB4 的误差 2 范数是 0.004809，而 AB4-AM4 降到 0.000920，差距五倍多。

4.4 改进版多出来的那一步：251/270 从哪来

改进的 AB4-AM4 在预测值和校正值之间再做一次加权组合：

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

这两个系数不是凑的，是把预测值和校正值的误差主项做抵消推出来的。251 + 19 = 270，权重归一，保证常值解不被破坏。校正值权重远大于预测值，说明这一步是“以校正值为主、拿预测值做小幅修正”。

function [x, y] = ab4am4_improved(f, x0, xend, y0, h) % 改进的 AB4-AM4 预测校正，前三个点由 RK4 提供 x = (x0:h:xend)'; N = length(x); y = zeros(N,1); y(1) = y0; % 用 RK4 起步，步长必须一致 y(1:4) = rk4(f, x0, x0+3*h, y0, h); F = f(x, y); % 预先算好各点 f 值，避免重复调用 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)); y(i+1) = 251/270 * yc + 19/270 * yp; % 加权组合 F(i+1) = f(x(i+1), y(i+1)); % 更新 f 序列，供下一步用 end end

F这个数组是关键优化：AB4 和 AM4 都要回头取前几个点的 f 值，每次都重新调f是浪费。预先存成一个向量，循环里直接索引，代码也短。注意F(i+1)要在算完 y(i+1) 之后更新，顺序反了就会用到旧的 f 值。

y(1:4) = rk4(...)这一行返回的是 4 个点，正好覆盖 i=0 到 i=3，从 i=4 开始 AB4 才有足够的四个历史点可用。这里如果返回 3 个点，后面索引会从 4 开始但 F(1) 是空的，MATLAB 不报错、直接拿 0 参与运算，结果全错。

4.5 四种方法跑同一组数据，误差 2 范数排序

统一步长 h=0.1，x 从 0 到 1.5，共 16 个点。误差用和精确解之差的 2 范数：

xe = (0:0.1:1.5)'; ye = fx(xe); [~, y1] = rk4(@fxy, 0, 1.5, 3, 0.1); [~, y2] = ab4(@fxy, 0, 1.5, 3, 0.1); % 内部用 RK4 起步 [~, y3] = ab4am4(@fxy, 0, 1.5, 3, 0.1); [~, y4] = ab4am4_improved(@fxy, 0, 1.5, 3, 0.1); fprintf('RK4=%.6f AB4=%.6f AB4-AM4=%.6f 改进=%.6f\n', ... norm(y1-ye), norm(y2-ye), norm(y3-ye), norm(y4-ye));

结果对照如下：

精度从高到低：RK4 > 改进 AB4-AM4 > AB4-AM4 > AB4。这个排序和理论完全吻合。RK4 虽然是单步法，但每一步只依赖当前点，误差不累积历史信息，而且四阶精度在 h=0.1 这种不算小的步长下非常稳。AB4 是显式多步法，用到前四个点的 f 值外推，误差主项系数大（55/24 量级），加上稳定性区域窄，误差一路累积。AB4-AM4 靠隐式校正把误差压下来五倍，改进版再靠加权组合压到 RK4 的十四倍以内。

值得注意的一点：AB4-AM4 比 AB4 准，但计算量是两倍（每次迭代多算一次 f）。如果只看“精度/计算量”，RK4 在这道题上是完胜的——这也是为什么工程上 RK4 用得远比 Adams 系列广。Adams 系列的优势在于右端函数 f 特别贵、且步长可以取得很小时，多步法每步只算一次新的 f，总调用次数少于 RK4。

5. 上机避坑：五个反复出现的现象、原因和改法

下面这五条，都是上机过程中真实出现过、而且在新一批作业里还会再出现的。每一条都按“看到什么现象 → 为什么会这样 → 怎么改”写，照着核对能省不少时间。

5.1 公式与数据层面的坑

坑一：样条曲线整体平滑，但节点处一阶导明显断裂。

现象：把 S(x) 在 [0,10] 上密集取点画图，曲线看着连续，但放大到节点附近能看出“折角”，和理论上的 C² 连续对不上。

原因：三弯矩方程里 μ 和 λ 的挂靠位置写反了。μ_j 乘的是 M_{j-1}，λ_j 乘的是 M_{j+1}，两者数值上互补（和为 1），写成同一行时肉眼看不出区别。方程组照样有唯一解，解出来的 M 序列整体错位。

解决：等距情况下 μ=λ=0.5，这个特例下写反了程序也能跑对，所以不能只拿等距数据测。造一组不等距的 x（比如 [0 1 3 6 10]），手算第一行的 μ、λ、d，和程序打印的中间变量逐项对。另外在回代后加一条检查：max(abs(diff(eval_spline(x,y,M,x))))应该逼近两端给定的导数值，不逼近就是方程组错了。

坑二：三次样条第一个区间结果离谱，后面全部正常。

现象：S(0.5) 的值大得离谱，S(1.5) 之后基本合理。

原因：型值点录入时2.51和3.30被粘成了13.30，第二个点跳高。三弯矩方程的第一行右端项直接由这个点算出来，误差被放大进 M_0 和 M_1。

解决：数据录入后先算一阶差分，检查有没有异常跳变。等距情况下差分应该缓慢变化，突然出现一个 +9 的跳变百分之百是录入错误。另外把 y 的均值、极差打出来看一眼，车门曲线这种物理量的极差应该在合理范围内。

坑三：AB4 的误差大得离谱，比理论值差好几个量级。

现象：RK4 误差 3.8e-5，AB4 误差 4.8e-3，差了两个数量级，比预期的大。

原因：AB4 公式那一行的因子抄成了 h/12。分母 24 是四阶 AB4 的标配，分母写成 12 相当于把整步长放大了一倍。还有一个隐藏原因：前三个点如果不用同一步长的 RK4 提供，f 序列的步长不一致，递推公式不成立。

解决：用系数和检查——55 − 59 + 37 − 9 = 24，等于分母。这一条对任何 Adams 类公式都适用，写之前先加一遍。前三个点严格用同一个 RK4 函数、同一个 h 算，不做任何“先小步长跑准一点”的优化。

坑四：二重积分结果比预期大一倍或小一半。

现象：算出来的积分值总是和参考答案差个整数倍，二分次数多了也收敛不到正确答案。

原因：复化梯形的权重分配错了。二重积分的四个角点权重是 1/4，四条边上的点是 1/2，内部点是 1。常见错误是只在外层考虑权重、内层全用 1，或者反过来，结果正好差一个因子。

解决：拿 f(x,y)=1 做单元测试，积分值应该严格等于 (b-a)(d-c)。这个测试五秒钟，能立刻暴露权重错误。然后再拿 f(x,y)=x+y 测，积分值应该等于 (b-a)(d-c)(a+b+c+d)/2，能把内层外层的变量顺序也一并验出来。

5.2 程序结构与调用层面的坑

坑五：函数文件明明存在，调用却报“未定义函数”。

现象：myfun.m就在当前目录，主程序里myfun(0,0)报错找不到。

原因：MATLAB 要求文件名和函数名严格一致，且函数文件不能放在脚本内部（除非用局部函数语法）。另一种常见情况是主程序本身是个脚本，函数定义写在了脚本中间，旧版本 MATLAB 会报错。

解决：文件名、函数名、调用名三者统一，大小写也要一致。检查当前工作目录是不是文件所在目录，用which myfun能直接看到 MATLAB 实际找到的是哪个文件——这个命令我每次都先用一遍，比猜快得多。

坑六：改进的 AB4-AM4 结果和 AB4-AM4 完全一样。

现象：两个方法输出的误差一模一样，改进似乎没生效。

原因：加权组合那一行写成了y(i+1) = yc，或者系数写成了251/270 + 19/270 = 1但把 yp 和 yc 的位置写反了。位置写反的情况下结果会明显变差，改成“几乎一样”通常是系数只写了一半。

解决：在加权那一行后面加一句if i==5, fprintf('yp=%.6f yc=%.6f y=%.6f\n', yp, yc, y(i+1)); end，把第一次迭代的三个数打出来，手算一遍 251/270yc + 19/270yp 对不对。两三行调试代码能定位的问题，不要靠盯着公式找。

6. 误差排序之后：按精度选方法，按习惯保结果

三道士机题走完，最有价值的结论来自第三题那张误差表：RK4 0.000038、改进 AB4-AM4 0.000531、AB4-AM4 0.000920、AB4 0.004809。这个排序不是偶然，它把“单步 vs 多步”“显式 vs 隐式”“是否校正”三组设计选择的代价和收益摆在了明面上。单步法误差不累积历史、显式法稳定性差、隐式校正能压误差但计算量翻倍——这三条判断依据，比记住任何一个公式都管用。

真正要带走的一个技巧是：用水滴法验证阶数。把步长从 h 减到 h/2，看误差是不是按 2⁴ = 16 倍下降。RK4 从 h=0.1 换到 h=0.05，误差应该降到 1/16 左右；如果只降到 1/4，说明实际是二阶方法，代码里某个地方阶数掉了。这个检查对任何数值方法都适用，不依赖参考答案：

% 阶数验证：误差比应接近 2^p，p 为方法阶数 for h = [0.1 0.05 0.025] xe = (0:h:1.5)'; [~, yn] = rk4(@fxy, 0, 1.5, 3, h); err = norm(yn - fx(xe)) / sqrt(length(xe)); % 均方根误差，消除点数影响 fprintf('h=%.3f rmse=%.3e\n', h, err); end

注意这里用的是均方根误差而不是原始 2 范数，因为不同步长下点数不一样，直接比 2 范数会混进点数的影响。连着三个 h 打印出来，误差比稳定在 16 附近，才敢说这段 RK4 写对了。

回过头看，我在这三道题上栽的跟头，一半是公式抄写，一半是数据结构没想清楚就开始敲。后来养成的习惯是：任何数值程序，写之前先拿一个能手工验算的小规模输入（等距、两三个点、简单函数），把中间变量打印出来逐项对，对上了再换真题数据。这条习惯救过我不止一次——三次样条那次，就是因为等距特例下 μ 和 λ 都是 0.5，写反了也看不出来，才拖到别人帮忙跑的时候才发现。从那以后我每次写三对角方程组，都强制先用不等距的 5 个点跑一遍手算对照，再也不敢拿等距数据当验收标准。希望帮到你。

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