从三次样条到RK4：局部信息组合出高阶全局精度

简介：这是一份面向学习数值分析课程、需要完成上机作业的高校学生与自学者的资料；它基于东南大学出版社教材，借助 MATLAB 编程实现三类典型数值计算问题：编制第一型三次样条插值函数的通用程序，并处理车门曲线型值点算例；用复化梯形公式配合逐次二分和外推思想计算二重积分；用 RK4、AB4、AB4-AM4 预测校正等方法求常微分方程初值问题，并将结果与精确解对比误差。压缩包共 1 个 docx 文件，约 566KB，按上机题四、五、六分节记录题意、算法推导、程序结果与感悟反思，便于对照教材章节顺序查阅。目前已有 905 人学习。读者可据此获得能直接改写复用的 MATLAB 子函数框架、边界条件与步长的处理细节，以及四种方法精度排序等结论性对照，适合作为课程作业参考与考前复盘材料。

1. 三道题共用同一条主线：用局部信息组合出全局精度

数值分析上机题最容易遇到的问题，是公式抄对了、程序也跑通了，可打印出来的 S(i+0.5) 跟同学差一个数量级。三道题——三次样条插值、重积分计算、常微分方程初值问题数值解——看起来是三个互不相干的算法，实际上共用同一条底层逻辑：用低阶的局部信息加权组合，组合出高阶的全局精度。三次样条把 4n 个多项式系数压缩成 n+1 个弯矩，靠的是二阶导连续；复化梯形用 4/3 和 −1/3 一配，h² 误差项当场消失；AB4 预测、AM4 校正、再按 251/270 与 19/270 混合一次，还是同一个思路。看懂这条线，MATLAB 里那些 55、−59、37、−9 就不再是需要硬背的常数。这份资源把三题的通用程序、边界条件处理、误差对比都整理成了可直接运行的代码，下面从样条的未知数计数开始拆。

2. 三次样条插值：为什么“分段三次 + 二阶连续”最终会形成一个三对角方程组

2.1 先数未知数：4n 个系数，只有 4n−2 个条件

把区间 [x₀, xₙ] 切成 n 段，每段一个三次多项式，每段 4 个待定系数，一共 4n 个未知数。约束条件呢？每段两端函数值必须等于型值点，2n 个；内节点处左右两段的一阶导相等，n−1 个；二阶导也相等，n−1 个。加起来 4n−2，比未知数少 2 个。

这 2 个缺口就是所谓的边界条件，三次样条必须补上两个端点约束，否则解不唯一。常见的三类：自然边界（两端二阶导为 0，M₀ = Mₙ = 0）、第一型边界（给定端点一阶导）、周期边界（首尾导数与二阶导都相等）。车门曲线那道题给的是 y₀′ = 0.8、y₁₀′ = 0.2，属于第一型——端点斜率是已知的，不允许你取自然边界那种简单处理。

提示：少补两个条件不是“程序能跑就行”的小事。强行令 M₀ = Mₙ = 0 跑出来的曲线在两端会明显翘起来，越靠近端点误差越大，因为自然边界等价于假设端点处没有弯矩，跟真实的车门轮廓对不上。

2.2 换未知数：把 4n 个系数压缩成 n+1 个弯矩

直接解 4n 个系数既不方便编写代码又不方便调试。弯矩法换了个角度：三次多项式的二阶导是一次函数，所以整个 S″(x) 是分段线性的。记 Mⱼ = S″(xⱼ)，那么在第 j 段上

S″(x) = Mⱼ·(x_{j+1} − x)/hⱼ + M_{j+1}·(x − xⱼ)/hⱼ

对它积分两次，并用 S(xⱼ) = yⱼ、S(x_{j+1}) = y_{j+1} 这两个端点值把积分常数定下来，整段 S(x) 就被 Mⱼ、M_{j+1} 完全确定了。也就是说，插值条件自动满足，剩下的只有一个约束没用：内节点处一阶导连续。把它写出来，就得到经典的弯矩方程

μⱼ·M_{j−1} + 2·Mⱼ + λⱼ·M_{j+1} = dⱼ，j = 1, 2, …, n−1

其中 μⱼ = h_{j−1}/(h_{j−1}+hⱼ)，λⱼ = hⱼ/(h_{j−1}+hⱼ) = 1 − μⱼ，右端 dⱼ = 6·f[x_{j−1}, xⱼ, x_{j+1}]，也就是三点二阶差商乘 6。二阶差商对二次函数正好等于 f″/2，所以 6f[·,·,·] 在量纲上跟 2M 对得上——这个细节能帮你记住右端到底该乘 6 还是乘 3。

每个方程只牵涉 M_{j−1}、Mⱼ、M_{j+1} 三个未知量，整张系数矩阵是三对角的。这就是为什么 n 再大也不怕：三对角系统的解法是 O(n)，而稠密消元是 O(n³)。

2.3 第一型边界：端点导数怎么变成第一行和最后一行

内点方程有了 n−1 个，还差两个，正好由端点条件提供。第一型边界把 M₀ 和 M₁ 的关系写出来：

2·M₀ + M₁ = (6/h₀)·( (y₁ − y₀)/h₀ − y₀′ )

右端那项是“一阶差商减去端点导数”，物理含义很清楚：如果第一段用直线去拟，端点的斜率就应该是差商；真实斜率比它大还是小，差值决定了端点弯矩的大小。右端符号写反是这类作业里最高频的错误之一，结果就是整条曲线在左端系统性偏移，但中点值看起来“还行”，很容易漏过去。右端点同理：

M_{n−1} + 2·Mₙ = (6/h_{n−1})·( yₙ′ − (yₙ − y_{n−1})/h_{n−1} )

注意右端两项的顺序是“给定的导数减去差商”，和左端是相反的——因为左端差商是向前指的，右端是向后指的。

2.4 解三对角方程组与任意点求值

把三部分拼起来就是完整的程序。下面这段是核心，MATLAB 下标从 1 开始，而公式从 0 开始，代码里的 j 对应公式里的 j−1，这个平移是主要难点：

function [M, h] = spline_moment_1st(x, y, dyL, dyR) % 第一型三次样条：给定左右端点一阶导，返回各节点弯矩 M % x, y 长度均为 n+1；dyL/dyR 为 y'(x0)、y'(xn) x = x(:); y = y(:); % 强制列向量，防止行/列混用让 diff 方向错掉 n = numel(x) - 1; % 区间个数 h = diff(x); % h(j) = x(j+1) - x(j)，j = 1..n A = zeros(n+1); d = zeros(n+1,1); A(1,1) = 2; A(1,2) = 1; % 左端第一型边界 d(1) = 6/h(1) * ((y(2)-y(1))/h(1) - dyL); for j = 2:n % 内点，对应公式 j = 1..n-1 mu = h(j-1) / (h(j-1) + h(j)); lam = h(j) / (h(j-1) + h(j)); 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 A(n+1, n) = 1; A(n+1, n+1) = 2; % 右端第一型边界 d(n+1) = 6/h(n) * (dyR - (y(n+1)-y(n))/h(n)); M = A \ d; % 规模小时直接解；几千点以上换追赶法 end

参数说明：x、y是型值点，车门那道题是 11 个点、x = 0…10、y 从 2.51 到 5.80；dyL、dyR传 0.8 和 0.2。M返回长度为 n+1 的弯矩向量，一阶、二阶导的信息全在里面。用A \ d而不是自己手写追赶法，对 11 个点的作业完全够用，但要知道它内部走的是带选主元的 LU，复杂度不是最优的；节点上千时手写 Thomas 算法能省一个数量级的常数。

拿到 M 之后还要能算 S(i+0.5)。直接用弯矩形式的分段表达式最稳妥：

function S = spline_eval(x, y, M, xs) % 用弯矩向量 M 在任意点 xs（向量）处求插值函数值 h = diff(x); S = zeros(size(xs)); for k = 1:numel(xs) i = min(find(xs(k) <= x, 1) - 1, numel(h)); % 定位落在第 i 段 i = max(i, 1); % 防止 xs 等于 x(1) 时 i 变成 0 t1 = (x(i+1) - xs(k)) / h(i); t2 = (xs(k) - x(i)) / h(i); S(k) = M(i)*h(i)^2/6 * (t1^3 - t1) ... + M(i+1)*h(i)^2/6 * (t2^3 - t2) ... + y(i)*t1 + y(i+1)*t2; end end

这个式子的结构值得多看一眼：前两项是弯矩贡献，(t³ − t)在端点处恰好为 0，所以不会破坏插值条件；后两项是一次线性插值，负责把两端的函数值接上。两种效应叠加，才既过型值点又保持二阶导连续。求 i+0.5 时，直接把xs = (0:9) + 0.5传进去，一次拿到 10 个值。

3. 重积分：复化梯形逐次二分再加外推，把 h² 误差项消掉

3.1 二重积分先分解为两层一维积分

I(f) = ∫_c^d (∫_a^b f(x,y) dx) dy 这种形式，处理时不需要发明新的二重积分公式。固定一个 y，内层就是一个普通的一维积分，用复化梯形算出结果；外层再对这些“关于 y 的函数值”做一次一维复化梯形。总误差是两层误差的和，量级 O(h²) + O(k²)，h 是 x 方向步长，k 是 y 方向步长。

这也解释了一个常见困惑：为什么题目要求“逐次二分步长”。因为误差主项是 h² 和 k² 的线性组合，两个方向必须同步加密，否则你只把 h 减半、k 不动，总误差被 k² 卡住，怎么迭代都降不下去；程序看起来在收敛，其实早就进入平台期了。

3.2 4/3 与 −1/3 的来历：误差展开后的一次消元

复化梯形公式的误差有渐近展开：

T(h) = I + C₁h² + C₂h⁴ + …

把步长减半再算一次：T(h/2) = I + C₁h²/4 + C₂h⁴/16 + …

两个式子一联立，把 C₁h² 消掉：4/3·T(h/2) − 1/3·T(h) = I + O(h⁴)。这就是 Richardson 外推，也正好等价于复化 Simpson 公式。所以题目里那句“按公式 (1) 做线性组合得到 Sₙ(f)”，本质上就是用两次梯形值换一次四阶精度——多算的点没白算，信息被复用了一次。

这一点是整道题的核心：单纯把步长减半，误差降 4 倍；做一次外推，误差降 16 倍。代价只是多一行减法。

3.3 二分时只算新点：复化梯形的递推

朴素写法是每次二分都从头把所有点重算一遍，函数求值次数直接翻倍再翻倍。实际上二分后的新节点正好是原区间的中点，复化梯形有个递推关系：

T(2n) = T(n)/2 + (h/2)·Σ f(新中点)

原来的点全部复用，只算新加的那一半点。函数 f 求值次数从 2n 降到 n，这是这类程序里最容易忽略、收益却最大的优化。

function [S, T, n] = trapz_richardson(f, a, b, tol, maxit) % 复化梯形 + 逐次二分 + 一次 Richardson 外推（结果等价于复化 Simpson） h = b - a; n = 1; T = h/2 * (f(a) + f(b)); % T_1，两个端点 Sold = inf; for k = 1:maxit xodd = a + h/2 : h : b - h/2; % 只取新增的中点 T2 = T/2 + h/2 * sum(f(xodd)); % 递推得到 T_{2n} S = 4/3*T2 - 1/3*T; % 消掉 h^2 项，进入四阶 if abs(S - Sold) <= tol % 停机：相邻两次外推值之差 T = T2; n = 2*n; return end Sold = S; T = T2; n = 2*n; h = h/2; end warning('达到最大迭代次数，tol 可能设得比舍入误差还小'); end

参数说明：f是被积函数句柄，a/b是积分区间，tol是精度要求，maxit是迭代上限（必须有，否则浮点舍入误差会让差值永远降不到 tol 以下，循环停不下来）。T返回最后一次的梯形值，S返回外推值，n是最终的等分份数。

二重积分就是把这个一维过程套两层：外层每个 y 节点上调用一次内层积分，外层自身也走二分 + 外推。题目给的 ε = 1.2×10⁻⁵ 是个不算苛刻的精度，实际跑下来外推三五次就够了。

3.4 ε 设定的隐含前提：被积函数在区间内要足够光滑

题目的被积函数是 tan(x² + y²)，x ∈ [0, π/3]，y ∈ [0, π/6]。x² + y² 在这个矩形上的最大值约 1.37 rad，离 π/2 还有距离，tan 在区间内光滑、没有奇点，所以复化梯形的误差展开成立，外推才有意义。

如果换成 x² + y² 能接近 π/2 的函数，情况就完全变了：被积函数在积分区间内出现陡峭上升，高阶导数急剧增大，渐近展开的主项被高次项淹没，外推会表现得像随机数——有时收敛，有时越算越差。这时该做的是换积分方法（比如自适应分割或者变量替换消奇点），而不是把 ε 设置得更小。理解“外推依赖光滑性”这个前提，比记住公式的系数重要得多。

4. 常微分方程初值问题：RK4 起步、AB4 外推、AM4 校正的精度阶梯

4.1 单步与多步的分工：为什么必须先用 RK4 启动

初值问题 y′ = −x²y²，y(0) = 3，精确解 y = 3/(1 + x³)。题目要求四种方法都用同一个步长 h = 0.1，并且互相比较。

RK4 是单步法：算 y_{i+1} 只需要 y_i，代价是每步要 4 次 f 求值，但精度是四阶。AB4 多步法是显式公式，每步只需要 1 次 f 求值，但算 y_{i+1} 要往前用四个历史点 y_i、y_{i−1}、y_{i−2}、y_{i−3}。问题来了：i = 0 的时候没有历史点。所以 AB4 必须由 RK4 提供前四步的初值，从第 4 步开始才切到多步公式。这个“起步段”的精度直接决定了后面整条轨迹的精度上限，不是可以随便用欧拉法凑合的地方。

4.2 AB4 的 55/−59/37/−9 是插值积分出来的

AB4 的形式是

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

系数看着奇怪，来源其实很直接：把 f 在 x_i、x_{i−1}、x_{i−2}、x_{i−3} 四个点上做三次拉格朗日插值，然后在 [x_i, x_{i+1}] 上对这个插值多项式积分，得到的积分权重乘以 h 就是这组系数。分母 24 来自 4 个点插值积分后通分的结果。

这也顺带解释了它为什么叫“显式”：用的是历史点上的 f 值，y_{i+1} 直接算出来，不需要解方程。代价是精度上限被插值外推限制住——把三次多项式外推到当前区间之外去积分，误差常数 251/720·h⁵ 比 RK4 大得多。所以理论上同为四阶，AB4 的实际误差往往差一个数量级，后面那张表会证实这一点。

4.3 AM4 是隐式的：右端含有待求的 y_{i+1}

AM4 的公式是

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

注意第一项用的是 x_{i+1} 处的 f，而那个点的 y 正是你要求的未知量。这就是隐式法的本质：等式两边都有 y_{i+1}，不能直接算，得用迭代求解或做预测校正。

预测校正的思路是：先用 AB4 把 y_{i+1} 算出来当预测值 y^{(p)}，代入右端算出 f(x_{i+1}, y^{(p)})，再代回 AM4 得到校正值 y^{(c)}。整个过程每步 2 次 f 求值，比 RK4 少一半，精度却比纯 AB4 高得多——因为 AM4 的误差常数只有 −19/720·h⁵，比 AB4 小一个数量级。

4.4 再混合一次：251/270 与 19/270 的作用

“带改进的预测校正”在预测值和校正值之上又做了一次加权：

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

两个系数之和正好是 1，所以这是一个凸组合，结果一定落在预测值和校正值之间，不会发散。它的作用是把残留在校正值里的低阶误差项再消一次——预测值和校正值的误差主项符号相反、大小成比例，按这个比例混合能让主项相互抵消，等效阶数进一步提高。

function [x, y] = ode_ab4_am4(f, x0, x1, y0, h, mode) % mode: 'ab4' | 'am4pc' | 'am4imp'（改进的预测校正） x = (x0:h:x1)'; N = numel(x); y = zeros(N,1); y(1) = y0; % --- 起步段：RK4 提供前四个点 --- for i = 1:3 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 % --- 多步段：i 从 4 开始，历史点用 i、i-1、i-2、i-3 --- for i = 4:N-1 f0 = f(x(i), y(i)); f1 = f(x(i-1), y(i-1)); f2 = f(x(i-2), y(i-2)); f3 = f(x(i-3), y(i-3)); yp = y(i) + h/24*(55*f0 - 59*f1 + 37*f2 - 9*f3); % AB4 预测 if strcmp(mode, 'ab4') y(i+1) = yp; continue end fp = f(x(i+1), yp); yc = y(i) + h/24*(9*fp + 19*f0 - 5*f1 + f2); % AM4 校正 if strcmp(mode, 'am4pc') y(i+1) = yc; else y(i+1) = 251/270*yc + 19/270*yp; % 改进：再加权一次 end end end

参数说明：f是右端函数句柄，本题是@(x,y) -x.^2.*y.^2；x0/x1是区间端点 0 和 1.5；y0是初值 3；h取 0.1；mode切换三种多步方案。起步段固定跑 3 次 RK4 得到 y(1)…y(4)，循环里 i 从 4 开始，历史点索引 i、i−1、i−2、i−3 全部合法——这是最容易被下标搞错的地方，写之前先在纸上标一遍。

4.5 误差表怎么读：同为四阶，差了两个数量级

用精确解 y = 3/(1+x³) 逐点对比，取 2-范数：

排序是 RK4 > 改进 AB4-AM4 > AB4-AM4 > AB4，跟理论完全对得上：RK4 用了 4 次求值换来最小的误差常数；AB4 只用 1 次求值，误差常数是 RK4 的几十倍；AM4 校正把误差砍到 AB4 的五分之一；再加权一次又砍掉四成。但要注意这张表是在 h = 0.1、这个特定的方程上测出来的，换个刚性方程或者换个步长，比例关系会变。把“某次实验的排序”当成“方法的固有排名”，是这类作业里最容易带出去的误解。

5. 排查：程序正常跑完但答案不对，先查这几处

5.1 下标基准错位与边界行写反

现象：样条曲线整体能过型值点，但端点附近明显翘起，中点值跟别人的差在小数点后第二位。

原因：公式用 0 基写的，MATLAB 用 1 基，把 M_{j−1} 直接写成M(j-1)时少平移了一位；或者第一型边界右端写成了dyL + 差商，符号反了。后者最隐蔽，因为曲线看上去还是光滑的，只是整体偏移。

解决：解出 M 之后不要立刻往下算，先验证两件事——把A*M - d打出来，范数应该是 1e-14 量级；再用 M 回代检查端点二阶导是否满足边界方程。两行disp能省掉半小时排查。

5.2 停机准则选错、迭代上限缺失

现象：重积分程序跑了几十次迭代还在继续，或者报出“相邻两次之差反而变大”。

原因有两种。一是停机准则用了 |T(2n) − T(n)| 而不是外推后的 |S_new − S_old|，前者是二阶量、后者是四阶量，用错了会提前停机或者白跑。二是 tol 设到了机器精度以下（比如 1e-18），浮点舍入误差形成的噪声底永远达不到，差值在一个小范围内来回抖。

解决：必须带maxit，并且在外推值不再单调下降时主动退出。判断方法很直接——把abs(S - Sold)打印出来，如果它在某个量级上下振荡而不是持续下降，就说明已经已经达到舍入误差的底限了，加密没有意义。

5.3 起步段和历史点索引错位

现象：RK4 单独跑没问题，一切到 AB4 之后前几步误差突然放大，后面又平稳下来。

原因：多步法的起步段没对齐。要么 RK4 只跑了 3 步却让多步循环从 i = 4 开始（历史点 i−3 = 1 是合法的，但 i = 4 处的 f 还没准备好），要么循环里的 f1、f2、f3 对应错了节点。这种错误不会报错，只会让前几个点偏离，然后误差沿着递推链一路传下去。

解决：单独打印前 5 步的数值解与精确解之差。正常情况起步段误差应该跟 RK4 同量级（1e-6 左右），如果第一步就跳到 1e-3，那就是索引或者公式系数抄错了，别再往后跑。

5.4 二分时重复计算旧节点

现象：结果完全正确，但运行时间随迭代次数指数增长，跑一次几十秒。

原因：每次二分都把整个区间从头求值一遍，函数求值次数变成 2 + 4 + 8 + …，而实际上新增的只有中点。被积函数如果是 tan 这种超越函数，单次求值不便宜，浪费会非常明显。

解决：改用递推式 T(2n) = T(n)/2 + (h/2)·Σf(新中点)，只对新节点求值。改完之后迭代次数不变，耗时通常掉到原来的三分之一以下。

5.5 把“某次实验的误差排序”当成通用结论

现象：报告里写“RK4 精度最高，所以以后一律用 RK4”。

原因：只在一道题、一个步长上做了对比。AB4 每步只要 1 次 f 求值，如果 f 的计算极其昂贵（比如每步要解一个子问题），在相同计算预算下多步法反而可能赢——精度不是唯一指标。

解决：对比时同时记录“达到给定精度所需的 f 求值总次数”，而不只是最终误差。真正有意义的比较是等计算成本下的精度，或者等精度下的成本。

6. 精度与代价的取舍：最高阶的方法为什么不是默认选项

四种方法放在一起看，会得到一个不太符合直觉的结论：精度最高的方法，往往不是应该默认选用的那个。

RK4 误差 3.8×10⁻⁵，代价是每步 4 次 f 求值；改进的 AB4-AM4 误差 5.31×10⁻⁴，代价是每步 2 次。如果 f 只是一行代数表达式，RK4 的额外开销无所谓，直接用最省心的。但如果 f 每次求值要解一个隐式子问题、查一次大表、或者跑一次仿真，成本差异就会被放大十几倍——这时候拿精度换成本是划算的，把步长减半把误差补回来，总成本可能还更低。

验证也不能只看一张误差表。更严格的做法是反算收敛阶：把步长减半，看误差缩小了几倍，四阶方法应该缩小约 16 倍。

% 用已知精确解反算收敛阶：p = log2(e(h)/e(h/2)) for h = [0.2 0.1 0.05 0.025] [~, y1] = ode_ab4_am4(f, 0, 1.5, 3, h, 'am4imp'); e = norm(y1 - exact(0:h:1.5)', 2); % 与精确解的 2-范数误差 fprintf('h=%.3f err=%.3e\n', h, e); end % 相邻两行误差之比应接近 2^4 = 16；只接近 2^2 = 4 说明实际只有二阶

如果算出来比值接近 4 而不是 16，基本可以断定是起步段精度不够，或者某处公式被降阶了（比如把 AM4 的 9 写成了 1）。这个自检比对着答案抠数字有用得多，因为它检验的是方法的阶，而不是某一次的具体数值。

最后一个容易被忽略的取舍是外推的边界。Richardson 外推每次能提高两阶，但每外推一次，参与运算的两个量就越接近，相减时的有效位数损失就越多。迭代到七八层之后，舍入误差会盖过被消掉的截断误差，结果反而更差。我一般在程序里固定一个迭代上限（比如 6 层），到顶就停，宁可接受当前精度，也不让浮点噪声进来。

从那以后我每次写完这类数值程序，都会先拿一个能写出精确解的算例反算收敛阶，对不上就不往下算。这个习惯帮我拦下过好几次“结果看着差不多、其实公式少了个系数”的问题。希望帮到你。

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