简介：面向数值分析课程学习与上机实践的MATLAB作业文档，对应东南大学出版教材中的典型上机题，适合数学、计算机及工科学生理解数值计算方法的编程实现。内容围绕三次样条插值、重积分计算和常微分方程初值问题数值解展开，含汽车门曲线型值点的样条插值、复化梯形公式与外推、RK4、AB4、AB4-AM4及改进预测校正等方法的实现与结果比较。文档按题目、算法、运行结果与感悟反思组织，呈现插值条件与连接条件、积分精度控制、初值问题步长设置等细节。压缩包共1个docx文件，约566KB，以文字、公式和结果整理为主，便于对照课程进度查阅。已有905人学习下载。读者可从中获得三类上机题的建模思路、算法公式、MATLAB程序组织方式、误差范数对比结论及排错借鉴，适合课程作业参考与方法复盘。

1. 数值分析上机题拆开看：三道题、四套程序、一份能对着改的 MATLAB 骨架

数值分析上机题这类作业，翻车往往不在推导，而在临门一脚：三弯矩方程抄对了吗，二重积分二分到第八次为什么还不收敛，RK4 和 AB4 的误差为什么差两个数量级。这份资料把三块内容串成一条完整的线——3 次样条插值函数（第一型边界、车门曲线型值点、逐个打印 S(i+0.5)）、重积分的计算（复化梯形 + 逐次二分 + 外推）、常微分方程初值问题数值解（RK4 / AB4 / AB4-AM4 / 改进预测校正），配套的是能直接跑的 MATLAB M 文件而不是伪代码。适合正在补数值分析上机、要交一份跑得通又讲得清的结果的人，也适合把 MATLAB 当日常工具但没系统写过数值程序的人。下面按三道题的顺序拆，每个公式都落到具体行号和参数上，最后给一套用观测阶反推代码正确性的验证办法。

2. 三次样条插值：从三弯矩方程装到车门曲线 S(i+0.5)

2.1 为什么这里的端点条件必须是第一型

三次样条在每个小区间上是三次多项式，要把它唯一确定下来，除了节点函数值和一阶、二阶导数连续这两类连接条件，还差两个端点约束，这就是所谓的第一型、第二型、周期型三选一。第二型（自然边界）直接令两端二阶导为零，写起来最省事，但代价是端点附近的弯矩被强行摁成 0，曲线在靠近端点的一两个区间内容易出现轻微的外翻，做车门轮廓这种端点斜率有明确物理含义的曲线时会很难看。周期型只对首尾闭合的型值点成立，这里 x 从 0 到 10 是一条开口曲线，用不上。

题目给的端点条件是 y0' = 0.8、y10' = 0.2，属于典型的第一型（固支）边界：直接指定端点一阶导数，让样条在端点处的走向和设计意图一致。这一类边界的方程组仍然是三对角、严格对角占优，直接左除就能解，不需要迭代。选型的判断标准其实就一句话——端点的一阶导你有没有可信来源。有，就用第一型；没有且曲线是自然延伸的，退而求其次用第二型；闭合轮廓才轮到周期型。

还有一处符号上的麻烦要提前说清楚：不同教材里 μ 和 λ 的记法不统一，有的书写μ_j M_{j-1} + 2M_j + λ_j M_{j+1} = d_j，有的书把后面那个系数印成 α。这份笔记里用的是 α，α_j = 1 - μ_j，含义完全一样。抄公式的时候别被字母带偏，认准它是"右邻居的权重"就行。

2.2 三弯矩矩阵怎么装：h、μ、α、d 四个量逐个对

三弯矩法的思路是先把每个节点上的二阶导 M_j 解出来，再回代成每个区间上的三次多项式。核心方程对内部节点 j = 1, 2, …, n-1 成立：

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

四个量的定义和车门数据（x = 0..10，共 11 个点，等距）下的取值如下：

等距是个很舒服的特殊情况：μ 和 α 全部退化成 0.5，矩阵变成[0.5, 2, 0.5]的常系数三对角，只有首尾两行被边界条件改写。d_j 里的 f[·,·,·] 是二阶差商，等距时等价于(y_{j+1} - 2y_j + y_{j-1}) / (2h^2)，乘 6 之后就是3(y_{j+1} - 2y_j + y_{j-1})/h^2。写成差商形式更通用，因为一旦换成不等距型值点，这个式子不用改。

第一型边界的首尾两行是这么来的：

• 左端：2M_0 + M_1 = (6/h_0)·((y_1 - y_0)/h_0 - y_0')
• 右端：M_{n-1} + 2M_n = (6/h_{n-1})·(y_n' - (y_n - y_{n-1})/h_{n-1})

手工算几个数当校验用，比跑完再看曲线靠谱得多。y = [2.51, 3.30, 4.04, 4.70, 5.22, 5.54, 5.78, 5.40, 5.57, 5.70, 5.80]：

• d_0 = 6·(0.79 - 0.8) = -0.06
• d_1 = 6·(4.04 - 2×3.30 + 2.51)/2 = -0.15
• d_9 = 6·(5.80 - 2×5.70 + 5.57)/2 = -0.09
• d_10 = 6·(0.2 - 0.10) = 0.60

如果程序打印出来的 d 向量头尾不是这四个数，基本可以断定边界行的系数或者差商分母出了问题，不用再往下看结果。这一步花两分钟，能省掉后面半小时的自我怀疑。

2.3 完整 MATLAB 实现：装配、回代、分段求值

装配和求解写在一个函数里，求值单独拆出来，这样后面做自检、画图、换数据都方便。注意 MATLAB 下标从 1 开始，x(1) 对应 x_0，所有公式里的 j 都要整体平移一位。

function [M, S] = spline3_fixed(x, y, d0, dn) % SPLINE3_FIXED 第一型（固支）三次样条插值 % x, y : 型值点，长度 n+1，等距或不等距都可以 % d0, dn: 左右端点一阶导数 y'(x0)、y'(xn) % M : 弯矩向量，长度 n+1 % S : 函数句柄，输入标量或向量，返回插值结果 n = numel(x) - 1; h = diff(x); % h(j) = x(j+1) - x(j)，长度 n A = zeros(n + 1, n + 1); rhs = zeros(n + 1, 1); % ---- 左端点：2M0 + M1 = (6/h1) * ((y2-y1)/h1 - d0) ---- A(1,1) = 2; A(1,2) = 1; rhs(1) = 6 / h(1) * ((y(2) - y(1)) / h(1) - d0); % ---- 内部节点 j = 2..n，对应数学下标 1..n-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; % 二阶差商乘 6，直接展开成三点形式，省掉差商表 rhs(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} + 2Mn = (6/h_n) * (dn - (y_{n+1}-y_n)/h_n) ---- A(n+1, n) = 1; A(n+1, n+1) = 2; rhs(n+1) = 6 / h(n) * (dn - (y(n+1) - y(n)) / h(n)); M = A \ rhs; % 三对角正定，直接左除足够稳 S = @(t) spline_eval(t, x, y, M, h); end function v = spline_eval(t, x, y, M, h) % 在 t 上逐点定位区间并求值，t 可以是向量 v = zeros(size(t)); n = numel(h); for k = 1:numel(t) xi = t(k); j = find(x <= xi, 1, 'last'); % 落在哪一段 j = min(max(j, 1), n); % 右端点处钳位，防止越界 hj = h(j); a = (x(j+1) - xi) / hj; % 左端基函数 b = (xi - x(j)) / hj; % 右端基函数 % 等价于教材的分段表达式，但写成 a^3-a 的形式数值更稳 v(k) = M(j) * hj^2 / 6 * (a^3 - a) ... + M(j+1) * hj^2 / 6 * (b^3 - b) ... + y(j) * a + y(j+1) * b; end end

分段求值那三行是整个函数里最容易写错的地方。教材上的标准形式是

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^2/6)(x_{j+1}-x)/h_j + (y_{j+1} - M_{j+1} h_j^2/6)(x-x_j)/h_j

四项展开合并之后，含 y 的两项正好是y_j·a + y_{j+1}·b，含 M 的两项是M_j h_j^2/6·(a^3-a)。合并写法少两次乘法和一次减法，在x接近x_j时也不会因为两个大数相减丢精度。

参数上还有两个可调项值得留意。A \ rhs走的是 LU 分解，对 11 阶矩阵来说是杀鸡用牛刀，真要在几千个节点上跑，换成追赶法（Thomas 算法）能把复杂度从 O(n^3) 压到 O(n)，不过上机题这个规模完全没必要。spline_eval里的find是线性扫描，对每个 t 都重扫一遍，向量化之后可以先histc或discretize一次性算出所有点的区间号；节点数上三位数以后再考虑这个优化。

2.4 S(i+0.5) 打印前的两道自检

程序跑通不等于结果对。样条有两类必须满足的硬约束，用它们做断言比看曲线敏感得多。

第一道是节点插值条件：S(x_i) = y_i必须严格成立（差到 1e-12 量级）。这一条不过，说明分段求值的区间定位或基函数写反了。第二道是端点一阶导：用中心差分数值求 S'(x_0) 和 S'(x_n)，看是否逼近 0.8 和 0.2。差分的步长取 1e-4 就够，再小会被舍入误差淹没。这两条都过了，再打印 S(i+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]; [M, S] = spline3_fixed(x, y, 0.8, 0.2); % 自检 1：节点处必须严格等于型值点 assert(norm(S(x(:))' - y) < 1e-12, '节点插值条件不满足'); % 自检 2：端点一阶导用中心差分粗验 hs = 1e-4; dL = (S(x(1) + hs) - S(x(1))) / hs; dR = (S(x(end)) - S(x(end) - hs)) / hs; fprintf('左端斜率 %.6f（目标 0.8），右端斜率 %.6f（目标 0.2）\n', dL, dR); % 正式输出：每个区间的中点 tm = (0:9) + 0.5; fprintf('%8s %12s\n', 'x', 'S(x)'); fprintf('%8.1f %12.6f\n', [tm; S(tm)]);

打印出来的表就是题目要的 S(i+0.5)，i = 0..9。这里有个容易被忽略的细节：S(i+0.5) 落在第 i+1 个区间内，而三次样条只在节点处保证函数值精确，区间内部是有截断误差的，量级取决于四阶导数和 h^4。所以这十个数的意义是"光滑逼近"，不是"精确值"，写报告的时候别把它们当成真值去和解析解对。要判断这组数靠不靠谱，真正的依据是前面两道自检，加上 2.2 节里那几个手算的 d 值。

3. 重积分与外推：复化梯形逐次二分，把 O(h²) 压到 O(h⁴)

3.1 Richardson 外推为什么能白捡两阶

复化梯形公式的误差有渐近展开：T(h) = I + c₁h² + c₂h⁴ + …。步长减半，误差主项直接缩到四分之一。把两个相邻的梯形值按 4:1 反号组合：

S = (4·T(h/2) - T(h)) / 3

代入展开式，c₁h²项正好被消掉，剩下的是 O(h⁴)——这个S就是复化 Simpson 公式。教材上写的S_n = 4/3·T_{2n} - 1/3·T_n是同一个式子的另一种排列，两者只差一个写法。

这个操作的本质是"用两次低阶结果构造一次高阶结果"，代价是梯形值必须配对：T(h)和T(h/2)得对应同一套网格的嵌套关系。放到二重积分上，这一点就变成了硬要求——x 方向和 y 方向的步长必须同步减半。二维复化梯形的误差是c₁(h² + k²)量级，如果只把 h 减半、k 不动，主项里的k²不会被消掉，外推出来的结果仍然是二阶精度，看起来"精度不涨"，其实是外推前提被破坏了。这是我在这一章里最想强调的一条边界。

被积函数tan(x² + y²)在积分区域x ∈ [0, π/3]、y ∈ [0, π/6]上没有奇点：x² + y² 的最大值取在角点，等于π²/9 + π²/36 = 5π²/36 ≈ 1.3708，而 tan 的奇点在 π/2 ≈ 1.5708，还差着一大截。函数在区域上的最大值约等于 tan(1.3708) ≈ 4.99，被积函数有界且光滑，复化梯形不会发散，外推的渐近展开也成立。这一点在做题时最好写进报告，它解释了为什么这套方法在这里能用。

3.2 getT 与被积函数：两个子函数的分工

结构上按笔记里的拆法是合理的：把"在给定网格上算一次复化梯形值"写成独立子函数（原文叫 getT.m），把被积函数单独写一个（myfun.m），主程序只负责二分和外推的控制流。这样换一道积分题，只需要动 myfun.m 一个文件。

二维复化梯形的权重是张量积形式：角点权重 1，边上非角点权重 2，内部点权重 4，最后乘h·k/4。用meshgrid一次性算出所有网格点的函数值，再和权重矩阵点乘，比双重 for 循环快一个数量级，也短得多。

function [val, tab] = dbl_trap(f, ax, bx, cy, dy, tol, maxit) % 二重积分：复化梯形 + 逐次二分 + Richardson 外推（等价于复化 Simpson） % f : 被积函数句柄 f(x, y)，要求支持向量化输入 % ax,bx : 内层 x 的积分区间 % cy,dy : 外层 y 的积分区间 % tol : 相邻两次外推近似值之差的容限 ε % maxit : 最大二分次数，防止死循环 % tab : [T_n, S_n, |S_n - S <p> <a href="https://download.csdn.net/download/qq_35293123/11071944" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>