MATLAB数值分析实战：样条插值、重积分外推与ODE求解

简介：一份围绕数值分析课程上机作业整理的 docx 文档，面向正在学习数值分析、需要完成 MATLAB 编程实验的本科生与研究生。内容以三次样条插值函数、重积分计算和常微分方程初值问题数值解三条主线展开：既给出第一型三次样条通用程序的构造思路与车门曲线算例，也梳理复化梯形公式配合逐次二分与外推求重积分的实现方式，并对比 RK4、AB4、AB4-AM4 预测校正及改进格式在初值问题中的计算精度。压缩包共 1 个 docx 文件，约 566KB，篇幅紧凑，适合作为课程作业对照与算法复习材料。已有 905 人学习，读者可从中获得公式推导、程序组织方式、结果验证与误差排序等参考，便于理解东南大学版教材相关上机题的求解流程，并据此搭建自己的 MATLAB 实验脚本。

1. 从一份上机作业说起：三道题串起数值分析最容易翻车的三条链路

数值分析的上机题，文档看着像三份互不相关的作业，实际是三条会在真实计算里反复出现的链路：插值、积分、微分方程初值问题。这份资料用 MATLAB 把三道题做完了——三次样条插值求车门曲线、复化梯形加外推算二重积分、RK4 与 AB4-AM4 解常微分方程——每个题都给公式推导、程序结构和运行结果截图，末尾还附了误差对比和方法排序。适合两类人：正在啃数值分析课程、需要照着把作业跑出正确数字的同学；以及手上有一堆 MATLAB 老脚本、想看清楚「主程序 + 子函数」这套写法怎么组织的人。它最大的价值不是答案本身，而是把公式怎么落到索引、边界条件怎么塞进三对角方程组、精度不够时该怀疑哪一行代码，这些没人写在教材里的事摆了出来。

2. 运行链路拆解：exp5.m 主脚本、getT.m 子函数与重积分的外推循环

2.1 目录与文件职责划分

这份资料的目录组织没有统一工程模板，是典型的课程作业式布局：按章分文件夹，每道题一个主脚本，通用逻辑抽成子函数。我按原文提到的 M 文件梳理了一下，合理的结构大致是这样：

这种拆法的好处在改题的时候特别明显：换一个被积函数只动 myfun.m，换一个右端函数只动 fxy.m，求解框架一行不用碰。原文里写「对于不同的定积分的形式，也写作单独的子函数，保证其修改的灵活性」，说的就是这个意思。反过来，如果所有东西都塞进一个脚本，第二次做同类题基本只能复制粘贴再改，改错了还找不出来。

提示：MATLAB 要求函数文件名与函数名一致，myfun.m 里第一行必须是 function z = myfun(x,y)，否则调用会报「未定义函数或变量」。这一点新手最容易忽略，尤其是从别人那里拷来文件后随手重命名。

2.2 主脚本怎么调子函数：以重积分为例

重积分这道题的要求是：给定 I(f)=∫_c^d (∫_a^b f(x,y)dx)dy，取初始步长 h、k 和精度 ε，用复化梯形公式做逐次二分，并应用外推思想，算到相邻两次近似值之差不超过 ε。核心思路是先用步长 (h,k) 算一个 T_n，再把两个方向步长都减半算 T_2n，用线性组合得到 Simpson 级别的近似：

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

主脚本负责的是循环和判据，数值计算全部下沉到 getT.m：

% exp5.m 重积分主程序：二分 + 外推 a = 0; b = pi/3; % 内层 x 的积分区间 c = 0; d = pi/6; % 外层 y 的积分区间 eps = 0.5e-5; % 题面精度要求 maxit = 20; % 二分次数上限，防止死循环 h = (b-a)/2; k = (d-c)/2; % 初始步长，先粗分两段 T1 = getT(a,b,c,d,h,k); % 粗网格的复化梯形值 Sold = T1; % 保存上一次外推结果 for it = 1:maxit h = h/2; k = k/2; % 两个方向同时二分 T2 = getT(a,b,c,d,h,k); % 细网格的复化梯形值 Snew = 4/3*T2 - 1/3*T1; % 外推组合得到 Simpson 级精度 fprintf('iter=%2d T2n=%.10f S=%.10f diff=%.2e\n', ... it, T2, Snew, abs(Snew-Sold)); if abs(Snew-Sold) < eps % 相邻两次近似值之差满足精度 break; end Sold = Snew; T1 = T2; % 滚动前进，注意 T1 要更新成 T2 end fprintf('最终结果 I = %.10f，二分次数 %d\n', Snew, it);

这段代码里有两个容易写错的地方。一是 T1 的滚动更新：外推需要的是「二分前后两个梯形值」，所以每轮结束必须把 T2 赋给 T1，漏了这一句，第二轮开始外推用的还是最初那个粗网格值，结果会稳定地偏。二是 maxit 上限，题面只说「算到满足精度为止」，但没说不满足怎么办，加上限之后至少能看到每轮的 diff 是在收敛还是在震荡。

getT.m 是双层复化梯形，外层对 y、内层对 x，各自端点权重减半：

function T = getT(a,b,c,d,h,k) % 双层复化梯形：外层 y 步长 k，内层 x 步长 h ny = round((d-c)/k); % 外层分段数 T = 0; for iy = 0:ny y = c + iy*k; if iy == 0 || iy == ny wy = 0.5; % 复化梯形端点权重 1/2 else wy = 1; end T = T + wy * innerTrap(a,b,h,y); end T = T * k; % 外层乘步长 k end function s = innerTrap(a,b,h,y) % 固定 y 之后对内层 x 做复化梯形 nx = round((b-a)/h); s = 0; for ix = 0:nx x = a + ix*h; if ix == 0 || ix == nx wx = 0.5; % 内层端点同样减半 else wx = 1; end s = s + wx * myfun(x,y); end s = s * h; % 内层乘步长 h end

参数含义直白：a、b 是内层积分上下限，c、d 是外层，h、k 是两个方向的当前步长。分段数用 round((b-a)/h) 而不是直接 (b-a)/h，是为了避免浮点除法带来的 3.9999999 被当成 3 段，这个细节在步长反复二分之后特别重要，血泪经验。

2.3 被积函数独立成文件与运行方式

myfun.m 只有一行有效代码：

function z = myfun(x,y) % 被积函数：tan(x^2 + y^2) z = tan(x.^2 + y.^2); end

写成独立函数文件，而不是在主脚本末尾用匿名函数 @(x,y) tan(x.^2+y.^2)，主要是为了调试时能单独测：在命令行敲 myfun(0.1,0.2)，看返回值是否合理，比在整个二重循环里打断点快得多。

跑起来的步骤不复杂，但顺序要对：

1. 打开 MATLAB，把当前工作目录切到 ch5_integral 文件夹，或者在「设置路径」里把该文件夹加进去；
2. 命令行输入 exp5，回车；
3. 观察每轮打印的 diff 是否单调下降，最后一行会给出最终积分值和二分次数；
4. 想改精度就改 exp5.m 里的 eps，想换被积函数就改 myfun.m，主程序不用动。

原文档里的运行结果截图显示积分值收敛到小数后若干位稳定，具体数值随 eps 和初始步长略有差异。核对结果时可以拿一个能解析积分的函数替换 myfun.m 试跑，比如 f = x.*y，它在一象限矩形上的二重积分有闭式解 (b²−a²)(d²−c²)/4，用这个基准去验证你的外推循环写没写对，比盯着 tan 的结果瞎猜靠谱得多。

3. 三次样条插值：μ/α/d 三对角组装、追赶法求 M 与车门曲线 S(i+0.5)

3.1 从公式到三对角方程组

第一型三次样条的意思是给定两个端点的一阶导数值 y₀′ 和 yₙ′，其余节点只给函数值。设节点 x₀ < x₁ < … < xₙ，步长 hⱼ = xⱼ₊₁ − xⱼ，二阶导数（弯矩）Mⱼ 满足：

• 内点：μⱼMⱼ₋₁ + 2Mⱼ + αⱼMⱼ₊₁ = dⱼ，j = 1…n−1
• 其中 μⱼ = hⱼ₋₁/(hⱼ₋₁+hⱼ)，αⱼ = 1 − μⱼ = hⱼ/(hⱼ₋₁+hⱼ)
• dⱼ = 6·f[xⱼ₋₁,xⱼ,xⱼ₊₁] = 6·((yⱼ₊₁−yⱼ)/hⱼ − (yⱼ−yⱼ₋₁)/hⱼ₋₁)/(hⱼ₋₁+hⱼ)

两个边界条件补两行：

• 左端：2M₀ + M₁ = (6/h₀)·((y₁−y₀)/h₀ − y₀′)
• 右端：Mₙ₋₁ + 2Mₙ = (6/hₙ₋₁)·(yₙ′ − (yₙ−yₙ₋₁)/hₙ₋₁)

一共 n+1 个未知数、n+1 个方程，系数矩阵是三对角且严格对角占优，用追赶法解最省事。原文档里提到「公式中的 lamda 没有找到对应符号，用 α 代替」，其实就是这个 αⱼ，数学书写里常记作 λⱼ，代码里叫什么无所谓，别和 μ 弄混就行。

把上面的推导落成 MATLAB，第一段是组装：

function M = spline3(x, y, dy0, dyn) % 第一型三次样条：返回各节点的二阶导数 M n = length(x) - 1; % 区间数，节点数为 n+1 h = diff(x); % h(j) = x(j+1) - x(j)，长度 n % 差商：f[x(j-1),x(j)] 与 f[x(j),x(j+1)] f1 = diff(y) ./ h; % 一阶差商，长度 n % 组装三对角：a 下对角、b 主对角、c 上对角、r 右端 a = zeros(n+1,1); b = zeros(n+1,1); c = zeros(n+1,1); r = zeros(n+1,1); % 左边界：2*M0 + M1 = rhs b(1) = 2; c(1) = 1; r(1) = 6/h(1) * (f1(1) - dy0); % 内点 j = 1..n-1（MATLAB 下标 2..n） for j = 2:n mu = h(j-1) / (h(j-1) + h(j)); alp = h(j) / (h(j-1) + h(j)); d = 6 * (f1(j) - f1(j-1)) / (h(j-1) + h(j)); a(j) = mu; b(j) = 2; c(j) = alp; r(j) = d; end % 右边界：M(n-1) + 2*Mn = rhs a(n+1) = 1; b(n+1) = 2; r(n+1) = 6/h(n) * (dyn - f1(n)); M = thomas(a, b, c, r); % 追赶法求解 end

逻辑上有两处必须对齐下标：h 的长度是 n（区间数），f1 长度也是 n，而 M 的长度是 n+1（节点数）。内点循环跑 j = 2…n，正好覆盖 M₁ 到 Mₙ₋₁，边界两行分别占 M₀ 和 Mₙ 的位置。参数 dy0、dyn 就是题面的 y₀′=0.8、y₁₀′=0.2。

3.2 追赶法：不调 \ 也能解

MATLAB 里直接写 M = A\b 是最省事的，但作业的意图是让你把追赶法写出来。三对角方程组的追赶法就是先消元再回代，O(n) 复杂度：

function x = thomas(a, b, c, r) % 追赶法解三对角方程组：a 下对角、b 主对角、c 上对角、r 右端 n = length(r); cp = zeros(n,1); dp = zeros(n,1); cp(1) = c(1) / b(1); dp(1) = r(1) / b(1); % 追：前向消元 for i = 2:n m = b(i) - a(i) * cp(i-1); % 主元 cp(i) = c(i) / m; dp(i) = (r(i) - a(i) * dp(i-1)) / m; end % 赶：回代 x = zeros(n,1); x(n) = dp(n); for i = n-1:-1:1 x(i) = dp(i) - cp(i) * x(i+1); end end

需要注意 a(1) 和 c(n) 虽然用不到，但数组必须存在，传零进去更安全。这段代码和直接 A\b 的结果应当在 1e-12 量级内一致，写完先拿一个 3×3 的小系统手算对照，通过之后再上 11 个节点的车门数据。

3.3 回代求 S(x) 与打印 S(i+0.5)

拿到 M 之后，每个小区间 [xⱼ, xⱼ₊₁] 上的三次多项式写成：

S(x) = Mⱼ(xⱼ₊₁−x)³/(6hⱼ) + Mⱼ₊₁(x−xⱼ)³/(6hⱼ) + (yⱼ − Mⱼhⱼ²/6)(xⱼ₊₁−x)/hⱼ + (yⱼ₊₁ − Mⱼ₊₁hⱼ²/6)(x−xⱼ)/hⱼ

对应的求值函数用一个循环定位区间再套公式：

function s = evalS(x, y, M, xq) % 在查询点 xq 处计算三次样条 S(xq) n = length(x) - 1; s = zeros(size(xq)); for t = 1:numel(xq) xv = xq(t); j = find(x <= xv, 1, 'last'); % 定位所在区间 if j > n, j = n; end % 右端点保护 h = x(j+1) - x(j); A = (x(j+1) - xv)^3 / (6*h); B = (xv - x(j))^3 / (6*h); C = (y(j) - M(j) *h^2/6) * (x(j+1) - xv) / h; D = (y(j+1) - M(j+1)*h^2/6) * (xv - x(j)) / h; s(t) = M(j)*A + M(j+1)*B + C + D; end end

主程序把题面数据填进去，按 0.5 的间隔求值并打印：

x = 0:10; y = [2.51 3.30 4.04 4.70 5.22 5.54 5.78 5.40 5.57 5.70 5.80]; dy0 = 0.8; dyn = 0.2; % 端点导数条件 M = spline3(x, y, dy0, dyn); % 求弯矩 xq = 0.5:1:9.5; % i+0.5, i=0..9 sq = evalS(x, y, M, xq); for i = 1:numel(xq) fprintf('S(%.1f) = %.6f\n', xq(i), sq(i)); end

跑完先做两个自检：把 xq 换成原节点 x，算出的 S(xⱼ) 应当逐个等于 yⱼ，误差在 1e-10 以内；把端点导数换成有限差分近似，结果应当和现在相差不大。两条都过了，再去看 S(i+0.5) 的数值，那才是题面要交的输出。

4. 常微分方程初值问题：RK4 起步、AB4-AM4 预测校正与四种格式的误差排序

4.1 右端函数与精确解各占一个文件

题面给的初值问题是 y′ = −x²y²，0 ≤ x ≤ 1.5，y(0) = 3，步长 h = 0.1，精确解 y(x) = 3/(1+x³)。把右端函数和精确解拆成两个文件，是为了让四种方法共用同一套输入：

function z = fxy(x, y) % 右端函数 f(x,y) = -x^2 * y^2 z = -x.^2 .* y.^2; end function z = fx(x) % 精确解，用于误差统计 z = 3 ./ (1 + x.^3); end

这看起来是小事，但四种方法各写一遍 −x²y²，改题的时候就要改四个地方，漏一个就会出现「三种方法对、一种方法离谱」的假象，排查半天发现是抄错了。

4.2 RK4：四阶单步方法

经典四级四阶 Runge-Kutta：

yᵢ₊₁ = yᵢ + h/6·(k₁ + 2k₂ + 2k₃ + k₄)

k₁ = f(xᵢ, yᵢ)，k₂ = f(xᵢ + h/2, yᵢ + h/2·k₁)，k₃ = f(xᵢ + h/2, yᵢ + h/2·k₂)，k₄ = f(xᵢ + h, yᵢ + h·k₃)

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

单步方法的好处是自启动，只要给 y(0) 就能一路算到底，不需要额外提供前几个点的值。AB4 系列恰恰相反，四阶格式要用到前四步的 f 值，必须先用 RK4 把前 4 个点算出来。

4.3 AB4 显式与 AM4 隐式的预测校正

四阶 Adams-Bashforth 显式格式：

yᵢ₊₁ = yᵢ + h/24·(55fᵢ − 59fᵢ₋₁ + 37fᵢ₋₂ − 9fᵢ₋₃)

四阶 Adams-Moulton 隐式格式：

yᵢ₊₁ = yᵢ + h/24·(9fᵢ₊₁ + 19fᵢ − 5fᵢ₋₁ + fᵢ₋₂)

隐式里带 fᵢ₊₁ = f(xᵢ₊₁, yᵢ₊₁)，未知数在右边，所以先用显式预测一个 ŷᵢ₊₁，再代入隐式做一次校正，这就是预测-校正法（PECE）：

function [x, y] = ab4am4(f, a, b, y0, h, mode) % mode: 'ab4' 纯显式 | 'pc' 预测校正 | 'ipc' 改进的预测校正 x = (a:h:b)'; N = length(x); y = zeros(N,1); y(1) = y0; % 用 RK4 提供前 4 个点 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 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)); % 预测：AB4 显式 yp = y(i) + h/24*(55*f0 - 59*f1 + 37*f2 - 9*f3); switch mode case 'ab4' y(i+1) = yp; case 'pc' % 校正：AM4 隐式，用 yp 代 f(i+1) yc = y(i) + h/24*(9*f(x(i+1), yp) + 19*f0 - 5*f1 + f2); y(i+1) = yc; case 'ipc' % 改进：预测值与校正值加权 yc = y(i) + h/24*(9*f(x(i+1), yp) + 19*f0 - 5*f1 + f2); y(i+1) = 251/270*yc + 19/270*yp; end end end

三处细节值得拎出来说。第一，h/24 这个系数：四阶 Adams 公式的标准形式分母就是 24，写成 h/12 再配 55 的系数会让步长整体偏大一倍，程序不会报错，只是结果差了量级，这种坑在避坑章里还会再提。第二，前 4 个点用 RK4 算出来之后，循环从 i = 4 开始（MATLAB 从 1 开始计数，对应数学上的第 4 步），此时 f₀、f₁、f₂、f₃ 都拿得到。第三，改进格式里的权重 251/270 和 19/270 是题面给的，校正式权重约 0.93，预测式约 0.07，本质上是把校正值再往预测值方向拉一点，抑制预测校正系统的振荡。

4.4 误差统计与方法排序

主程序把四种方法都跑一遍，跟精确解对比，用 2-范数衡量整体误差：

h = 0.1; a = 0; b = 1.5; y0 = 3; [x, ye] = odeEulerGrid(a, b, h); % 生成网格并算精确解 ye = fx(x); [~, y1] = rk4 ( @fxy, a, b, y0, h); [~, y2] = ab4am4(@fxy, a, b, y0, h, 'ab4'); [~, y3] = ab4am4(@fxy, a, b, y0, h, 'pc'); [~, y4] = ab4am4(@fxy, a, b, y0, h, 'ipc'); e = @(y) norm(y - ye); % 2-范数误差 fprintf('RK4 err=%.6e\n', e(y1)); fprintf('AB4 err=%.6e\n', e(y2)); fprintf('AB4-AM4 err=%.6e\n', e(y3)); fprintf('改进型 err=%.6e\n', e(y4));

原文档给出的误差分别是 3.8e-5、4.809e-3、9.2e-4、5.31e-4，排序是 RK4 > 改进型 AB4-AM4 > AB4-AM4 > AB4。这个排序很符合理论预期：RK4 是四阶单步，每步都重新评估斜率，稳定性最好；AB4 是显式多步，误差累积快；加一步校正能把误差压下去一个量级；再对预测值和校正值加权，还能再降一点。把这张表和误差数字一起看，比只背「RK4 精度高」有用得多。

5. 常见问题排查：五个不报错但结果全错的坑

5.1 坑一：Adams 公式系数分母写成 12

现象：AB4 和 AB4-AM4 的结果都偏大，跟精确解的趋势一致但数值系统性偏高，RK4 却完全正常。

原因：四阶 Adams-Bashforth 的标准形式是 yᵢ₊₁ = yᵢ + h/24·(55fᵢ − 59fᵢ₋₁ + 37fᵢ₋₂ − 9fᵢ₋₃)。有些推导把括号内系数整体乘 2 写成 h/12·(110fᵢ − 118fᵢ₋₁ + 74fᵢ₋₂ − 18fᵢ₋₃)，也等价。但只改了分母、没改括号内系数，等价于把步长放大了一倍，误差自然上升一个量级，而且程序不会报任何错。

解决：写完公式先拿一个能解析求解的简单方程对一遍，比如 y′ = y、y(0) = 1，精确解是 eˣ，跑十步看误差量级是否落在 O(h⁴) 附近。对不上就老老实实回教材核对每一项系数，别凭记忆。

5.2 坑二：样条回代时 M 的下标错位

现象：S(i+0.5) 打出来的数值在部分区间明显跳变，尤其在 i = 6、7 附近，曲线看起来「折了一下」。

原因：回代公式里同时出现 Mⱼ 和 Mⱼ₊₁，MATLAB 里对应 M(j) 与 M(j+1)。如果定位区间时用了 find(x <= xv, 1, 'last') 却把结果直接当成 j+1 用，或者循环里 j 从 1 跑到 n 而没有对右端点做保护，区间和弯矩就会整体错一位，只有端点附近碰巧对得上。

解决：写完求值函数先做自检——把查询点设成原节点，S(xⱼ) 必须等于 yⱼ。这一条能抓住 90% 的下标错位。再检查 xq 打在 x(n+1) 上时有没有越界，加一句 if j > n, j = n; end 兜底。

5.3 坑三：二分循环没有上限，程序卡死

现象：exp5.m 跑起来之后命令行一直刷 iter，几分钟都不停，最后只能 Ctrl+C。

原因：外推判据用的是相邻两次近似值之差小于 eps，如果 eps 设得比机器精度还小（比如 1e-16），或者被积函数在积分区域内有接近奇点的行为，差值永远降不到阈值以下。

解决：循环里加 maxit 上限，并且每轮把 diff 打印出来。正常收敛时 diff 应当每轮下降大约一个数量级（复化梯形二分一次误差降 4 倍，外推后更快）；如果 diff 在某个量级来回震荡，说明已经到浮点精度极限，直接把当前值当结果，别硬等。

5.4 坑四：复化梯形端点权重漏乘 1/2

现象：积分结果和参考值差一个小量，而且二分次数增加后误差不按预期下降，一直卡在 1e-3 左右。

原因：复化梯形公式里首尾两个节点的权重是 1/2，中间节点是 1。循环写 for ix = 0:nx 时如果忘记判断端点，全部按 1 累加，等价于多加了两端各半个函数值，误差量级正好是 O(h)。

解决：在内外两层循环里都显式写端点判断，或者更稳妥地用向量化写法：w = ones(nx+1,1); w([1 end]) = 0.5; s = h * sum(w .* fvals)。写完之后拿 f = x.*y 这种有闭式解的函数验证一遍，误差应该随 h 减小而快速下降。

5.5 坑五：AB4 起步的索引从 0 还是从 1 开始

现象：AB4 的结果从第二步开始就明显偏离，而 RK4 完全正常。

原因：数学推导里 x₀ 是第一个点，f₀ 是第一步的斜率；MATLAB 数组从 1 开始，f(1) 对应数学上的 f₀。如果循环写 for i = 4:N-1 却把 f1 取成 f(i-3)、f2 取成 f(i-2)，整个差商窗口就整体右移了一位，前几步看起来还行，越往后偏得越多。

解决：写循环之前先在纸上列一张对照表，把数学下标 0、1、2、3 和 MATLAB 下标 1、2、3、4 写清楚，再动手。跑完之后单独打印前五步的 y 值，跟 RK4 的结果逐点比，前四步必须完全一致（因为就是 RK4 算的），第五步开始才允许有微小差异。

6. 二次开发：把三份作业拼成一个可复用的小型数值算法库

三份作业各自能跑，但接口风格不统一：样条那边传 x、y、dy0、dyn，ODE 这边传函数句柄和步长，重积分又传四个边界加两个步长。想把这套东西变成以后能反复用的小工具，第一件事是统一调用方式。

我的做法是保留底层的单个算法文件，再加一层薄薄的调度脚本，用字符串选择方法，返回值固定成结构体：

function res = numlab(task, opt) % 统一入口：task 取 'spline' | 'integral' | 'ode' switch task case 'spline' M = spline3(opt.x, opt.y, opt.dy0, opt.dyn); res.M = M; res.sq = evalS(opt.x, opt.y, M, opt.xq); % 查询点上的函数值 case 'integral' res.I = integrate2d(opt.f, opt.ab, opt.cd, opt.eps, opt.maxit); case 'ode' res.y = odeSolve(opt.f, opt.span, opt.y0, opt.h, opt.method); otherwise error('未知任务类型: %s', task); end end

odeSolve 内部再用 switch 分发到 rk4、ab4am4 的不同 mode，这样以后加一个新方法只需要在一个地方加分支。第二件事是把误差统计做成表格，四种方法跑完直接出一张对照表，比看一屏 fprintf 直观：

这张表比单纯背结论有用：RK4 精度最高但每步要算四次 f，AB4 只要一次，预测校正两次。实际工程里如果 f 的计算很贵（比如每次都要解一个子问题），AB4-AM4 这类多步法反而更划算，精度和 RK4 差一个量级，代价只有一半。第三件事是给每个函数加一个自检开关，默认关掉，调试时打开就跑一遍基准算例：样条用线性函数 y = 2x+1 验证（三次样条对线性函数应当精确复现），积分用 x.*y 验证，ODE 用 y′ = y 验证。基准过了再上真实数据。

从那以后我每次写完数值代码，都强制先跑一遍能解析求解的基准算例，再去看真实问题的结果。翻车的次数没少，但至少能立刻分清是公式抄错了还是数据本身有问题，排查时间从半天缩到十几分钟。这份上机资料里的三个主脚本和对应子函数都在包里，接手之后建议先跑通 exp5 看输出格式，再换自己的函数进去试，最后再回头看样条和 ODE 那两块，顺序反了容易被一堆下标绕晕。希望帮到你。

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