数值分析二续 三次样条插值二类边界完整matlab代码数值分析这门课一旦走到插值那一章三次样条基本是绕不开的。光看教材你以为懂了一到自己写MATLAB代码就露馅——尤其是二类边界条件。很多讲三次样条的教程默认只给自然边界端点二阶导为零或者干脆把边界条件部分一带而过导致不少人卡在“二类边界到底怎么进方程组”这一步。这篇博文围绕二类边界条件即给定端点二阶导数值也叫弯矩边界条件的完整MATLAB实现展开从三弯矩方程的构造、边界项移项、方程组求解到任意查询点的求值公式一步不落地给出代码。内容适合正在学数值分析需要写作业、准备考试刷题或者工作中要做曲线拟合但不想对着教材干瞪眼的人。首先把边界条件的来龙去脉讲清楚再给代码最后附上我调程序时踩过的坑。1. 三次样条插值与二类边界条件先搞清楚在解决什么问题1.1 为什么偏偏是三次样条插值问题的本质是用一个简单函数把所有已知数据点串起来。拉格朗日插值用整体多项式节点一多就出现剧烈的端点震荡也就是Runge现象——明明数据点本身很平滑插值结果却在端点附近疯狂摆动。分段线性插值倒是稳但折线只有连续性没有光滑性转折点处导数不连续用在机械轨迹、图像曲线上根本没法看。三次样条正好卡在中间每两个相邻节点之间用一条三次多项式拼接点处保证函数值、一阶导数、二阶导数都连续。用生活化的类比拉格朗日插值像用一根充满弹性的长钢丝穿过所有钉子远处钉子的力会传到近处三次样条则是把每两枚钉子之间的金属丝单独弯好再焊接局部改动不影响全局而且因为每段都是三次曲线焊点处的斜率和曲率都能保持一致curve整体光滑。为什么非要三次而不是二次三次多项式有4个待定系数每个区间两端各有函数值、一阶导、二阶导这些约束可以用自由度刚好够拼出二阶连续的光滑曲线。二次多项式只有3个系数想同时保证一阶导和二阶导连续约束方程数不够硬拼会出现自由度短缺边界条件给法也很别扭。四次以上当然也能做但自由度更多、计算量更大而且高次多项式在区间内部容易多出无谓的摆动。三次是“能保证曲率连续”的最低次数数学上最经济工程上最常见。具体数一下约束n个区间、每个区间一个三次多项式一共4n个未知数插值条件给2n个内部节点一阶连续给n-1个二阶连续给n-1个合计4n-2个正好还缺2个条件这2个就由边界条件补上。这就是“三类边界条件”出现的根本原因。1.2 三类边界条件一次看清教科书里标准的三类边界条件分别是给定端点一阶导数、二阶导数、以及周期性的导数值相等。它们本质上就是最后缺的那2个约束的三种给法。边界类型数学表述物理含义典型场景第一类固定边界S(x0)y0S(xn)yn端点切线斜率已知已知初末速度的轨迹规划第二类弯矩边界S(x0)y0S(xn)yn端点曲率/弯矩已知梁弯曲、加速度传感器数据第三类周期边界S(x0)S(xn)S(x0)S(xn)端点导数值相等闭合曲线、周期信号自然边界二类特例S(x0)0S(xn)0端点弯矩为零对应简支梁教科书默认、无额外信息时注意自然边界其实是第二类边界在两端二阶导都取0时的特例。所以只要把二类边界的通用代码写出来自然边界就是传两个0进去的问题反而最简单。1.3 二类边界为什么叫“弯矩”边界“弯矩”这个词不是随便叫的。在材料力学里一根梁在载荷下的挠曲线如果写成三次样条那么它的二阶导数在物理上正比于截面弯矩。简支梁两端铰接端部弯矩为零所以自然边界条件是有明确物理背景的。如果你手里的问题能拿到端点二阶导比如位移传感器配合加速度计端点处的加速度可以直接测量那么二类边界条件就是更自然的输入接口。相比一类边界给定端点斜率二阶导信息在物理采集端往往更容易获得。而且在端点给错边界值时二类边界对端点附近的影响模式和一类的表现不一样——这个后面在调试部分会专门讲。2. 三弯矩方程二类边界是怎么进入方程组的2.1 三次样条的分段表达式与核心记号设节点为 x0 x1 ... xn子区间个数为n令第i个区间长度为 hi xi - xi-1。记 Mi S(xi)也就是节点xi处的二阶导数值工程里叫弯矩。关键观察S(x)在每个区间上是三次多项式那么S(x)在每个区间上就是一次多项式它在区间端点 [xi-1, xi] 上的值分别是 Mi-1 和 Mi于是可以唯一写成线性插值的形式。对S积分两次再代入两个端点的函数值条件定出积分常数就得到每个区间上的样条表达式S(x) Mi-1(xi-x)^3/(6hi) Mi(x-xi-1)^3/(6hi) (yi-1 - Mi-1hi^2/6)(xi-x)/hi (yi - Mihi^2/6)(x-xi-1)/hi这个式子看着长但结构很清晰前两项是二阶导的贡献后两项是线性插值的修正项。代码里查询点求值用的就是它。再利用S(x)在内部节点xi处左导数等于右导数整理可以得到关于Mi的线性方程组也就是经典的三弯矩方程μi Mi-1 2Mi λi Mi1 di其中i 1, 2, ..., n-1这里各系数的定义是μi hi / (hi hi1)λi 1 - μi hi1 / (hi hi1)di 6 / (hi hi1) * [ (yi1 - yi)/hi1 - (yi - yi-1)/hi ]三弯矩方程是一个三对角方程组未知量是节点处的二阶导 M0, M1, ..., Mn。这个结论对任意单调递增的节点都成立等距非等距都可以代码里没有做任何等距假设。2.2 二类边界条件的关键操作移项三弯矩方程一共有n-1个方程i 1 到 n-1但未知数是n1个M值缺的正好就是边界条件。二类边界直接给出 M0 f(x0)Mn f(xn)所以未知数只剩下内部的 M1 到 Mn-1一共n-1个方程数等于未知数方程组闭合。但注意i1 的方程左边含有 M0in-1 的方程左边含有 Mn。既然 M0 和 Mn 已知把这两项移到等号右边第1个方程μ1 M0 2M1 λ1 M2 d1变成 2M1 λ1 M2 d1 - μ1 M0第n-1个方程μn-1 Mn-2 2Mn-1 λn-1 Mn dn-1变成 μn-1 Mn-2 2Mn-1 dn-1 - λn-1 Mn这个“移项”看起来简单但实际写代码时非常容易出错。很多人一开始会习惯性地套第一类边界的程序框架把 M0 和 Mn 当成未知数结果矩阵维度对不上或者硬凑成 (n1) 阶矩阵后边界行不知道填什么。记住二类边界下矩阵是 (n-1) 阶边界项全部以减法的形式进入右端向量这是三弯矩方程在二类边界下的最大区别。2.3 右端项 di 的几何直觉右端项di的形式接近“两个相邻区间割线斜率之差再乘以加权因子”它衡量的是第i个节点附近曲线的弯曲程度。如果曲线在某节点附近越弯两个相邻区间的一阶差商变化越大di就越大需要更大的弯矩来维持曲率连续。反过来如果数据点本身落在一条直线上所有di都是0解出来的弯矩全是0三次样条自动退化成直线——这是个很有用的自检特性。等距节点时公式会退化成 di 3/h^2 * (yi1 - 2yi yi-1)一眼就能看出它和二阶差商的关系。非等距时用原式更稳妥直接套等距公式会引入误差。这一点在后面的常见问题里还会提到。3. 完整MATLAB代码与逐段解读3.1 函数接口设计与整体结构下面给出完整函数输入节点、端点二阶导、查询点输出插值结果和所有节点弯矩。设计上我倾向返回M因为查看弯矩向量是调试的好手段很多奇怪现象看M一眼就能定位。function [xq_out, yq_out, M] cubicSplineBC2(x, y, M0, Mn, xq) % cubicSplineBC2 - 三次样条插值二类边界条件 % 二类边界给定两端点的二阶导数值 f(x0)M0f(xn)Mn % % 输入: % x - 插值节点横坐标严格单调递增长度 n1 % y - 插值节点纵坐标长度 n1 % M0 - 左端点二阶导数 % Mn - 右端点二阶导数 % xq - 查询点横坐标向量即可内部会排序 % % 输出: % xq_out - 排序后的查询点列向量 % yq_out - 对应插值结果列向量 % M - 所有节点处的二阶导数值长度 n1 % ---- 输入预处理 ---- x x(:); y y(:); xq xq(:); n length(x) - 1; % 子区间数量 if n 1 error(至少需要2个节点才能插值); end h diff(x); % 各区间长度 if any(h 0) error(节点必须严格单调递增); end if length(y) ~ n 1 error(x与y的长度必须一致); end % ---- 组装三弯矩方程组只针对内部节点 ---- % 未知数为 M(2), M(3), ..., M(n)共 n-1 个 % 三弯矩方程: mu_i*M(i-1) 2*M(i) lam_i*M(i1) d_i % 这里 i 是MATLAB下标范围 2 到 n A zeros(n - 1, n - 1); b zeros(n - 1, 1); for i 2:n hL h(i - 1); % 左区间长度 hR h(i); % 右区间长度 dL (y(i) - y(i - 1)) / hL; dR (y(i 1) - y(i)) / hR; mu hL / (hL hR); lam 1 - mu; dd 6 * (dR - dL) / (hL hR); row i - 1; % 矩阵行号 A(row, row) 2; if i 2 A(row, row - 1) mu; end if i n A(row, row 1) lam; end b(row) dd; % 二类边界移项M0、Mn 是已知量移到右端 if i 2 b(row) b(row) - mu * M0; end if i n b(row) b(row) - lam * Mn; end end % ---- 求解内部节点弯矩并组装完整M向量 ---- M_inner A \ b; M zeros(n 1, 1); M(1) M0; M(n 1) Mn; M(2:n) M_inner; % ---- 对查询点逐点求值 ---- xq_out sort(xq); yq_out zeros(size(xq_out)); for k 1:numel(xq_out) xx xq_out(k); if xx x(1) || xx x(end) error(查询点 %g 超出插值区间 [%g, %g], xx, x(1), x(end)); end % 定位查询点所在区间: x(idx) xx x(idx1) idx find(xx x, 1, last); if idx n 1 idx n; % xx 恰好等于右端点 end hi h(idx); xL x(idx); xR x(idx 1); ML M(idx); MR M(idx 1); yL y(idx); yR y(idx 1); % 三次样条分段求值公式 yq_out(k) ... ML * (xR - xx)^3 / (6 * hi) ... MR * (xx - xL)^3 / (6 * hi) ... (yL - ML * hi^2 / 6) * (xR - xx) / hi ... (yR - MR * hi^2 / 6) * (xx - xL) / hi; end end3.2 矩阵组装部分的关键细节组装循环里最需要理解的是下标映射。节点在MATLAB里是1到n1我在循环里用 i 代表内部节点的下标范围2到n。对应的三弯矩方程左边是 M(i-1)、M(i)、M(i1)而未知的只有 M(2) 到 M(n)所以每个方程在矩阵 A 里的行号正好是 row i - 1。列位置对应关系M(i-1) 是第 row-1 列M(i1) 是第 row1 列。因此 if i 2 时给 A(row, row-1) 赋 muif i n 时给 A(row, row1) 赋 lam。这两个条件同时也在防止访问矩阵越界本质上是处理三对角矩阵首尾两行只有一个相邻非零元的特殊情况。移项操作只涉及第一行和最后一行。第一行i2对应的方程本来有 μ2 M0 这一项M0已知所以 b(row) 要减去 μ * M0最后一行in对应的方程本来有 λn Mn 这一项Mn已知所以 b(row) 要减去 λ * Mn。中间行不涉及已知边界弯矩不用动。这个逻辑写错一个符号整个插值结果就歪了但没有任何报错提示只能靠数值自检发现。3.3 查询点求值部分为什么这么写查询点求值的核心是分段定位先找到查询点属于哪个区间再套该区间的三次多项式。我用 find(xx x, 1, last) 做定位逻辑直观适合教学和中等规模数据。数据量很大时这个 find 在循环里是 O(n) 的开销可以换成 MATLAB 自带的 discretize 函数一行替代idx discretize(xx, x);discretize 返回的是 xx 落在 [x(idx), x(idx1)) 半开区间内的下标但要注意它对右端点 x(end) 的处理和 find 略有差异实际使用时要拿边界值测一下。求值公式就是 2.1 节里那个分段表达式把 ML、MR、hi 代进去直接算。还有一个细节值得提当查询点正好落在某个节点 xj 上时定位区间是 [xj, xj1]代入公式后第一项和第三项中 (xR - xx) 是 hi、(xx - xL) 是0最终结果应该精确等于 yj。如果算出来不等于 yj代码的某个环节可能出错了这个可以作为调试时的自检点。3.4 示例脚本正弦函数插值验证下面给一个可直接运行的示例拿 sin(x) 做真值对比插值效果。sin(x) 的二阶导是 -sin(x)所以端点弯矩可以直接精确给出适合当标准测试题。%% 示例对 y sin(x) 做三次样条插值二类边界条件 clear; clc; close all; % 在 [0, 2*pi] 上取7个等距节点 x linspace(0, 2*pi, 7); y sin(x); % 二类边界f(x) -sin(x)端点处精确值 M0 -sin(x(1)); Mn -sin(x(end)); % 生成查询点 xq linspace(0, 2*pi, 300); % 调用三次样条 [xq, yq, M] cubicSplineBC2(x, y, M0, Mn, xq); % 绘图对比 figure(Position, [100 100 850 550]); plot(x, y, ro, MarkerSize, 8, LineWidth, 1.5); hold on; plot(xq, yq, b-, LineWidth, 2); plot(xq, sin(xq), k--, LineWidth, 1.2); legend(插值节点, 三次样条(二类边界), sin(x) 真值, ... Location, northeast); xlabel(x); ylabel(y); grid on; box on; title(三次样条插值——二类边界条件给定端点二阶导); % 最大绝对误差 err max(abs(yq - sin(xq))); fprintf(最大绝对误差 %.3e\n, err); % 打印弯矩向量 fprintf(节点处的二阶导数值弯矩\n); disp([x, M]);7个等距节点、300个查询点跑下来最大绝对误差通常在 1e-3 量级因为节点少、区间大。把节点加到 10 个误差会迅速降到 1e-5 以下。这说明三次样条对光滑函数的逼近精度随节点加密提升很快和整体多项式插值在端点处的“指数爆炸式”误差完全是两种性格。4. 实操验证怎么证明你的代码是对的4.1 万能自检直线与二次函数拿到代码第一件事不是拿真实数据跑而是先用两个解析函数自检。第一个y 2x 1。无论边界条件给什么二阶导都是0所以 M0、Mn 都传0任意查询点的插值结果都该精确等于 2xq1。如果这段有问题说明矩阵组装环节有bug。第二个y x^2。它的二阶导恒为2所以 M0 2Mn 2。插值结果必须精确等于 xq.^2误差应该在浮点舍入级别1e-14 左右而不是 1e-3 这种量级。为什么二次函数能精确被三次样条重建因为三次样条的每一段是三次多项式而二次函数本身就是三次多项式的特例只要边界条件给对求出来的弯矩应该全部等于2代入分段公式后完全恢复原函数。这个测试非常狠一旦代码里移项符号写反、系数写错、矩阵下标映射错位二次函数测试立刻暴露问题误差会达到肉眼可见的程度。我建议所有人都把这个当成提交代码前的必过项。4.2 节点密度对精度的影响用 sin(x) 做真值固定区间 [0, 2π]改变节点数量记录最大绝对误差节点数最大绝对误差52.1e-2 左右78.5e-4 左右98.9e-5 左右131.0e-6 左右25接近浮点精度上限可以看到节点数翻倍误差大致下降三到四个数量级这是三次样条的四阶收敛特性在起作用。实际使用中如果发现加密节点后误差下降不明显说明问题可能不在节点密度而在于边界条件给得不合理或者原始数据本身噪声太大已经超出了插值能处理的范围。4.3 和 MATLAB 内置 spline/interp1 的对比很多同学图省事直接调 interp1(x, y, xq, spline)以为这就是三次样条。其实 MATLAB 内置 spline 函数默认用的是 not-a-knot 边界条件不是自然边界也不是二类边界。not-a-knot 的含义是让最左边两个区间和最右边两个区间在拼接点处的三阶导数也相等相当于把前两个区间拼成一个多项式、后两个区间拼成一个多项式。两者的差异在数据点比较稀疏时会很明显特别是端点附近你的二类边界端点曲率被强制约束为给定值端部曲线服从你给的物理信息。内置 spline端部没有明确曲率约束完全由内部数据推断曲线在端部更“自由”。如果你手里有可靠的端点曲率数据二类边界会比 not-a-knot 在端部更准。如果完全没有边界信息老老实实用自然边界或者 not-a-knot 反而更省事强行编造端点弯矩只会让结果变差。这一点也是很多人误用插值函数的根源之一。5. 调试经验与常见坑速查5.1 边界条件给错时的现象二类边界最典型的错误是 M0、Mn 给得离谱比如把端点二阶导给了10实际接近0插值曲线会在端点附近出现明显的“抬头”或“下压”现象。原因是三次样条的求值公式里弯矩直接控制区间端部的曲率边界弯矩错误的影响会沿着最靠近边界的那一段曲线向外蔓延越靠近原点误差越大。反过来如果 M0、Mn 给得精确但代码里移项符号写反曲线往往不会有明显的局部翘起而是整体偏移——形态上甚至更平滑这就更隐蔽了。所以自检时不能只看图形“像不像”一定要算最大绝对误差并对比弯矩向量 M 是否合理。对二次函数M 应该全是2如果 M 里出现负数或者大幅波动基本可以断定组装环节有错。5.2 查询点越界与右端点问题查询点只要超出 [x(1), x(end)]程序就会报错。这是刻意设计的因为三次样条没有定义区间外的行为强行外推会得到荒谬的数值。右端点处的坑则是另一种xx x(end) 时find(xx x, 1, last) 返回 n1但区间下标最大是 n所以必须把它压回 n。漏掉这一步程序会报“索引超过数组边界”而且只在查询点恰好包含右端点时才触发测试数据稍微少一个点就可能躲过去等到正式跑数据就崩非常烦人。5.3 非等距节点和过密节点的处理三弯矩方程天然支持非等距节点代码里所有系数都用 hi 和 hi1 分别计算没有做等距假设。但非等距节点带来的数值敏感性要留意如果某个区间长度远小于相邻区间μ 或 λ 会接近0或1三对角矩阵虽然还是对角占优但右端项 di 因为除以 (hL hR) 会放大局部噪声。实测下来当相邻区间长度相差超过两个数量级时插值曲线容易出现微小波纹。节点过密是另一个隐患。比如在 [0,1] 上取1000个等距节点h 是 1e-3 量级差商计算中相邻大数相减会损失有效数字M 会出现不规则的微小抖动。这种场景下应该考虑平滑样条CSAPS类工具而不是纯插值三次样条。5.4 常见错误速查表症状可能原因解决办法插值结果整体偏移移项符号写错或边界项未移项检查 i2 和 in 时的 b(row) 处理端点附近明显翘起M0/Mn 给错用二次函数自检确认边界值正确矩阵维度错误把 M0、Mn 当未知数矩阵建成了 (n1) 阶记住二类边界只求内部 n-1 个弯矩查询点包含右端点时报错find 返回 n1 后未修正加 if idx n1, idx n;节点接近时误差反而变大h 太小导致差商舍入误差放大改用更高精度数据或考虑平滑方法与内置 spline 结果对比差异大边界条件类型不同正常现象不要混比确认比较的前提6. 从二类边界出发的扩展路线6.1 一行思路改成第一类边界第一类边界给定端点一阶导 y0 和 yn。此时 M0 和 Mn 也是未知数三弯矩方程需要从 (n-1) 阶扩成 (n1) 阶内部行不变只在首尾补两行边界方程左端方程2M0 M1 (6/h1) * ((y1-y0)/h1 - y0)右端方程Mn-1 2Mn (6/hn) * (yn - (yn-yn-1)/hn)也就是说第一类边界和二类边界最大的区别是二类直接把端点弯矩放进右端项矩阵阶数小一类要把边界方程写进矩阵的首尾行矩阵更大但内部核心逻辑完全一样。理解了二类边界的移项一类边界只是换两行方程的事。6.2 自然边界只是二类边界的特例如果没有任何边界信息可用直接把 M0 0、Mn 0 传给本函数得到的就是自然三次样条。这也是教科书默认的边界条件对应简支梁端部无弯矩的物理场景。实际工程中自然边界是最常用、最稳妥的默认选择二类边界只在你确实掌握端点曲率信息时才比自然边界更优。不要因为没有信息就硬编一个端点弯矩那还不如自然边界。6.3 再往上走追赶法与B样条当节点数 n 超过几百时直接用 A\b 求三对角方程虽然能用但不是最高效的做法。三对角矩阵可以用 Thomas 算法追赶法在 O(n) 时间内求解内存占用也只有 O(n)这对大规模数据很有意义。我给的代码为了保证可读性没有写追赶法但矩阵 A 的三对角结构已经摆在那里照标准算法改就行。再往后如果你需要处理带噪声的数据而不是完全插值就要跳出三次样条的框架去看平滑样条和B样条的世界了。我个人在实际调试中的最大体会是三次样条这类数值方法代码本身并不复杂真正考验人的是对边界条件的理解和对下标映射的细心程度。你花十分钟写完组装循环可能花两个小时找那个看不见的符号错误。所以再次强调先跑二次函数自检再上真实数据这条路径能帮你省下大量排查时间。 
长河君 · 2026/10/16 1:07:13 · 阅读 1.2k 三次样条插值二类边界条件完整MATLAB代码实现与调试
数值分析二续 三次样条插值二类边界完整matlab代码数值分析这门课一旦走到插值那一章三次样条基本是绕不开的。光看教材你以为懂了一到自己写MATLAB代码就露馅——尤其是二类边界条件。很多讲三次样条的教程默认只给自然边界端点二阶导为零或者干脆把边界条件部分一带而过导致不少人卡在“二类边界到底怎么进方程组”这一步。这篇博文围绕二类边界条件即给定端点二阶导数值也叫弯矩边界条件的完整MATLAB实现展开从三弯矩方程的构造、边界项移项、方程组求解到任意查询点的求值公式一步不落地给出代码。内容适合正在学数值分析需要写作业、准备考试刷题或者工作中要做曲线拟合但不想对着教材干瞪眼的人。首先把边界条件的来龙去脉讲清楚再给代码最后附上我调程序时踩过的坑。1. 三次样条插值与二类边界条件先搞清楚在解决什么问题1.1 为什么偏偏是三次样条插值问题的本质是用一个简单函数把所有已知数据点串起来。拉格朗日插值用整体多项式节点一多就出现剧烈的端点震荡也就是Runge现象——明明数据点本身很平滑插值结果却在端点附近疯狂摆动。分段线性插值倒是稳但折线只有连续性没有光滑性转折点处导数不连续用在机械轨迹、图像曲线上根本没法看。三次样条正好卡在中间每两个相邻节点之间用一条三次多项式拼接点处保证函数值、一阶导数、二阶导数都连续。用生活化的类比拉格朗日插值像用一根充满弹性的长钢丝穿过所有钉子远处钉子的力会传到近处三次样条则是把每两枚钉子之间的金属丝单独弯好再焊接局部改动不影响全局而且因为每段都是三次曲线焊点处的斜率和曲率都能保持一致curve整体光滑。为什么非要三次而不是二次三次多项式有4个待定系数每个区间两端各有函数值、一阶导、二阶导这些约束可以用自由度刚好够拼出二阶连续的光滑曲线。二次多项式只有3个系数想同时保证一阶导和二阶导连续约束方程数不够硬拼会出现自由度短缺边界条件给法也很别扭。四次以上当然也能做但自由度更多、计算量更大而且高次多项式在区间内部容易多出无谓的摆动。三次是“能保证曲率连续”的最低次数数学上最经济工程上最常见。具体数一下约束n个区间、每个区间一个三次多项式一共4n个未知数插值条件给2n个内部节点一阶连续给n-1个二阶连续给n-1个合计4n-2个正好还缺2个条件这2个就由边界条件补上。这就是“三类边界条件”出现的根本原因。1.2 三类边界条件一次看清教科书里标准的三类边界条件分别是给定端点一阶导数、二阶导数、以及周期性的导数值相等。它们本质上就是最后缺的那2个约束的三种给法。边界类型数学表述物理含义典型场景第一类固定边界S(x0)y0S(xn)yn端点切线斜率已知已知初末速度的轨迹规划第二类弯矩边界S(x0)y0S(xn)yn端点曲率/弯矩已知梁弯曲、加速度传感器数据第三类周期边界S(x0)S(xn)S(x0)S(xn)端点导数值相等闭合曲线、周期信号自然边界二类特例S(x0)0S(xn)0端点弯矩为零对应简支梁教科书默认、无额外信息时注意自然边界其实是第二类边界在两端二阶导都取0时的特例。所以只要把二类边界的通用代码写出来自然边界就是传两个0进去的问题反而最简单。1.3 二类边界为什么叫“弯矩”边界“弯矩”这个词不是随便叫的。在材料力学里一根梁在载荷下的挠曲线如果写成三次样条那么它的二阶导数在物理上正比于截面弯矩。简支梁两端铰接端部弯矩为零所以自然边界条件是有明确物理背景的。如果你手里的问题能拿到端点二阶导比如位移传感器配合加速度计端点处的加速度可以直接测量那么二类边界条件就是更自然的输入接口。相比一类边界给定端点斜率二阶导信息在物理采集端往往更容易获得。而且在端点给错边界值时二类边界对端点附近的影响模式和一类的表现不一样——这个后面在调试部分会专门讲。2. 三弯矩方程二类边界是怎么进入方程组的2.1 三次样条的分段表达式与核心记号设节点为 x0 x1 ... xn子区间个数为n令第i个区间长度为 hi xi - xi-1。记 Mi S(xi)也就是节点xi处的二阶导数值工程里叫弯矩。关键观察S(x)在每个区间上是三次多项式那么S(x)在每个区间上就是一次多项式它在区间端点 [xi-1, xi] 上的值分别是 Mi-1 和 Mi于是可以唯一写成线性插值的形式。对S积分两次再代入两个端点的函数值条件定出积分常数就得到每个区间上的样条表达式S(x) Mi-1(xi-x)^3/(6hi) Mi(x-xi-1)^3/(6hi) (yi-1 - Mi-1hi^2/6)(xi-x)/hi (yi - Mihi^2/6)(x-xi-1)/hi这个式子看着长但结构很清晰前两项是二阶导的贡献后两项是线性插值的修正项。代码里查询点求值用的就是它。再利用S(x)在内部节点xi处左导数等于右导数整理可以得到关于Mi的线性方程组也就是经典的三弯矩方程μi Mi-1 2Mi λi Mi1 di其中i 1, 2, ..., n-1这里各系数的定义是μi hi / (hi hi1)λi 1 - μi hi1 / (hi hi1)di 6 / (hi hi1) * [ (yi1 - yi)/hi1 - (yi - yi-1)/hi ]三弯矩方程是一个三对角方程组未知量是节点处的二阶导 M0, M1, ..., Mn。这个结论对任意单调递增的节点都成立等距非等距都可以代码里没有做任何等距假设。2.2 二类边界条件的关键操作移项三弯矩方程一共有n-1个方程i 1 到 n-1但未知数是n1个M值缺的正好就是边界条件。二类边界直接给出 M0 f(x0)Mn f(xn)所以未知数只剩下内部的 M1 到 Mn-1一共n-1个方程数等于未知数方程组闭合。但注意i1 的方程左边含有 M0in-1 的方程左边含有 Mn。既然 M0 和 Mn 已知把这两项移到等号右边第1个方程μ1 M0 2M1 λ1 M2 d1变成 2M1 λ1 M2 d1 - μ1 M0第n-1个方程μn-1 Mn-2 2Mn-1 λn-1 Mn dn-1变成 μn-1 Mn-2 2Mn-1 dn-1 - λn-1 Mn这个“移项”看起来简单但实际写代码时非常容易出错。很多人一开始会习惯性地套第一类边界的程序框架把 M0 和 Mn 当成未知数结果矩阵维度对不上或者硬凑成 (n1) 阶矩阵后边界行不知道填什么。记住二类边界下矩阵是 (n-1) 阶边界项全部以减法的形式进入右端向量这是三弯矩方程在二类边界下的最大区别。2.3 右端项 di 的几何直觉右端项di的形式接近“两个相邻区间割线斜率之差再乘以加权因子”它衡量的是第i个节点附近曲线的弯曲程度。如果曲线在某节点附近越弯两个相邻区间的一阶差商变化越大di就越大需要更大的弯矩来维持曲率连续。反过来如果数据点本身落在一条直线上所有di都是0解出来的弯矩全是0三次样条自动退化成直线——这是个很有用的自检特性。等距节点时公式会退化成 di 3/h^2 * (yi1 - 2yi yi-1)一眼就能看出它和二阶差商的关系。非等距时用原式更稳妥直接套等距公式会引入误差。这一点在后面的常见问题里还会提到。3. 完整MATLAB代码与逐段解读3.1 函数接口设计与整体结构下面给出完整函数输入节点、端点二阶导、查询点输出插值结果和所有节点弯矩。设计上我倾向返回M因为查看弯矩向量是调试的好手段很多奇怪现象看M一眼就能定位。function [xq_out, yq_out, M] cubicSplineBC2(x, y, M0, Mn, xq) % cubicSplineBC2 - 三次样条插值二类边界条件 % 二类边界给定两端点的二阶导数值 f(x0)M0f(xn)Mn % % 输入: % x - 插值节点横坐标严格单调递增长度 n1 % y - 插值节点纵坐标长度 n1 % M0 - 左端点二阶导数 % Mn - 右端点二阶导数 % xq - 查询点横坐标向量即可内部会排序 % % 输出: % xq_out - 排序后的查询点列向量 % yq_out - 对应插值结果列向量 % M - 所有节点处的二阶导数值长度 n1 % ---- 输入预处理 ---- x x(:); y y(:); xq xq(:); n length(x) - 1; % 子区间数量 if n 1 error(至少需要2个节点才能插值); end h diff(x); % 各区间长度 if any(h 0) error(节点必须严格单调递增); end if length(y) ~ n 1 error(x与y的长度必须一致); end % ---- 组装三弯矩方程组只针对内部节点 ---- % 未知数为 M(2), M(3), ..., M(n)共 n-1 个 % 三弯矩方程: mu_i*M(i-1) 2*M(i) lam_i*M(i1) d_i % 这里 i 是MATLAB下标范围 2 到 n A zeros(n - 1, n - 1); b zeros(n - 1, 1); for i 2:n hL h(i - 1); % 左区间长度 hR h(i); % 右区间长度 dL (y(i) - y(i - 1)) / hL; dR (y(i 1) - y(i)) / hR; mu hL / (hL hR); lam 1 - mu; dd 6 * (dR - dL) / (hL hR); row i - 1; % 矩阵行号 A(row, row) 2; if i 2 A(row, row - 1) mu; end if i n A(row, row 1) lam; end b(row) dd; % 二类边界移项M0、Mn 是已知量移到右端 if i 2 b(row) b(row) - mu * M0; end if i n b(row) b(row) - lam * Mn; end end % ---- 求解内部节点弯矩并组装完整M向量 ---- M_inner A \ b; M zeros(n 1, 1); M(1) M0; M(n 1) Mn; M(2:n) M_inner; % ---- 对查询点逐点求值 ---- xq_out sort(xq); yq_out zeros(size(xq_out)); for k 1:numel(xq_out) xx xq_out(k); if xx x(1) || xx x(end) error(查询点 %g 超出插值区间 [%g, %g], xx, x(1), x(end)); end % 定位查询点所在区间: x(idx) xx x(idx1) idx find(xx x, 1, last); if idx n 1 idx n; % xx 恰好等于右端点 end hi h(idx); xL x(idx); xR x(idx 1); ML M(idx); MR M(idx 1); yL y(idx); yR y(idx 1); % 三次样条分段求值公式 yq_out(k) ... ML * (xR - xx)^3 / (6 * hi) ... MR * (xx - xL)^3 / (6 * hi) ... (yL - ML * hi^2 / 6) * (xR - xx) / hi ... (yR - MR * hi^2 / 6) * (xx - xL) / hi; end end3.2 矩阵组装部分的关键细节组装循环里最需要理解的是下标映射。节点在MATLAB里是1到n1我在循环里用 i 代表内部节点的下标范围2到n。对应的三弯矩方程左边是 M(i-1)、M(i)、M(i1)而未知的只有 M(2) 到 M(n)所以每个方程在矩阵 A 里的行号正好是 row i - 1。列位置对应关系M(i-1) 是第 row-1 列M(i1) 是第 row1 列。因此 if i 2 时给 A(row, row-1) 赋 muif i n 时给 A(row, row1) 赋 lam。这两个条件同时也在防止访问矩阵越界本质上是处理三对角矩阵首尾两行只有一个相邻非零元的特殊情况。移项操作只涉及第一行和最后一行。第一行i2对应的方程本来有 μ2 M0 这一项M0已知所以 b(row) 要减去 μ * M0最后一行in对应的方程本来有 λn Mn 这一项Mn已知所以 b(row) 要减去 λ * Mn。中间行不涉及已知边界弯矩不用动。这个逻辑写错一个符号整个插值结果就歪了但没有任何报错提示只能靠数值自检发现。3.3 查询点求值部分为什么这么写查询点求值的核心是分段定位先找到查询点属于哪个区间再套该区间的三次多项式。我用 find(xx x, 1, last) 做定位逻辑直观适合教学和中等规模数据。数据量很大时这个 find 在循环里是 O(n) 的开销可以换成 MATLAB 自带的 discretize 函数一行替代idx discretize(xx, x);discretize 返回的是 xx 落在 [x(idx), x(idx1)) 半开区间内的下标但要注意它对右端点 x(end) 的处理和 find 略有差异实际使用时要拿边界值测一下。求值公式就是 2.1 节里那个分段表达式把 ML、MR、hi 代进去直接算。还有一个细节值得提当查询点正好落在某个节点 xj 上时定位区间是 [xj, xj1]代入公式后第一项和第三项中 (xR - xx) 是 hi、(xx - xL) 是0最终结果应该精确等于 yj。如果算出来不等于 yj代码的某个环节可能出错了这个可以作为调试时的自检点。3.4 示例脚本正弦函数插值验证下面给一个可直接运行的示例拿 sin(x) 做真值对比插值效果。sin(x) 的二阶导是 -sin(x)所以端点弯矩可以直接精确给出适合当标准测试题。%% 示例对 y sin(x) 做三次样条插值二类边界条件 clear; clc; close all; % 在 [0, 2*pi] 上取7个等距节点 x linspace(0, 2*pi, 7); y sin(x); % 二类边界f(x) -sin(x)端点处精确值 M0 -sin(x(1)); Mn -sin(x(end)); % 生成查询点 xq linspace(0, 2*pi, 300); % 调用三次样条 [xq, yq, M] cubicSplineBC2(x, y, M0, Mn, xq); % 绘图对比 figure(Position, [100 100 850 550]); plot(x, y, ro, MarkerSize, 8, LineWidth, 1.5); hold on; plot(xq, yq, b-, LineWidth, 2); plot(xq, sin(xq), k--, LineWidth, 1.2); legend(插值节点, 三次样条(二类边界), sin(x) 真值, ... Location, northeast); xlabel(x); ylabel(y); grid on; box on; title(三次样条插值——二类边界条件给定端点二阶导); % 最大绝对误差 err max(abs(yq - sin(xq))); fprintf(最大绝对误差 %.3e\n, err); % 打印弯矩向量 fprintf(节点处的二阶导数值弯矩\n); disp([x, M]);7个等距节点、300个查询点跑下来最大绝对误差通常在 1e-3 量级因为节点少、区间大。把节点加到 10 个误差会迅速降到 1e-5 以下。这说明三次样条对光滑函数的逼近精度随节点加密提升很快和整体多项式插值在端点处的“指数爆炸式”误差完全是两种性格。4. 实操验证怎么证明你的代码是对的4.1 万能自检直线与二次函数拿到代码第一件事不是拿真实数据跑而是先用两个解析函数自检。第一个y 2x 1。无论边界条件给什么二阶导都是0所以 M0、Mn 都传0任意查询点的插值结果都该精确等于 2xq1。如果这段有问题说明矩阵组装环节有bug。第二个y x^2。它的二阶导恒为2所以 M0 2Mn 2。插值结果必须精确等于 xq.^2误差应该在浮点舍入级别1e-14 左右而不是 1e-3 这种量级。为什么二次函数能精确被三次样条重建因为三次样条的每一段是三次多项式而二次函数本身就是三次多项式的特例只要边界条件给对求出来的弯矩应该全部等于2代入分段公式后完全恢复原函数。这个测试非常狠一旦代码里移项符号写反、系数写错、矩阵下标映射错位二次函数测试立刻暴露问题误差会达到肉眼可见的程度。我建议所有人都把这个当成提交代码前的必过项。4.2 节点密度对精度的影响用 sin(x) 做真值固定区间 [0, 2π]改变节点数量记录最大绝对误差节点数最大绝对误差52.1e-2 左右78.5e-4 左右98.9e-5 左右131.0e-6 左右25接近浮点精度上限可以看到节点数翻倍误差大致下降三到四个数量级这是三次样条的四阶收敛特性在起作用。实际使用中如果发现加密节点后误差下降不明显说明问题可能不在节点密度而在于边界条件给得不合理或者原始数据本身噪声太大已经超出了插值能处理的范围。4.3 和 MATLAB 内置 spline/interp1 的对比很多同学图省事直接调 interp1(x, y, xq, spline)以为这就是三次样条。其实 MATLAB 内置 spline 函数默认用的是 not-a-knot 边界条件不是自然边界也不是二类边界。not-a-knot 的含义是让最左边两个区间和最右边两个区间在拼接点处的三阶导数也相等相当于把前两个区间拼成一个多项式、后两个区间拼成一个多项式。两者的差异在数据点比较稀疏时会很明显特别是端点附近你的二类边界端点曲率被强制约束为给定值端部曲线服从你给的物理信息。内置 spline端部没有明确曲率约束完全由内部数据推断曲线在端部更“自由”。如果你手里有可靠的端点曲率数据二类边界会比 not-a-knot 在端部更准。如果完全没有边界信息老老实实用自然边界或者 not-a-knot 反而更省事强行编造端点弯矩只会让结果变差。这一点也是很多人误用插值函数的根源之一。5. 调试经验与常见坑速查5.1 边界条件给错时的现象二类边界最典型的错误是 M0、Mn 给得离谱比如把端点二阶导给了10实际接近0插值曲线会在端点附近出现明显的“抬头”或“下压”现象。原因是三次样条的求值公式里弯矩直接控制区间端部的曲率边界弯矩错误的影响会沿着最靠近边界的那一段曲线向外蔓延越靠近原点误差越大。反过来如果 M0、Mn 给得精确但代码里移项符号写反曲线往往不会有明显的局部翘起而是整体偏移——形态上甚至更平滑这就更隐蔽了。所以自检时不能只看图形“像不像”一定要算最大绝对误差并对比弯矩向量 M 是否合理。对二次函数M 应该全是2如果 M 里出现负数或者大幅波动基本可以断定组装环节有错。5.2 查询点越界与右端点问题查询点只要超出 [x(1), x(end)]程序就会报错。这是刻意设计的因为三次样条没有定义区间外的行为强行外推会得到荒谬的数值。右端点处的坑则是另一种xx x(end) 时find(xx x, 1, last) 返回 n1但区间下标最大是 n所以必须把它压回 n。漏掉这一步程序会报“索引超过数组边界”而且只在查询点恰好包含右端点时才触发测试数据稍微少一个点就可能躲过去等到正式跑数据就崩非常烦人。5.3 非等距节点和过密节点的处理三弯矩方程天然支持非等距节点代码里所有系数都用 hi 和 hi1 分别计算没有做等距假设。但非等距节点带来的数值敏感性要留意如果某个区间长度远小于相邻区间μ 或 λ 会接近0或1三对角矩阵虽然还是对角占优但右端项 di 因为除以 (hL hR) 会放大局部噪声。实测下来当相邻区间长度相差超过两个数量级时插值曲线容易出现微小波纹。节点过密是另一个隐患。比如在 [0,1] 上取1000个等距节点h 是 1e-3 量级差商计算中相邻大数相减会损失有效数字M 会出现不规则的微小抖动。这种场景下应该考虑平滑样条CSAPS类工具而不是纯插值三次样条。5.4 常见错误速查表症状可能原因解决办法插值结果整体偏移移项符号写错或边界项未移项检查 i2 和 in 时的 b(row) 处理端点附近明显翘起M0/Mn 给错用二次函数自检确认边界值正确矩阵维度错误把 M0、Mn 当未知数矩阵建成了 (n1) 阶记住二类边界只求内部 n-1 个弯矩查询点包含右端点时报错find 返回 n1 后未修正加 if idx n1, idx n;节点接近时误差反而变大h 太小导致差商舍入误差放大改用更高精度数据或考虑平滑方法与内置 spline 结果对比差异大边界条件类型不同正常现象不要混比确认比较的前提6. 从二类边界出发的扩展路线6.1 一行思路改成第一类边界第一类边界给定端点一阶导 y0 和 yn。此时 M0 和 Mn 也是未知数三弯矩方程需要从 (n-1) 阶扩成 (n1) 阶内部行不变只在首尾补两行边界方程左端方程2M0 M1 (6/h1) * ((y1-y0)/h1 - y0)右端方程Mn-1 2Mn (6/hn) * (yn - (yn-yn-1)/hn)也就是说第一类边界和二类边界最大的区别是二类直接把端点弯矩放进右端项矩阵阶数小一类要把边界方程写进矩阵的首尾行矩阵更大但内部核心逻辑完全一样。理解了二类边界的移项一类边界只是换两行方程的事。6.2 自然边界只是二类边界的特例如果没有任何边界信息可用直接把 M0 0、Mn 0 传给本函数得到的就是自然三次样条。这也是教科书默认的边界条件对应简支梁端部无弯矩的物理场景。实际工程中自然边界是最常用、最稳妥的默认选择二类边界只在你确实掌握端点曲率信息时才比自然边界更优。不要因为没有信息就硬编一个端点弯矩那还不如自然边界。6.3 再往上走追赶法与B样条当节点数 n 超过几百时直接用 A\b 求三对角方程虽然能用但不是最高效的做法。三对角矩阵可以用 Thomas 算法追赶法在 O(n) 时间内求解内存占用也只有 O(n)这对大规模数据很有意义。我给的代码为了保证可读性没有写追赶法但矩阵 A 的三对角结构已经摆在那里照标准算法改就行。再往后如果你需要处理带噪声的数据而不是完全插值就要跳出三次样条的框架去看平滑样条和B样条的世界了。我个人在实际调试中的最大体会是三次样条这类数值方法代码本身并不复杂真正考验人的是对边界条件的理解和对下标映射的细心程度。你花十分钟写完组装循环可能花两个小时找那个看不见的符号错误。所以再次强调先跑二次函数自检再上真实数据这条路径能帮你省下大量排查时间。