数值分析复习：样条插值

本篇文章适合个人复习翻阅，不建议新手入门使用

样条插值

1.Euler-Bernoulli梁

我们考虑水平长横梁发生微小形变的问题，主要研究 Euler-Bernoulli 梁，这一环节我们的总目标是证明 Euler-Bernoulli 梁的数学模型为弯矩等于挠度的二阶导

1.1 弹性力学的基本概念

定义： Euler-Bernoulli 梁
现有一个细长横梁，其均匀且各向同性，若横梁在弯曲过程中横截面始终垂直于中性轴（Neutral axis），则称这个横梁为 Euler-Bernoulli 梁

定义：挠度
假设横梁水平并纵向弯曲，横梁弯曲时一点产生的位移称为该点的挠度，记为 $y(x)$

称横梁形变时内部产生的力为内力，可以将内力分解如下：

定义：内力的分解
在平面杆件的任意截面上一般有三个内力分量：轴力 $F_N$，剪力 $F_Q$，弯矩 $M$；其定义如下

1. 截面上应力沿杆轴切线方向的合力，称为轴力 $F_N$，规定轴力以拉为正，以压为负
2. 截面上应力沿杆轴法线方向的合力，称为剪力 $F_Q$，规定剪力以绕微端隔离体顺时针旋转为正，反之为负
3. 截面上应力对截面形心的力矩称为弯矩 $M$，方向垂直于梁的中性轴，大小等于法应力在该截面上的合力乘以截面到中性轴的距离。在水平杆件中，规定当弯矩使杆件下部受拉时，弯矩为正，当弯矩使杆件上部受拉时，弯矩为负
   

接下来定义应力，即单位面积的内力

定义：应力
根据物体连续性的假设，可认为作用在微小面 $\Delta S$ 上的力是连续分布的，内力 $\Delta P$ 是这个分布力的合力，则称
$$F_{\nu }=\underset{\Delta S\to 0}{\lim }\frac{\Delta P}{\Delta S}$$为应力，$\nu$ 表示微分面的法向

定义：应力张量
记应力的各分量如下

1. x 方向上的正向应力 ${\sigma }_{xx}$
2. x 为法向的平面中，y 方向的剪切应力 ${\tau }_{xy}$
3. x 为法向的平面中，z 方向的剪切应力 ${\tau }_{xz}$
4. y 为法向的平面中，x 方向的剪切应力 ${\tau }_{yx}$
5. y 方向上的正向应力 ${\sigma }_{yy}$
6. y 为法向的平面中，z 方向的剪切应力 ${\tau }_{yz}$
7. z 为法向的平面中，x 方向的剪切应力 ${\tau }_{zx}$
8. z 为法向的平面中，y 方向的剪切应力 ${\tau }_{zy}$
9. z 方向上的正向应力 ${\sigma }_{zz}$

这些元素描述了在三个不同的坐标轴方向上的正向应力和剪切应力。称如下的矩阵为应力张量
$$\left(\begin{array}{ccc}{\sigma }_{xx}& {\tau }_{xy}& {\tau }_{xz}\\ {\tau }_{yx}& {\sigma }_{yy}& {\tau }_{yz}\\ {\tau }_{zx}& {\tau }_{zy}& {\sigma }_{zz}\end{array}\right)$$

下面我们用应变描述横梁弯曲时小微元的空间变化量

定义：应变
根据物体连续性的假设，要求变形前连续的物体在变形后仍保持为连续体，物体 $D$ 内一点 $P$ 经变形后移到 $D_1$ 内 $P_1$ 的位移在三个坐标轴上的分量分别用 $u,v,w$ 表示
正应变定义为棱边的伸长变化
$${\epsilon }_x=\frac{\partial u}{\partial x},{\epsilon }_y=\frac{\partial v}{\partial y},{\epsilon }_z=\frac{\partial w}{\partial z}$$剪应变定义为棱边间夹角的变化
$${\gamma }_{xy}=\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}$$$${\gamma }_{xz}=\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}$$$${\gamma }_{yz}=\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}$$

1.2 Euler-Bernoulli梁的应力与应变

我们取横梁中的一小部分，将其原位置和弯曲变形后的位置画图如下

由图看出，$w$ 是挠度，容易求出应变 ${\epsilon }_x$
$$u_1=u-z\frac{\mathrm{d}w}{\mathrm{d}x}$$$${\epsilon }_x=\frac{\partial u_1}{\partial x}=\frac{\mathrm{d}u}{\mathrm{d}x}-z\frac{{\mathrm{d}}^{2}w}{\mathrm{d}{x}^2}$$
由胡克定律，小微元上的应力与应变成正比，比值为杨氏模量 E，可得 $x$ 方向上的应力为
$${\sigma }_{xx}=E\cdot {\epsilon }_x=E\frac{\mathrm{d}u}{\mathrm{d}x}-Ez\frac{{\mathrm{d}}^{2}w}{\mathrm{d}{x}^2}$$

1.3 Euler-Bernoulli梁的弯矩

设小微元 $A$， $z$ 表示 $A$ 中的某点到中性轴的距离，根据弯矩的定义
$$\begin{array}{rl}M(x)& ={\int }_A{\sigma }_{xx}\cdot z\mathrm{d}A\\ & ={\int }_{A}Ez\frac{\mathrm{d}u}{\mathrm{d}x}-E{z}^2\frac{{\mathrm{d}}^{2}w}{\mathrm{d}{x}^2}\mathrm{d}A\\ & =-EI\frac{{\mathrm{d}}^{2}w}{\mathrm{d}{x}^2}\end{array}$$显然 ${\int }_{A}z\mathrm{d}A=0$，令 $I={\int }_A{z}^2\mathrm{d}A$ 称为面惯性矩

综上所述，我们得到结论$$M(x)=-EI{y}^{\prime \prime }(x)$$即在忽略常数的意义下，弯矩等于挠度的二阶导数

2. 样条函数

定义：弯曲应变能
梁在弯曲形变时会产生势能，称为弯曲应变能 $W$，设水平横梁定义在区间 $[a,b]$上，在忽略常数的意义下，其与弯矩的关系如下：
$$W={\int }_a^b{M}^2(x)\mathrm{d}x$$故在Euler-Bernoulli梁的假设下有
$$W={\int }_a^b[y^{\prime \prime }(x){]}^2\mathrm{d}x$$

接下来我们要刻画使得弯曲应变能最小的形变曲线，我们将其定义为三次样条

2.1 泛函极小解和三次样条函数

给定区间 $[a,b]$ 的分划：$$a=x_0<x_1<\cdots <x_n=b$$给定 $y_i,i=0,1,\dots ,n,y_0^{\prime },y_n^{\prime }$
设
Ω y = { f ( x ) : f ∈ C 2 [ a , b ] , f ( x i ) = y i , f ′ ( x 0 ) = y 0 ′ , f ′ ( x n ) = y n ′ } \Omega_y=\{f(x):f\in C^2[a,b],f(x_i)=y_i,f'(x_0)=y_0',f'(x_n)=y_n'\} Ωy​={f(x):f∈C2[a,b],f(xi​)=yi​,f′(x0​)=y0′​,f′(xn​)=yn′​}

Ω 0 = { η ( x ) : η ∈ C 2 [ a , b ] , η ( x i ) = 0 , η ′ ( a ) = η ′ ( b ) = 0 } \Omega_0=\{\eta(x):\eta\in C^2[a,b],\eta(x_i)=0,\eta'(a)=\eta'(b)=0\} Ω0​={η(x):η∈C2[a,b],η(xi​)=0,η′(a)=η′(b)=0}
$$J(f)={\int }_a^b(f^{\prime \prime }(x){)}^2\mathrm{d}x$$考虑如下的泛函极小问题：
求 $S(x)\in {\Omega }_y$ 满足 $J(S)=\underset{f\in {\Omega }_y}{\min }J(f)$
则 $S(x)\in {\Omega }_y$ 是泛函极小问题的解当且仅当
$$\forall \eta (x)\in {\Omega }_0,{\int }_a^b{\eta }^{\prime \prime }(x)S^{\prime \prime }(x)\mathrm{d}x=0$$

证明：（相当经典）

充分性：任取 $f\in {\Omega }_y$ ，令 $\eta (x)=f(x)-S(x)$
$$\begin{array}{rl}& {\int }_a^b(f^{\prime \prime }(x){)}^2\mathrm{d}x\\ =& {\int }_a^b({\eta }^{\prime \prime }(x)+S^{\prime \prime }(x){)}^2\mathrm{d}x\\ =& {\int }_a^b({\eta }^{\prime \prime }(x){)}^2+(S^{\prime \prime }(x){)}^2\mathrm{d}x\\ \ge & {\int }_a^b(S^{\prime \prime }(x){)}^2\mathrm{d}x\end{array}$$

必要性：任取 $\eta \in {\Omega }_0$，则有$$J(S+\alpha \eta )\ge J(S)$$
记 $g(\alpha )=J(S+\alpha \eta )$，则 $\alpha =0$ 是 $g(\alpha )$ 的极小值点
由于$$\begin{array}{rl}& \frac{\mathrm{d}}{\mathrm{d}\alpha }\left({\int }_a^b(S^{\prime \prime }(x)+\alpha {\eta }^{\prime \prime }(x){)}^2\mathrm{d}x\right)\\ =& {\int }_a^b\left(2\alpha ({\eta }^{\prime \prime }(x){)}^2+2S^{\prime \prime }(x){\eta }^{\prime \prime }(x)\right)\mathrm{d}x\end{array}$$令 $\alpha =0$ ，即得 $${\int }_a^b{\eta }^{\prime \prime }(x)S^{\prime \prime }(x)\mathrm{d}x=0$$

2.2 $S(x)$ 的结构

设 $S(x)$ 即为前文所述泛函极小问题的解，且令 $S(x)\in C^4(x_i,x_{i+1}),i=0,1,\dots ,n-1$，在每个 $x_i$ 处有直到四阶的左、右导数
在每个小区间上使用两次分部积分公式，得
$$\begin{array}{rl}& {\int }_a^b{\eta }^{\prime \prime }(x)S^{\prime \prime }(x)\mathrm{d}x\\ =& \sum _{i=0}^{n-1}{\int }_{x_i}^{x_{i+1}}{\eta }^{\prime \prime }(x)S^{\prime \prime }(x)\mathrm{d}x\\ =& \sum _{i=0}^{n-1}{\eta }^{\prime }(x)S^{\prime \prime }(x){|}_{x_i^{+}}^{x_{i+1}^{-}}-{\int }_{x_i}^{x_{i+1}}{\eta }^{\prime }(x)S^{\prime \prime \prime }(x)\mathrm{d}x\\ =& \sum _{i=0}^{n-1}{\eta }^{\prime }(x)S^{\prime \prime }(x){|}_{x_i^{+}}^{x_{i+1}^{-}}-\eta (x)S^{\prime \prime \prime }(x){|}_{x_i^{+}}^{x_{i+1}^{-}}+{\int }_{x_i}^{x_{i+1}}\eta (x)S^{\prime \prime \prime \prime }(x)\mathrm{d}x\end{array}$$由 $S(x)\in C^4(x_i,x_{i+1})$，将上式合并同类项

$$\begin{array}{rl}& {\int }_a^b{\eta }^{\prime \prime }(x)S^{\prime \prime }(x)\mathrm{d}x\\ =& [{\eta }^{\prime }(b)S^{\prime \prime }(b)-{\eta }^{\prime }(a)S^{\prime \prime }(a)+\sum _{i=1}^{n-1}{S}^{\prime \prime }(x_i)({\eta }^{\prime }(x_i^{+})-{\eta }^{\prime }(x_i^{-}))]\\ & -[\eta (b)S^{\prime \prime \prime }(b)-\eta (a)S^{\prime \prime \prime }(a)+\sum _{i=1}^{n-1}{S}^{\prime \prime \prime }(x_i)(\eta (x_i^{+})-\eta (x_i^{-}))]\\ & +\sum _{i=0}^{n-1}{\int }_{x_i}^{x_{i+1}}\eta (x)S^{\prime \prime \prime \prime }(x)\mathrm{d}x\\ =& \sum _{i=0}^{n-1}{\int }_{x_i}^{x_{i+1}}\eta (x)S^{\prime \prime \prime \prime }(x)\mathrm{d}x\\ =& 0\end{array}$$由 $\eta (x)$ 的任意性，可得 $S^{\prime \prime \prime \prime }(x)$ 在每个小区间为 0

结论
$S(x)$ 是具有下列性质的函数

• 在每个小区间上是三次多项式
• 在整个区间上二阶连续可导
• 全体 $S(x)$ 形成有限维空间 $Dim=4n-3(n-1)=n+3$

2.3 三次样条插值

定义：样条函数

• 在 $n+1$ 个插值节点 $x_0,x_1,\dots ,x_n$ 处与给定的函数值相同
• 在每个小区间 $[x_i,x_{i+1}]$ 上是不超过 $k$ 次的多项式
• 在整个区间上 $k-1$ 阶连续可导

称满足这些条件的 $S_k(x)$ 为 $k$ 次样条函数，其中常用的是三次样条函数，也叫做三次样条插值

定义：边界条件
可以发现，我们给出的全体 $S(x)$ 形成的空间是 $n+3$ 维的，但根据样条函数的定义，却发现只有 $n+1$ 个插值条件，所以还需要两个边界处的约束条件以确定样条函数，称这两个条件为边界条件

常见的边界条件

1. 固支边界条件（第一型边界条件）：给出边界处的一阶导数
2. 自然边界条件（第二型边界条件）：给出边界处的二阶导数，且二阶导数为 0
3. 周期型边界条件（第三型边界条件）：分别令左右边界的函数值，一阶导数值，二阶导数值相等，并给出这三个值
4. Not-A-Knot 边界条件：
   $$S_3^{\prime \prime \prime }(x_1^{+})=S_3^{\prime \prime \prime }(x_1^{-}),S_3^{\prime \prime \prime }(x_{n-1}^{+})=S_3^{\prime \prime \prime }(x_{n-1}^{-})$$
5. 横梁固支端：挠度和截面转角为0
   $$S(x_0)=0,S^{\prime }(x_0)=0$$
   
6. 横梁简支端：挠度和弯矩为0
   $$S(x_0)=0,S^{\prime \prime }(x_0)=0$$
   
7. 横梁自由端：弯矩和剪力为0
   $$S^{\prime \prime }(x_0)=0,S^{\prime \prime \prime }(x_0)=0$$

3. 三次样条函数的构造方法

3.1 三转角

思想：
基于节点构造 $S_3(x)$ 的分段三次Hermite插值，求二阶导使其二阶连续可导

命题：三转角构造三次样条的公式
设节点处的一阶导数值 $m_i$，构造三次hermite插值
$$S_3(x)=\sum _{i=0}^n[y_i{h}_i(x)+m_i{\hat{h}}_i(x)]$$其中 $h,\hat{h}$ 为Hermite插值的基函数，则由三转角构造方法，$m_i$ 满足
$$(1-{\lambda }_i)m_{i-1}+2m_i+{\lambda }_i{m}_{i+1}={\mu }_i$$其中
$i=1,2,\dots ,n-1$

${\lambda }_i=\frac{\Delta x_{i-1}}{\Delta x_{i-1}+\Delta x_i}$

${\mu }_i=3((1-{\lambda }_i)f[x_{i-1},x_i]+{\lambda }_{i}f[x_i,x_{i+1}])$

证明思路
写出每个小区间上的 $S_3(x)$ 的形式，求出节点处（不含区间端点）的左右二阶导数值，令其相等并整理即得

命题：三转角结合边界条件的方程

• 三转角固支边界条件：$$m_0=f^{\prime }(x_0),m_n=f^{\prime }(x_n)$$
• 三转角自然边界条件：
  $$2m_0+m_1=3f[x_0,x_1]={\mu }_0$$$$m_{n-1}+2m_n=3f[x_{n-1},x_n]={\mu }_n$$

3.2 三弯矩

思想：
基于节点构造 $S^{\prime \prime }(x)$ 为连续的分段线性函数，积分两次后代入插值条件 $S_3(x_i)=y_i$，对得到的 $S_3(x)$ 进行求导，使其一阶导数也连续

命题：三弯矩构造三次样条的公式
设 $S_3^{\prime \prime }(x_i)=M_i$，基于节点构造分段线性插值
$$S_3^{\prime \prime }(x)=M_{i-1}\frac{x_i-x}{\Delta x_{i-1}}+M_i\frac{x-x_{i-1}}{\Delta x_{i-1}}$$由三弯矩构造方法，$M_i$ 满足
$${\lambda }_i{M}_{i-1}+2M_i+(1-{\lambda }_i)M_{i+1}=d_i$$

其中
$i=1,2,\dots ,n$

$d_i=6f[x_{i-1},x_i,x_{i+1}]$

证明思路
写出每个小区间上的 $S_3(x)$ 的形式，代入插值条件求解，再求出节点处（不含区间端点）的左右一阶导数值，令其相等并整理即得

命题：三弯矩结合边界条件的方程

• 三弯矩固支边界条件：
  $$2M_0+M_1=d_0=6f[x_0,x_0,x_1]$$$$M_{n-1}+2M_n=d_n=6f[x_{n-1},x_n,x_n]$$

3.3 待定系数法

思想
基于三次样条是多项式的事实，用待定系数法设分段函数在每个小区间上是不同的三次多项式，使其分别满足函数值连续，一阶导数值连续，二阶导数值连续

命题：待定系数法构造三次样条的公式
基于 $n$ 个节点 $(x_i,y_i),i=1,2,\dots ,n$ 构造如下的三次多项式组 $S_i(x),x\in [x_i,x_{i+1}],i=1,2,\dots ,n-1$
$$\{\begin{array}{l}{S}_1(x)=y_1+b_1(x-x_1)+c_1(x-x_1{)}^2+d_1(x-x_1{)}^3\\ S_2(x)=y_2+b_2(x-x_2)+c_2(x-x_2{)}^2+d_2(x-x_2{)}^3\\ ⋮\\ S_{n-1}(x)=y_{n-1}+b_{n-1}(x-x_{n-1})+c_{n-1}(x-x_{n-1}{)}^2+d_{n-1}(x-x_{n-1}{)}^3\end{array}$$使其满足
$$\{\begin{array}{ll}(1):S_i(x_i)=y_i,S_i(x_{i+1})=y_{i+1},& i=1,\dots ,n-1\\ (2):S_i^{\prime }(x_{i+1})=S_{i+1}^{\prime }(x_{i+1}),& i=1,2,\dots ,n-2\\ (3):S_i^{\prime \prime }(x_{i+1})=S_{i+1}^{\prime \prime }(x_{i+1}),& i=1,2,\dots ,n-2\end{array}$$则待定系数 $c_i$ 满足如下方程
$${\delta }_i{c}_i+2({\delta }_i+{\delta }_{i+1})c_{i+1}+{\delta }_{i+1}{c}_{i+2}=3(\frac{{\Delta }_{i+1}}{{\delta }_{i+1}}-\frac{{\Delta }_i}{{\delta }_i})$$其中
$i=1,2,\dots ,n-3$
${\delta }_i=x_{i+1}-x_i,{\Delta }_i=y_{i+1}-y_i$
从而可解 $b_i,d_i(i=1,\dots ,n-2)$
$$b_i=\frac{{\Delta }_i}{{\delta }_i}-\frac{{\delta }_i}{3}(2C_i+C_{i+1})$$$$d_i=\frac{C_{i+1}-C_i}{3{\delta }_i}$$

证明思路
从多项式满足的性质（3）可解出 $d_i$，从性质（1）的后半句（前半句已经成立）解出 $b_i$，将 $b_i,d_i$ 代入性质（2）解得 $c_i$

可以发现以上待定系数仅有 $b_{n-1},d_{n-1}$ 还未解出，它们只需额外依赖 $c_n$ 的值，故我们设 $c_n=\frac{S_{n-1}^{\prime \prime }(x_n)}{2}$

综上所述，有以下命题
命题
在前文的待定系数方程组中，待定系数 $a_i,b_i,c_i,i=1,2,\dots ,n-1$ 由以下方程给出
$${\delta }_i{c}_i+2({\delta }_i+{\delta }_{i+1})c_{i+1}+{\delta }_{i+1}{c}_{i+2}=3(\frac{{\Delta }_{i+1}}{{\delta }_{i+1}}-\frac{{\Delta }_i}{{\delta }_i})$$

对 $i=1,2,\dots ,n-2$ 成立
$$b_i=\frac{{\Delta }_i}{{\delta }_i}-\frac{{\delta }_i}{3}(2C_i+C_{i+1})$$

和
$$d_i=\frac{C_{i+1}-C_i}{3{\delta }_i}$$

对 $i=1,2,\dots ,n-1$ 成立

命题：待定系数法结合边界条件的方程

1. 待定系数法固支边界条件：$b_1=y_1^{\prime },b_n=y_n^{\prime }$
2. 待定系数法自然边界条件：$c_1=0,c_n=0$
   它给出了关于 $c_1,\dots ,c_n$ 的三对角线性方程组
   $$\left(\begin{array}{cccccc}1& 0& & & & \\ {\delta }_1& 2({\delta }_1+{\delta }_2)& {\delta }_2& & & \\ & \ddots & \ddots & \ddots & & \\ & & {\delta }_{n-2}& 2({\delta }_{n-2}+{\delta }_{n-1})& {\delta }_{n-1}& \\ & & & 0& 1& \end{array}\right)\left(\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_{n-1}\\ c_n\end{array}\right)=\left(\begin{array}{c}0\\ 3(\frac{{\Delta }_2}{{\delta }_2}-\frac{{\Delta }_1}{{\delta }_1})\\ ⋮\\ 3(\frac{{\Delta }_{n-1}}{{\delta }_{n-1}}-\frac{{\Delta }_{n-2}}{{\delta }_{n-2}})\\ 0\end{array}\right)$$
3. Not A Knot边界条件：$d_1=d_2,d_{n-1}=d_n$

4. 三次样条的截断误差

命题：三次样条的截断误差
设被插函数 $f(x)\in C^4[a,b]$，$S_3(x)$ 为相应的一或二型插值样条，则有
$$||S_3^{(k)}-f^{(k)}|{|}_{\infty }\le C_k||f^{(4)}|{|}_{\infty }{h}^{4-k},k=0,1,2,\dots$$其中 $h=\underset{i}{\max }\Delta x_i,C_0=\frac{24}{384}(\frac{5}{384}),C_1=\frac{1}{2}(\frac{1}{24}),C_2=\frac{1}{2}(\frac{3}{8})$

证明略

5.代码示例

示例1：待定系数法计算自然三次样条

首先展示使用matlab进行待定系数法求解自然三次样条，创建一个m文件，输入以下代码以定义函数。注意：这里的节点下标是从 1 到 n 的

%待定系数法计算三次样条
%算法思路，通过解三对角方程组解出c_i，进而解出b_i和d_i
%输入：数据点的x,y向量
%输出：系数矩阵b1,c1,d1;b2,c2,d2;...
function coeff=splinecoeff(x,y)
n=length(x);
A=zeros(n,n);%初始化矩阵A,r
r=zeros(n,1);
for i=1:n-1      %定义Delta_i(记为dy(i))和delta_i(记为dx(i))
    dx(i)=x(i+1)-x(i);dy(i)=y(i+1)-y(i);
end
for i=2:n-1
    A(i,i-1:i+1)=[dx(i-1) 2*(dx(i-1)+dx(i)) dx(i)];
    r(i)=3*(dy(i)/dx(i)-dy(i-1)/dx(i-1));
end
%根据需要编写边界条件
A(1,1)=1;      %自然样条，默认r_1=r_n=0（初始化r时赋的值）
A(n,n)=1;

coeff=zeros(n,3);  %初始化解矩阵
coeff(:,2)=A\r;    %求解Ac=r,得系数c
for i=1:n-1        %求解b和d
    coeff(i,3)=(coeff(i+1,2)-coeff(i,2))/(3*dx(i));
    coeff(i,1)=dy(i)/dx(i)-dx(i)*(2*coeff(i,2)+coeff(i+1,2))/3;
end
coeff=coeff(1:n-1,1:3);

接下来我们便可调用这个函数了，在相同文件夹内创建一个新的 m 文件，输入以下代码

clc,clear all
x=[0,1,2];
y=[3,-2,1];
splinecoeff(x,y)

输出结果如下

接下来我们尝试使用plot方法画出图像，继续调用刚刚写好的 splinecoeff 函数，可以得到样条函数的系数，在每个小区间上指定 $k$ 个等距点，代入样条函数求出函数值，依次用线段连接各点，即可得到样条函数的近似图像。

先写好画样条函数图像的代码，以便后续调用

%画三次样条图像（待定系数法）
%输入：x,y数据点向量，每个小区间内指定等距点个数k
%输出：所有点的横坐标向量x1,纵坐标向量y1,以及连接这些点的图像
function [x1,y1]=splineplot(x,y,k)
n=length(x);
coeff=splinecoeff(x,y);   %计算三次样条插值的系数
x1=[];y1=[];              %初始化用于存储插值结果的空数组 x1 和 y1
for i=1:n-1
    xs=linspace(x(i),x(i+1),k);%使用 linspace 函数生成 k 个均匀分布的插值点 xs
    dx=xs-x(i);
    ys=coeff(i,3)*dx;      %使用嵌套乘法（秦九韶算法）求多项式的值
    ys=(ys+coeff(i,2)).*dx
    ys=(ys+coeff(i,1)).*dx+y(i);
    x1=[x1;xs(1:k-1)'];y1=[y1;ys(1:k-1)'];%将每次循环新生成的结果追加到原来的结果上
end
x1=[x1;x(end)];y1=[y1;y(end)];  %将最后一个节点补充上
plot(x,y,'o',x1,y1)

创建一个新的 m 文件进行画图，首先是 $k=3$ 的情形

clc,clear all
x=[0,1,2,3,4,5];
y=[3,1,4,1,2,0];
splineplot(x,y,3)

图像显示如下

然后是 $k=15$ 的情形

clc,clear all
x=[0,1,2,3,4,5];
y=[3,1,4,1,2,0];
splineplot(x,y,15)

图像显示如下

示例2：matlab内置三次样条插值函数spline

具体可参见 Matlab文档：函数-插值-spline
spline函数默认自然边界条件，代码示例如下

clc,clear all
X = [0 1 3 4 6 7];
Y = sin(X);
x = linspace(0,7,200);
plot(x,ppval(spline(X,Y),x))

结果如下

若要将边界条件改为固支型（即规定一阶导数），只需在Y左右两侧添加对应导数值，示例如下

clc,clear all
X = [0 1 3 4];
Y = [0 0 2 2];
x = linspace(0,4,200);
plot(x,ppval(spline(X,[0 Y 0]),x))  %指定左右端点一阶导数均为0

结果如下

参考书籍：《数值分析》李庆扬 王能超 易大义 编