做数值分析课程设计或者工程数据预处理的时候我经常看到有人一遇到插值就打开interp1默认用linear而一旦数据结构复杂一点比如要算一阶导数、二阶导数线段之间的折角就会让导数图像非常难看。这时候就该上三次样条了。而在三次样条的三种边界条件里二类边界又是一个特别实用、但参考资料经常讲不透的方向给定两个端点处的二阶导数值其中当两端二阶导都为0时就是常说的自然边界条件。这篇文章我直接把完整可用的Matlab代码贴出来配合三弯矩方程的推导、边界条件处理方式和使用注意事项希望能让第一次接触这个问题的同学少走弯路。这篇文章适合几类人正在写数值分析作业的学生、需要在数据平滑和求导场景里用样条的工程师、想搞懂三弯矩方程而不是只会调工具箱的初学者。我默认你至少会基本的Matlab语法也见过diff、zeros这类函数如果你完全没接触过矩阵求解也没关系追赶法部分我会拆开讲。1. 三次样条插值为什么值得用二类边界又是干什么的1.1 从“一根被数据点钉住的弹性尺子”说起想象你有一块木板上面钉了一排钉子钉子高度和数据点的y值一致。现在拿一根又细又长的弹性木条从钉子之间穿过去它会在每个钉子位置被卡住其余地方自然弯曲成一条连续、光滑的曲线。三次样条就是这根木条的数学版本在相邻两个钉子之间它是一段三次多项式在两个钉子交界处曲线不仅高度一样连斜率和弯曲程度都完全连续。这种二阶导数连续的特性让三次样条在视觉上和物理上都显得非常“自然”。相比之下如果你硬用一个高次多项式去穿过所有钉子数据点一多两端就会剧烈摆动这就是常说的龙格现象如果退而求其次用分段线性插值钉子处会形成一个明显的折角导数不连续后续求速度、算曲率时会变得一塌糊涂。三次样条正好站在中间每一段足够简单整体又足够光滑代价只是需要解一个不大的线性方程组。这也是为什么数值分析课程里三次样条几乎是必讲内容的根本原因。1.2 三类边界条件的分工光有“每段三次多项式”和“内部连续”还不足以唯一确定样条因为整体未知量比方程数多两个必须由端点处的额外条件来补足。工程和教材里最常见的就是下面三类边界类型数学形式典型场景第一类固定边界给定S(x(1))和S(x(n))端点处的斜率/速度已知第二类自由/弯矩边界给定S(x(1))和S(x(n))端点处的弯矩/加速度已知若都为0就是自然边界周期边界S(x(1))S(x(n))S(x(1))S(x(n))闭合曲线、周期信号二类边界的数学表述是给出样条在左右两个端点处的二阶导数值记为M(1)M0M(n)Mn。当M0Mn0时这条样条被称为自然边界三次样条此时样条在端点处“没有任何额外弯曲”类似于一根两端松开的弹性细梁只靠内部数据点撑起来。很多教材默认先讲自然边界因为它最容易理解数值上也最稳定但工程里很多时候端点的二阶导并不是0例如机构启动瞬间末端加速度不为零这时候就必须用一般的二类边界而不是自然边界。1.3 什么时候优先选二类边界第一类边界需要知道端点的一阶导数这在理论题里很好办导数表达式直接给但在实测数据里端点斜率往往要拿差分近似噪声会被放大给出来的边界条件本身就不靠谱。二类边界需要的端点二阶导在很多物理问题里反而有明确含义梁端弯矩、运动末端的角加速度、曲线的曲率约束。比如你处理一段机器人关节角度轨迹起止时刻的角加速度如果已知通常由动力学方程给出那就直接用二类边界。所以我的建议是课程作业里如果明确给了端点二阶导或者题目条件适合用自然边界就放心用二类边界如果题目只给端点一阶导那请老老实实写第一类边界的版本不要硬套本文代码。边界条件选错样条照样能画出来但曲线形态和真实情况可能差得十万八千里。后面第4部分我会专门讲怎么自查。2. 三弯矩方程把光滑条件变成可解的线性方程组2.1 为什么需要“三弯矩”这个中间量如果直接设每一段三次多项式为S_i(x) a_i*x^3 b_i*x^2 c_i*x d_i一段有4个未知数n段就是4n个未知数。约束条件算下来也正好是4n个理论上可以硬解但矩阵会是块结构写起来很痛苦而且数值上不稳定。更聪明的做法是先换一组参数用每个节点处的函数值y(i)和二阶导数值M(i)来表达样条。一旦这样表达函数值连续条件和二阶导连续条件就被“自动焊死”了只剩下“一阶导连续”需要额外满足最后只要求解n个未知数M。这个M就是三弯矩方程里的“弯矩”。名字来自材料力学梁的弯曲程度和弯矩成正比而样条的二阶导数恰好对应弯曲程度。所以M(i)可以理解为第i个数据点处的“数学弯矩”这也是三弯矩方程名称的由来。2.2 从分段三次表达式到内部节点方程在每个区间[x(i), x(i1)]上用h(i)表示区间长度三次样条可以写成S_i(x) M(i)*(x(i1)-x)^3/(6*h(i)) M(i1)*(x-x(i))^3/(6*h(i)) (y(i)-M(i)*h(i)^2/6)*(x(i1)-x)/h(i) (y(i1)-M(i1)*h(i)^2/6)*(x-x(i))/h(i)这个式子看起来长但结构很清楚前两项是二阶导数对三次项的贡献后两项是端点函数值做线性插值后再减去弯矩修正。把xx(i)代进去能立刻得到y(i)求二阶导能得到M(i)所以插值条件和二阶导连续自动满足。接下来需要保证一阶导连续。对上面的表达式求一阶导在内部节点x(i)处让左区间和右区间的斜率相等化简后得到一个三系数递推关系mu(i)*M(i-1) 2*M(i) lambda(i)*M(i1) g(i)其中lambda(i) h(i) / (h(i-1) h(i)) mu(i) h(i-1) / (h(i-1) h(i)) g(i) 6 / (h(i-1) h(i)) * ( (y(i1)-y(i))/h(i) - (y(i)-y(i-1))/h(i-1) )这里的i取2到n-1也就是所有内部节点。注意不同教材里lambda和mu的符号可能互相调换但本质都是“左区间长度占比”和“右区间长度占比”两者之和为1。如果你不想记符号记住一句话就行M(i-1)的系数是左边区间长度占总长度的比例M(i1)的系数是右边区间长度占比中间M(i)的系数固定是2。等距节点时所有h(i)相等系数会简化成0.5、0.5右边g(i)变成3*(y(i1)-2*y(i)y(i-1))/h^2。这个简化形式在很多参考书里出现频率极高但只适用于等距节点后面会提醒大家别在非等距情况瞎套。2.3 二类边界的矩阵化M0和Mn直接进方程组上面的内部节点方程有n-2个未知数M(1)..M(n)有n个差的两个方程正好由二类边界条件补上第一行M(1) M0最后一行M(n) Mn。中间的每一行就是2.2节里的递推关系。于是整体是三对角线性方程组A*M d形如[ 1 0 0 0 ] [ M0 ] [ μ2 2 λ2 0 ] [ g2 ] [ 0 μ3 2 λ3 ] * M [ g3 ] [ 0 0 0 1 ] [ Mn ]这里我随意写了4个节点的示意。实际n个节点矩阵是n x n。主对角线上除了首尾行是1中间全部是2次对角线上的系数加起来等于1。这个矩阵严格对角占优所以是良态的用追赶法或者直接用Matlab的反斜杠都能稳解。注意边界条件写进方程组时是“硬约束”M0和Mn并不是通过内部方程推算出来的而是直接被赋值。可以理解为样条在端点处的弯曲程度由外部物理条件给定而不是由数据自动决定。这也是二类边界和自然边界在实现上的唯一区别自然边界只是把M0和Mn都设成0。2.4 追赶法求解三对角系统三对角线性方程组不需要动用通用高斯消元。追赶法Thomas算法专门利用三对角结构把复杂度从O(n^3)降到O(n)在节点上千时差距非常明显。思路分成两步追从第一行开始用消元法把下对角线消掉同时把主对角线和右端项更新成新的值赶从最后一行开始回代把每个未知数一个一个解出来。追赶法的前提是矩阵对角占优三弯矩方程天然满足这个条件。但程序里还是要留一手一旦发现消元过程中主元的绝对值接近0立刻报错而不是把NaN带进结果。我在下面提供的thomas_solve函数里就加了这么一道防线。如果偷懒也可以直接写M A \ dMatlab会自动选择合适算法不过为了讲清楚原理我这里还是把追赶法完整贴出来毕竟数值分析作业里老师往往就要求手写追赶法。3. 完整Matlab代码三弯矩方程从推导到实战3.1 主函数cubic_spline_2nd带边界值的二类样条function [S, M] cubic_spline_2nd(x, y, M0, Mn, xi) % CUBIC_SPLINE_2ND 三次样条插值第二类边界条件 % 输入: % x 节点横坐标必须严格单调递增行/列向量均可 % y 节点函数值与x长度一致 % M0 左端点二阶导数 S(x(1)) % Mn 右端点二阶导数 S(x(end)) % xi 待插值点标量或向量 % 输出: % S 插值结果列向量若需要与原xi方向一致可自行转置 % M 各节点处的二阶导数值即三弯矩方程的解 validateattributes(x, {double}, {vector, finite}); validateattributes(y, {double}, {vector, finite}); if numel(y) ~ numel(x) error(x和y长度必须一致); end if any(diff(x) 0) error(x 必须严格单调递增); end if numel(x) 3 error(三次样条至少需要 3 个节点); end x x(:); y y(:); n numel(x); h diff(x); % h(i) x(i1)-x(i) % ---------- 组装三弯矩方程组 A*M d ---------- A zeros(n, n); d zeros(n, 1); A(1,1) 1; % 左端点二阶导数条件 d(1) M0; A(n,n) 1; % 右端点二阶导数条件 d(n) Mn; for i 2 : n-1 hiPrev h(i-1); hi h(i); mu hiPrev / (hiPrev hi); lambda hi / (hiPrev hi); A(i, i-1) mu; A(i, i) 2; A(i, i1) lambda; d(i) 6 / (hiPrev hi) * ... ( (y(i1)-y(i)) / hi - (y(i)-y(i-1)) / hiPrev ); end % ---------- 用追赶法解三对角方程组 ---------- M thomas_solve(A, d); % ---------- 对待插值点逐点计算样条函数值 ---------- xi xi(:); S zeros(numel(xi), 1); for k 1 : numel(xi) xx xi(k); if xx x(1) || xx x(end) error(插值点 %.4f 超出数据范围 [%.4f, %.4f], xx, x(1), x(end)); end idx find(x(1:n-1) xx, 1, last); % 定位区间左端点 if isempty(idx), idx 1; end if idx n, idx n - 1; end hh h(idx); xL x(idx); xR x(idx1); yL y(idx); yR y(idx1); ML M(idx); MR M(idx1); S(k) ML*(xR-xx)^3/(6*hh) MR*(xx-xL)^3/(6*hh) ... (yL - ML*hh^2/6)*(xR-xx)/hh ... (yR - MR*hh^2/6)*(xx-xL)/hh; end end代码注释已经把每个模块说清楚了。组装矩阵时内部节点循环从2到n-1和数学推导一节里的i对应边界行直接放在首尾。求解完成后M就是所有节点处的二阶导数值这也是三弯矩方程的名字由来。如果你只需要节点处的M而不需要插值调用[~, M] cubic_spline_2nd(x, y, M0, Mn, [])即可xi给空数组时会得到空S但M照常返回。3.2 追赶法子函数thomas_solve稳定又好调试function X thomas_solve(A, d) % THOMAS_SOLVE 追赶法解三对角方程组 A*X d n numel(d); if n 1 X d / A(1,1); return; end a diag(A); % 主对角线 c diag(A, 1); % 上对角线c(i) A(i,i1) b diag(A, -1); % 下对角线b(i) A(i1,i) beta zeros(n, 1); gamma zeros(n-1, 1); y zeros(n, 1); beta(1) a(1); gamma(1) c(1) / beta(1); y(1) d(1) / beta(1); for i 2 : n-1 beta(i) a(i) - b(i-1) * c(i-1) / beta(i-1); if abs(beta(i)) eps error(追赶法主元接近0矩阵可能非对角占优或数据异常); end gamma(i) c(i) / beta(i); y(i) (d(i) - b(i-1)*y(i-1)) / beta(i); end beta(n) a(n) - b(n-1) * c(n-1) / beta(n-1); if abs(beta(n)) eps error(追赶法主元接近0矩阵可能非对角占优或数据异常); end y(n) (d(n) - b(n-1)*y(n-1)) / beta(n); X zeros(n, 1); X(n) y(n); for i n-1 : -1 : 1 X(i) y(i) - gamma(i) * X(i1); end end这个子函数可以单独拿出来当工具用只要矩阵是三对角的就能解。我在消元过程中加了主元接近0的报错平时可能用不上但一旦节点数据异常它会第一时间告诉你问题出在矩阵上而不是让你面对一堆NaN无从下手。实际测试下来只要x严格递增且数据不含NaN这个函数在几千个节点上运行依旧稳定。3.3 三个拿来即用的调用示例示例1自然边界处理正弦函数。正弦函数在[0, 2*pi]两个端点的真实二阶导数都是0所以自然边界条件恰好就是精确边界插值效果理应非常好。x linspace(0, 2*pi, 9); y sin(x); xi linspace(0, 2*pi, 500); [S, M] cubic_spline_2nd(x, y, 0, 0, xi); plot(x, y, ko, xi, S, b-, xi, sin(xi), r--); legend(样本点, 三次样条, sin(x), Location, best);运行后可以看到三条线几乎重合最大误差大约在1e-5量级。如果节点数增加到17个误差会继续下降约16倍这是三次样条典型的收敛行为。示例2抛物线精确边界自检。取x[0,1,2]y[0,1,4]这组数据来自yx^2。它的真实二阶导数处处为2所以如果边界传M02Mn2样条应当精确重现抛物线。x [0; 1; 2]; y x.^2; xi 0:0.01:2; [S, ~] cubic_spline_2nd(x, y, 2, 2, xi); err max(abs(S - xi.^2)); % 期望约1e-15这个实验跑通说明你的矩阵组装和求解逻辑没有错。我每次改完代码都用这个例子做回归测试比看十遍推导都管用。示例3对带噪声的离散数据做平滑。实际数据处理里x不一定是等距的y也可能带有测量噪声。假设你采集了一段等距时间戳但数值有毛刺的轨迹t (0:10); p [0.00, 0.12, 0.28, 0.35, 0.55, 0.70, 0.93, 1.12, 1.38, 1.65, 2.01]; ti linspace(0, 10, 300); [S, ~] cubic_spline_2nd(t, p, 0, 0, ti); plot(t, p, ko, ti, S, b-);注意样条插值是一种插值方法它会让曲线精确穿过每个样本点所以对噪声大的数据并不做“平滑滤波”而是提供一条光滑的连续曲线方便后续计算速度和加速度。如果你做的是真平滑应该先滤波再插值这个顺序不要搞反。3.4 自然边界一键封装把M0和Mn默认成0二类边界最常用的情况就是自然边界。为了调用方便可以在主函数下面加一个小封装function [S, M] cubic_spline_natural(x, y, xi) [S, M] cubic_spline_2nd(x, y, 0, 0, xi); end如果你希望更灵活一点还可以把主函数的签名改成允许缺省M0和Mnif nargin 4 || isempty(M0) M0 0; end if nargin 5 || isempty(Mn) Mn 0; end这样cubic_spline_2nd(x, y, xi)也会被解释成自然边界。但要注意如果M0传了空、Mn没传Matlab里参数位置会乱套我个人的习惯是保留显式传参宁可多打两个0也不在接口上搞太多隐式逻辑。4. 常见问题与调试实录我踩过的高频坑4.1 下标错位h(i-1)和h(i)别搞反三弯矩方程的下标是整个代码最容易写错的地方。我当时第一次手写时把内部节点方程里的h(i-1)和h(i)写反了结果矩阵凑出来系数和不为1解出来的M剧烈振荡插值曲线在个别区间出现大波浪。排查方法很简单构造一个只有3个节点的抛物线例子手算一遍A矩阵和右端向量再和程序里打出来的A、d逐项对照。A矩阵内部行两个非对角系数之和必须等于1主对角线固定是2这是最直观的检查项。4.2 边界值传反、传成第一类边界第一类边界和第二类边界的区别只在端点条件前者给一阶导后者给二阶导。如果你手上题目给的是f(x0)...却照搬本文代码传M0解出来的样条会把端点“弯曲”成另一个样子数据中间的曲线也会被连带影响。读题时先确认条件类型。另外M0和Mn的位置也要小心传反之后曲线在两端都会出问题这类错误用抛物线自检一秒钟就能暴露出来。4.3 非等距节点不能照抄“0.5系数”等距节点时lambda和mu都等于0.5方程简化为0.5*M(i-1)2*M(i)0.5*M(i1)g(i)。有人图省事把0.5写死一旦x换成非等距结果就废了。务必从h向量动态计算。我在第2节给出的代码里就是按公式现算的所以无论等距还是非等距都能直接用。如果拿不准自己的x是否等距可以直接看diff(x)是否全部相等别靠肉眼判断。4.4 插值点定位与越界外推待插值点位于两个节点之间时要先定位区间。代码里用的是find(x(1:n-1) xx, 1, last)好处是当xx正好等于最后一个节点时不会越界。有些写法喜欢用sum(x xx)这时如果xx等于x(end)得到n再取y(idx1)就必然越界需要额外加一句idx min(idx, n-1)。还有一个原则性问题二类边界的样条只保证在插值区间内有意义区间外部没有任何约束。你把xi取到数据范围之外代码虽然在报错前可以继续算但结果会沿着端点的斜率直线飞出去毫无物理意义。所以我的函数干脆直接报错。如果你确实需要外推请单独设计外推策略不要指望样条自己会“猜”出区间外面的趋势。4.5 三个自检实验跑之前先验炮给一张速查表我把它当成这个函数的出厂测试序号测试数据边界条件期望结果1yx^2x[0,1,2]M02, Mn2插值结果与真实抛物线误差接近机器精度2ysin(x)xlinspace(0,2*pi,9)M00, Mn0最大误差约1e-5量级加密节点后误差下降3任意数据与成熟工具箱里的二阶导边界样条对比相同点处插值结果一致第一个实验验证的是代码逻辑第二个验证的是自然边界的合理性第三个则是和成熟实现交叉验证。我建议你把这三个实验写成脚本改完代码就一键跑一遍每次改完都能快速知道有没有引入回归。最后说点我实际使用中的体会二类边界条件最容易被低估的就是“边界值的准确性”。我曾经处理一组实验轨迹数据图省事把所有情况都默认成自然边界曲线虽然光滑但在端点附近明显被“压下去”了一块和实测趋势对不上。后来把端点二阶导换成由动力学估算出的真实值曲线形态立刻正常。所以拿到一套代码之后先别急着套数据先花十分钟跑一遍上面的抛物线自检再想清楚你的端点二阶导到底应该是什么。样条算法本身不复杂复杂的是边界条件里藏着的物理意义搞懂了这一点三弯矩方程对你来说就再也不是黑箱了。