系列文章目录
数值分析——拉格朗日插值
数值分析——牛顿插值多项式
数值分析——埃尔米特（Hermit）插值
数值分析——分段低次插值
文章目录
- 系列文章目录
- 前言
- 一、理论推导
- 1.三次样条函数
- 2.三次样条插值函数的求解条件
- 3.三次样条插值函数的建立
- 二、MATLAB实现
- 1.一阶导数边界条件
- 2.二阶导数边界条件
- 3.测试程序
- 总结
- 参考文献
前言
在之前的博文中接受少了分段低次插值函数，该插值函数具有一致的收敛性，但是光滑性较差，对于像高速飞机的机翼形线，船体放样等型值线往往要求有二阶光滑度，即二阶导数连续。样条曲线是由分段三次曲线并接而成，在连接点上二阶导数连续。本文将介绍以此为理念设计的三次样条插值函数。
本文对一些复杂的流程与概念进行了简化与省略，如果想详细了解的小伙伴可以参考李庆扬老师的《数值分析 -第5版》的第2章（需要电子版请私戳）
一、理论推导
1.三次样条函数
设计一个函数
S
(
x
)
S\left( x \right)
S(x)属于在区间
[
a
,
b
]
[a,b]
[a,b]上二阶导数连续的函数空间，即
S
(
x
)
∈
C
2
[
a
,
b
]
S\left( x \right) \in C^2 \left[ a,b \right]
S(x)∈C2[a,b]。该函数在每个小区间
[
x
j
,
x
j
+
1
]
\left[ x_{j},x_{j+1} \right]
[xj,xj+1]上是三次多项式，其中
a
=
x
0
<
x
1
<
⋯
<
x
n
=
b
a=x_0<x_1<\cdots<x_n=b
a=x0<x1<⋯<xn=b，则称该
S
(
x
)
S\left( x \right)
S(x)是节点
x
0
,
x
1
,
⋯
,
x
n
−
1
,
x
n
x_0,x_1,\cdots,x_{n-1},x_{n}
x0,x1,⋯,xn−1,xn上的三次样条函数。如果在插值点上的函数值
f
0
,
f
1
,
f
2
,
⋯
,
f
n
−
1
,
f
n
f_0,f_1,f_2,\cdots,f_{n-1},f_n
f0,f1,f2,⋯,fn−1,fn是已知的，且：
S
(
x
j
)
=
f
(
x
j
)
=
y
j
j
=
0
,
1
,
⋯
,
n
−
1
,
n
(1)
S\left( x_j \right)=f\left( x_j \right)=y_j \quad j=0,1,\cdots,n-1,n \tag{1}
S(xj)=f(xj)=yjj=0,1,⋯,n−1,n(1)则称函数
S
(
x
)
S\left( x \right)
S(x)为三次样条插值函数。
2.三次样条插值函数的求解条件
由三次样条插值函数的定义可知，在每个区间上函数有四个未知量，共有n个小区间，故共需要求4n个未知量。现已知：
- 函数在区间分割点上的二阶导数连续：
S ( x j − 0 ) = S ( x j + 0 ) 、 S ′ ( x j − 0 ) = S ′ ( x j + 0 ) 、 S ′ ′ ( x j − 0 ) = S ′ ′ ( x j + 0 ) (2) S\left( x_j-0 \right)=S\left( x_j+0 \right)、S^{'}\left( x_j-0 \right)=S^{'}\left( x_j+0 \right)、S^{''}\left( x_j-0 \right)=S^{''}\left( x_j+0 \right) \tag{2} S(xj−0)=S(xj+0)、S′(xj−0)=S′(xj+0)、S′′(xj−0)=S′′(xj+0)(2) - 插值函数在插值点上与原函数相等：
S ( x j ) = f ( x j ) (3) S\left( x_j\right)=f\left( x_j\right) \tag{3} S(xj)=f(xj)(3)
上述总共给出了 3 ( n − 1 ) + ( n + 1 ) = 4 n − 2 3\left(n-1\right)+\left(n+1\right)=4n-2 3(n−1)+(n+1)=4n−2个条件，还需要两个边界条件：
- 已知两端的一阶导数值：
S ′ ( x 0 ) = f ′ ( x 0 ) 、 S ′ ( x n ) = f ′ ( x n ) (4) S^{'}\left( x_0\right)=f^{'}\left( x_0\right)、S^{'}\left( x_n\right)=f^{'}\left( x_n\right) \tag{4} S′(x0)=f′(x0)、S′(xn)=f′(xn)(4) - 已知两端的二阶导数值：
S ′ ′ ( x 0 ) = f ′ ′ ( x 0 ) 、 S ′ ′ ( x n ) = f ′ ′ ( x n ) (5) S^{''}\left( x_0\right)=f^{''}\left( x_0\right)、S^{''}\left( x_n\right)=f^{''}\left( x_n\right )\tag{5} S′′(x0)=f′′(x0)、S′′(xn)=f′′(xn)(5) - 原函数为以
[
x
0
,
x
n
]
\left[x_0,x_n\right]
[x0,xn]为周期的周期函数：
S ( x 0 + 0 ) = S ( x n − 0 ) 、 S ′ ( x 0 + 0 ) = S ′ ( x n − 0 ) 、 S ′ ′ ( x 0 + 0 ) = S ′ ′ ( x n − 0 ) (6) S\left( x_0+0 \right)=S\left( x_n-0 \right)、S^{'}\left( x_0+0 \right)=S^{'}\left( x_n-0 \right)、S^{''}\left( x_0+0 \right)=S^{''}\left( x_n-0 \right) \tag{6} S(x0+0)=S(xn−0)、S′(x0+0)=S′(xn−0)、S′′(x0+0)=S′′(xn−0)(6)如此确定的样条函数被称为周期样条函数。
3.三次样条插值函数的建立
满足上述条件的插值函数可以有多种形式，这里介绍只介绍一种由二阶导数进行推导得出的插值函数。首先由于插值函数
S
(
x
)
S\left( x \right)
S(x)在区间
[
x
j
,
x
j
+
1
]
\left[ x_j,x_{j+1} \right]
[xj,xj+1]中是三次多项式 ，可知二阶导数
S
′
′
(
x
)
S^{''}\left( x \right)
S′′(x)在区间
[
x
j
,
x
j
+
1
]
\left[ x_j,x_{j+1} \right]
[xj,xj+1]中是一次函数，故可设：
S
′
′
(
x
)
=
M
j
x
j
+
1
−
x
h
j
+
M
j
+
1
x
−
x
j
h
j
(7)
S^{''}\left( x \right)=M_j\frac{x_{j+1}-x}{h_j}+M_{j+1}\frac{x-x_j}{h_j} \tag{7}
S′′(x)=Mjhjxj+1−x+Mj+1hjx−xj(7)对上述公式进行积分：
S
(
x
)
=
M
j
(
x
j
+
1
−
x
)
3
6
h
j
+
M
j
+
1
(
x
−
x
j
)
3
6
h
j
+
C
1
x
+
C
2
(8)
S\left( x \right)=M_j\frac{\left(x_{j+1}-x\right)^3}{6h_j}+M_{j+1}\frac{\left(x-x_{j}\right)^3}{6h_j}+C_1x+C_2 \tag{8}
S(x)=Mj6hj(xj+1−x)3+Mj+16hj(x−xj)3+C1x+C2(8)其中
C
1
、
C
2
C_1、C_2
C1、C2分别是两个待定的常数。根据
S
(
x
j
)
=
y
j
S\left( x_j \right)=y_j
S(xj)=yj、
S
(
x
j
+
1
)
=
y
j
+
1
S\left( x_{j+1} \right)=y_{j+1}
S(xj+1)=yj+1可求得：
S
(
x
)
=
M
j
(
x
j
+
1
−
x
)
3
6
h
j
+
M
j
+
1
(
x
−
x
j
)
3
6
h
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
(9)
\begin{align} S\left( x \right)=&M_j\frac{\left(x_{j+1}-x\right)^3}{6h_j}+M_{j+1}\frac{\left(x-x_{j}\right)^3}{6h_j} \\ &+\left( y_j-\frac{M_jh^2_j}{6} \right)\frac{x_{j+1}-x}{h_j}+\left( y_{j+1}-\frac{M_{j+1}h^2_j}{6} \right)\frac{x-x_{j}}{h_j} \end{align} \tag{9}
S(x)=Mj6hj(xj+1−x)3+Mj+16hj(x−xj)3+(yj−6Mjhj2)hjxj+1−x+(yj+1−6Mj+1hj2)hjx−xj(9)其中
M
j
(
j
=
0
,
1
,
⋯
,
n
−
1
,
n
)
M_j\left( j=0,1,\cdots,n-1,n \right)
Mj(j=0,1,⋯,n−1,n)是未知的。对
S
(
x
)
S\left( x \right)
S(x)求导可得：
S
′
(
x
)
=
−
M
j
(
x
j
+
1
−
x
)
2
2
h
j
+
M
j
+
1
(
x
−
x
j
)
2
2
h
j
+
y
j
+
1
−
y
j
h
j
−
M
j
+
1
−
M
j
6
h
j
(10)
S^{'}\left( x \right)=-M_j\frac{\left( x_{j+1}-x \right)^2}{2h_j}+M_{j+1}\frac{\left( x-x_j \right)^2}{2h_j}+\frac{y_{j+1}-y_j}{h_j}-\frac{M_{j+1}-M_j}{6}h_j \tag{10}
S′(x)=−Mj2hj(xj+1−x)2+Mj+12hj(x−xj)2+hjyj+1−yj−6Mj+1−Mjhj(10)由
S
′
(
x
j
−
0
)
=
S
′
(
x
j
+
0
)
S^{'}\left( x_j-0 \right)=S^{'}\left( x_j+0 \right)
S′(xj−0)=S′(xj+0)可得：
μ
j
M
j
−
1
+
2
M
j
+
λ
j
M
j
+
1
=
d
j
,
j
=
1
,
2
,
⋯
,
n
−
2
,
n
−
1
(11)
\mu_jM_{j-1}+2M_j+\lambda_jM_{j+1}=d_j,\quad j=1,2,\cdots,n-2,n-1 \tag{11}
μjMj−1+2Mj+λjMj+1=dj,j=1,2,⋯,n−2,n−1(11)其中
μ
j
=
h
j
−
1
h
j
−
1
+
h
j
λ
j
=
h
j
h
j
−
1
+
h
j
d
j
=
6
f
[
x
j
,
x
j
+
1
]
−
f
[
x
j
−
1
,
x
j
]
h
j
−
1
+
h
j
=
6
f
[
x
j
−
1
,
x
j
,
x
j
+
1
]
j
=
1
,
2
,
⋯
,
n
−
2
,
n
−
1
(12)
\begin{align} \mu_j &=\frac{h_{j-1}}{h_{j-1}+h_{j}} \\ \lambda_j &=\frac{h_{j}}{h_{j-1}+h_{j}} \\ d_j &=6\frac{f\left[ x_j,x_{j+1} \right]-f\left[ x_{j-1},x_{j} \right]}{h_{j-1}+h_{j}}=6f\left[ x_{j-1},x_j,x_{j+1} \right] \\ j &=1,2,\cdots,n-2,n-1 \end{align} \tag{12}
μjλjdjj=hj−1+hjhj−1=hj−1+hjhj=6hj−1+hjf[xj,xj+1]−f[xj−1,xj]=6f[xj−1,xj,xj+1]=1,2,⋯,n−2,n−1(12)针对边界条件公式(4)可得：
[
2
λ
0
μ
1
2
λ
1
⋯
⋯
⋯
μ
n
−
1
2
λ
n
−
1
μ
n
2
]
[
M
0
M
1
:
M
n
−
1
M
n
]
=
[
d
0
d
1
:
d
n
−
1
d
n
]
(13)
\begin{bmatrix} 2 & \lambda_0 & & & \\ \mu_1 & 2 & \lambda_1 & & \\ & \cdots & \cdots & \cdots & \\ & & \mu_{n-1} & 2 & \lambda_{n-1}\\ & & & \mu_n & 2\\ \end{bmatrix} \begin{bmatrix} M_0\\ M_1\\ \colon\\ M_{n-1}\\ M_n\\ \end{bmatrix}= \begin{bmatrix} d_0\\ d_1\\ \colon\\ d_{n-1}\\ d_n\\ \end{bmatrix} \tag{13}
2μ1λ02⋯λ1⋯μn−1⋯2μnλn−12
M0M1:Mn−1Mn
=
d0d1:dn−1dn
(13)其中
λ
0
=
1
d
0
=
6
h
0
(
f
[
x
0
,
x
1
]
−
f
0
′
)
μ
n
=
1
d
n
=
6
h
n
−
1
(
f
n
′
−
f
[
x
n
−
1
,
x
n
]
)
(14)
\begin{align} \lambda_0 &=1 \\ d_0 &=\frac{6}{h_0}\left( f\left[ x_0,x_1 \right]-f^{'}_0 \right) \\ \mu_n &=1 \\ d_n &=\frac{6}{h_{n-1}}\left( f^{'}_n-f\left[ x_{n-1},x_n \right] \right) \end{align} \tag{14}
λ0d0μndn=1=h06(f[x0,x1]−f0′)=1=hn−16(fn′−f[xn−1,xn])(14)针对边界条件公式(5)可得：
[
2
λ
1
μ
1
μ
2
2
λ
2
⋯
⋯
⋯
μ
n
−
1
2
λ
n
−
1
λ
n
μ
n
2
]
[
M
1
M
2
:
M
n
−
1
M
n
]
=
[
d
1
d
2
:
d
n
−
1
d
n
]
(15)
\begin{bmatrix} 2 & \lambda_1 & & & \mu_1\\ \mu_2 & 2 & \lambda_2 & & \\ & \cdots & \cdots & \cdots & \\ & & \mu_{n-1} & 2 & \lambda_{n-1}\\ \lambda_n & & & \mu_n & 2\\ \end{bmatrix} \begin{bmatrix} M_1\\ M_2\\ \colon\\ M_{n-1}\\ M_n\\ \end{bmatrix}= \begin{bmatrix} d_1\\ d_2\\ \colon\\ d_{n-1}\\ d_n\\ \end{bmatrix} \tag{15}
2μ2λnλ12⋯λ2⋯μn−1⋯2μnμ1λn−12
M1M2:Mn−1Mn
=
d1d2:dn−1dn
(15)其中
λ
n
=
h
0
h
n
−
1
+
h
0
μ
n
=
h
n
−
1
h
n
−
1
+
h
0
d
n
=
6
f
[
x
0
,
x
1
]
−
f
[
x
n
−
1
,
x
n
]
h
0
+
h
n
−
1
(16)
\begin{align} \lambda_n &=\frac{h_{0}}{h_{n-1}+h_{0}} \\ \mu_n &=\frac{h_{n-1}}{h_{n-1}+h_{0}} \\ d_n &=6\frac{f\left[ x_0,x_1 \right]-f\left[ x_{n-1},x_n \right]}{h_{0}+h_{n-1}} \end{align} \tag{16}
λnμndn=hn−1+h0h0=hn−1+h0hn−1=6h0+hn−1f[x0,x1]−f[xn−1,xn](16)求解上述矩阵方程公式(13)或公式(15)即可得到插值函数公式(9)。
二、MATLAB实现
1.一阶导数边界条件
基于一阶导数已知的边界条件：
% 三次样条插值
function [yOutArr]=func_CubicSplineInterpolation1(xInArr,yInArr,dyInArr,xQuery)
% ************************************************************
% 识别误差
% ************************************************************
xInLength=length(xInArr);
yInLength=length(yInArr);
dyInLength=length(dyInArr);
xQueryLength=length(xQuery);
% 插值数组不匹配
if xInLength~=yInLength
error('Error:The xInArr=%d and the yInArr=%d are not equal in length.',xInLength,yInLength);
end
% 插值数组不匹配
if xInLength~=dyInLength
error('Error:The xInArr=%d and the dyInArr=%d are not equal in length.',xInLength,dyInLength);
end
% 插值数组不为递增数组
for i=1:1:xInLength-1
if xInArr(i)>xInArr(i+1)
error('Error:The xInArr should be a incremental array.');
end
end
% 查询数组不为递增数组
for i=1:1:xQueryLength-1
if xQuery(i)>xQuery(i+1)
error('Error:The xQuerry should be a incremental array.');
end
end
% 查询数组中的点超出插值数组的范围——查询数组有为一个点
if xQueryLength==1
% 查询数组中的点超出插值数组的范围
if xQuery(xQueryLength)>xInArr(xInLength)||xQuery(xQueryLength)<xInArr(1)
error('Error:The num of the query array xQuery should gainer than %d and lesser than %d.',xInArr(1),xInArr(xInLength));
end
% 查询数组中的点超出插值数组的范围——查询数组有多个点
else
% 查询数组中的点超出插值数组的范围
if xQuery(xQueryLength)>xInArr(xInLength)||xQuery(1)<xInArr(1)
error('Error:The num of the xQuery should gainer than %d and lesser than %d.',xInArr(1),xInArr(xInLength));
end
end
% ************************************************************
% 插值数组只有单个值时不存在区间
% ************************************************************
if xInLength==1
yOutArr=yInArr;
return;
end
% ************************************************************
% 插值数组只有两个点时直接使用三次Hermit插值
% ************************************************************
if xInLength==2
yOutArr=func_PiecewiseCubeHermitInterpolation(xInArr,yInArr,dyInArr,xQuery);
return;
end
% ************************************************************
% 插值数组至少有三个点时
% ************************************************************
prev=0;pros=0;
% 识别xQuerry查询数组在插值数组中的范围
for i=1:1:xInLength-1
% 查询数组为一个点
if xQueryLength==1
if (xQuery(xQueryLength)>xInArr(i)&&xQuery(xQueryLength)<xInArr(i+1))||...
(xQuery(xQueryLength)==xInArr(i))
prev=i;pros=i+1;
break;
end
% 查询数组有多个点
else
if(prev~=0&&pros~=0)
break;
end
if(prev==0)
if (xQuery(1)>xInArr(i)&&xQuery(1)<xInArr(i+1))||(xQuery(1)==xInArr(i))
prev=i;
end
end
if(pros==0)
if (xQuery(xQueryLength)>xInArr(xInLength-i)&&xQuery(xQueryLength)<xInArr(xInLength+1-i))||(xQuery(xQueryLength)==xInArr(xInLength+1-i))
pros=xInLength+1-i;
end
end
end
end
% ************************************************************
% 查询数组为一个点直接使用三次Hermit插值
% ************************************************************
if pros-prev==1
yOutArr=func_PiecewiseCubeHermitInterpolation(xInArr,yInArr,dyInArr,xQuery);
return;
end
% ************************************************************
% 查询数组有多个点使用三次样条插值
% ************************************************************
miu=zeros(1,pros-prev+1);
lambda=zeros(1,pros-prev+1);
M=zeros(pros-prev+1,1);
D=zeros(pros-prev+1,1);
h=zeros(1,pros-prev);
K=zeros(pros-prev+1,pros-prev+1);
% 计算λ、μ、d、h
for i=1:1:pros-prev+1
if i<=pros-prev
h(i)=xInArr(prev+i)-xInArr(prev-1+i);
end
if i==1
miu(i)=0;
lambda(i)=1;
D(i)=(6/h(i))*(func_DifferenceQuotient1(xInArr(prev-1+i:1:prev+i),yInArr(prev-1+i:1:prev+i),2)-dyInArr(prev-1+i));
K(i,i)=2;K(i,i+1)=lambda(i);
elseif i>1&&i<pros-prev+1
miu(i)=h(prev+i-2)/(h(prev+i-2)+h(prev-1+i));
lambda(i)=h(prev-1+i)/(h(prev+i-2)+h(prev-1+i));
D(i)=6*func_DifferenceQuotient1(xInArr(prev+i-2:1:prev+i),yInArr(prev+i-2:1:prev+i),3);
K(i,i-1)=miu(i);K(i,i)=2;K(i,i+1)=lambda(i);
elseif i==pros-prev+1
miu(i)=1;
lambda(i)=0;
D(i)=(6/h(i-1))*(dyInArr(prev-1+i)-func_DifferenceQuotient1(xInArr(prev+i-2:1:prev-1+i),yInArr(prev+i-2:1:prev-1+i),2));
K(i,i-1)=miu(i);K(i,i)=2;
end
end
% 计算M
M=K^(-1)*D;
% 计算yOutArr
syms x;
yOutArr=zeros(1,xQueryLength);
for i=1:1:xQueryLength
if i==1
S=M(1)*(xInArr(prev+1)-x)^3/(6*h(1))+...
M(2)*(x-xInArr(prev))^3/(6*h(1))+...
(yInArr(prev)-(M(1)*h(1)^2)/6)*(xInArr(prev+1)-x)/h(1)+...
(yInArr(prev+1)-(M(2)*h(1)^2)/6)*(x-xInArr(prev))/h(1);
yOutArr(i)=polyval(sym2poly(S),xQuery(i));
elseif i>1&&i<xQueryLength
for j=prev:1:pros-1
if (xQuery(i)>xInArr(j)&&xQuery(i)<=xInArr(j+1))||xQuery(i)==xInArr(j)
pos=j;
break;
end
end
S=M(pos-prev+1)*(xInArr(pos+1)-x)^3/(6*h(pos-prev+1))+...
M(pos-prev+2)*(x-xInArr(pos))^3/(6*h(pos-prev+1))+...
(yInArr(pos)-(M(pos-prev+1)*h(pos-prev+1)^2)/6)*(xInArr(pos+1)-x)/h(pos-prev+1)+...
(yInArr(pos+1)-(M(pos-prev+2)*h(pos-prev+1)^2)/6)*(x-xInArr(pos))/h(pos-prev+1);
yOutArr(i)=polyval(sym2poly(S),xQuery(i));
elseif i==xQueryLength
S=M(pros-prev)*(xInArr(pros)-x)^3/(6*h(pros-prev))+...
M(pros-prev+1)*(x-xInArr(pros-1))^3/(6*h(pros-prev))+...
(yInArr(pros-1)-(M(pros-prev+1)*h(pros-prev)^2)/6)*(xInArr(pros)-x)/h(pros-prev)+...
(yInArr(pros)-(M(pros-prev+1)*h(pros-prev)^2)/6)*(x-xInArr(pros-1))/h(pros-prev);
yOutArr(i)=polyval(sym2poly(S),xQuery(i));
end
end
end
2.二阶导数边界条件
基于二阶导数已知的边界条件：
% 三次样条插值
function [yOutArr]=func_CubicSplineInterpolation2(xInArr,yInArr,ddy1,ddy2,xQuery)
% ************************************************************
% 识别误差
% ************************************************************
xInLength=length(xInArr);
yInLength=length(yInArr);
xQueryLength=length(xQuery);
% 插值数组不匹配
if xInLength~=yInLength
error('Error:The xInArr=%d and the yInArr=%d are not equal in length.',xInLength,yInLength);
end
% 插值数组不为递增数组
for i=1:1:xInLength-1
if xInArr(i)>xInArr(i+1)
error('Error:The xInArr should be a incremental array.');
end
end
% 查询数组不为递增数组
for i=1:1:xQueryLength-1
if xQuery(i)>xQuery(i+1)
error('Error:The xQuerry should be a incremental array.');
end
end
% 查询数组中的点超出插值数组的范围——查询数组有为一个点
if xQueryLength==1
% 查询数组中的点超出插值数组的范围
if xQuery(xQueryLength)>xInArr(xInLength)||xQuery(xQueryLength)<xInArr(1)
error('Error:The num of the query array xQuery should gainer than %d and lesser than %d.',xInArr(1),xInArr(xInLength));
end
% 查询数组中的点超出插值数组的范围——查询数组有多个点
else
% 查询数组中的点超出插值数组的范围
if xQuery(xQueryLength)>xInArr(xInLength)||xQuery(1)<xInArr(1)
error('Error:The num of the xQuery should gainer than %d and lesser than %d.',xInArr(1),xInArr(xInLength));
end
end
% ************************************************************
% 插值数组只有单个值时不存在区间
% ************************************************************
if xInLength==1
yOutArr=yInArr;
return;
end
% ************************************************************
% 插值数组只有两个点时直接使用三次Hermit插值
% ************************************************************
if xInLength==2
yOutArr=func_PiecewiseCubeHermitInterpolation(xInArr,yInArr,dyInArr,xQuery);
return;
end
% ************************************************************
% 插值数组至少有三个点时
% ************************************************************
prev=0;pros=0;
% 识别xQuerry查询数组在插值数组中的范围
for i=1:1:xInLength-1
% 查询数组为一个点
if xQueryLength==1
if (xQuery(xQueryLength)>xInArr(i)&&xQuery(xQueryLength)<xInArr(i+1))||...
(xQuery(xQueryLength)==xInArr(i))
prev=i;pros=i+1;
break;
end
% 查询数组有多个点
else
if(prev~=0&&pros~=0)
break;
end
if(prev==0)
if (xQuery(1)>xInArr(i)&&xQuery(1)<xInArr(i+1))||(xQuery(1)==xInArr(i))
prev=i;
end
end
if(pros==0)
if (xQuery(xQueryLength)>xInArr(xInLength-i)&&xQuery(xQueryLength)<xInArr(xInLength+1-i))||(xQuery(xQueryLength)==xInArr(xInLength+1-i))
pros=xInLength+1-i;
end
end
end
end
% ************************************************************
% 查询数组在插值数组的一个区间中时直接使用三次Hermit插值
% ************************************************************
if pros-prev==1
yOutArr=func_PiecewiseCubeHermitInterpolation(xInArr,yInArr,dyInArr,xQuery);
return;
end
% ************************************************************
% 查询数组在插值数组的多个区间中时使用三次样条插值
% ************************************************************
miu=zeros(1,pros-prev);
lambda=zeros(1,pros-prev);
M=zeros(pros-prev+1,1);
D=zeros(pros-prev,1);
h=zeros(1,pros-prev);
K=zeros(pros-prev,pros-prev);
% 计算λ、μ、d、h
for i=1:1:pros-prev
h(i)=xInArr(prev+i)-xInArr(prev-1+i);
lambda(i)=h(1)/(h(i)+h(1));
miu(i)=1-lambda(i);
d(i)=6*(func_DifferenceQuotient1(xInArr(prev:1:prev+1),yInArr(prev:1:prev+1),2)-...
func_DifferenceQuotient1(xInArr(prev-1+i:1:prev+i),yInArr(prev-1+i:1:prev+i),2))/...
(h(i)+h(1));
K(i,i)=2;
K(i,max(1,mod(i+1,pros-prev+1)))=lambda(i);
K(i,max(1,mod(pros-prev-1+i,pros-prev+1)))=miu(i);
end
% 计算M
M(2:1:pros-prev+1)=K^(-1)*D;
M(1)=ddy1;
M(pros-prev+1)=ddy2;
% 计算yOutArr
syms x;
yOutArr=zeros(1,xQueryLength);
for i=1:1:xQueryLength
if i==1
S=M(1)*(xInArr(prev+1)-x)^3/(6*h(1))+...
M(2)*(x-xInArr(prev))^3/(6*h(1))+...
(yInArr(prev)-(M(1)*h(1)^2)/6)*(xInArr(prev+1)-x)/h(1)+...
(yInArr(prev+1)-(M(2)*h(1)^2)/6)*(x-xInArr(prev))/h(1);
yOutArr(i)=polyval(sym2poly(S),xQuery(i));
elseif i>1&&i<xQueryLength
for j=prev:1:pros-1
if (xQuery(i)>xInArr(j)&&xQuery(i)<=xInArr(j+1))||xQuery(i)==xInArr(j)
pos=j;
break;
end
end
S=M(pos-prev+1)*(xInArr(pos+1)-x)^3/(6*h(pos-prev+1))+...
M(pos-prev+2)*(x-xInArr(pos))^3/(6*h(pos-prev+1))+...
(yInArr(pos)-(M(pos-prev+1)*h(pos-prev+1)^2)/6)*(xInArr(pos+1)-x)/h(pos-prev+1)+...
(yInArr(pos+1)-(M(pos-prev+2)*h(pos-prev+1)^2)/6)*(x-xInArr(pos))/h(pos-prev+1);
yOutArr(i)=polyval(sym2poly(S),xQuery(i));
elseif i==xQueryLength
S=M(pros-prev)*(xInArr(pros)-x)^3/(6*h(pros-prev))+...
M(pros-prev+1)*(x-xInArr(pros-1))^3/(6*h(pros-prev))+...
(yInArr(pros-1)-(M(pros-prev+1)*h(pros-prev)^2)/6)*(xInArr(pros)-x)/h(pros-prev)+...
(yInArr(pros)-(M(pros-prev+1)*h(pros-prev)^2)/6)*(x-xInArr(pros-1))/h(pros-prev);
yOutArr(i)=polyval(sym2poly(S),xQuery(i));
end
end
end
3.测试程序
测试程序如下：
% 测试三次样条插值函数——func_CubicSplineInterpolation1、func_CubicSplineInterpolation2
clc;
clear;
close all;
x=-5.0:0.2:5.0;
y=zeros(1,length(x));
dy=zeros(1,length(x));
ddy=zeros(1,length(x));
for i=1:1:length(x)
y(i)=1/(1+x(i)^2);
dy(i)=-(2*x(i))/(x(i)^2 + 1)^2;
ddy(i)=(8*x(i)^2)/(x(i)^2 + 1)^3 - 2/(x(i)^2 + 1)^2;
end
xq=-5.0:0.05:5.0;
yq1=func_CubicSplineInterpolation1(x,y,dy,xq);
yq2=func_CubicSplineInterpolation2(x,y,ddy(1),ddy(length(x)),xq);
figure(1)
plot(x,y,'k*',xq,yq1,'r-');
xlabel('x');ylabel('y');
legend('原函数','插值函数')
figure(2)
plot(x,y,'k*',xq,yq2,'b-');
xlabel('x');ylabel('y');
legend('原函数','插值函数')
得到结果：
可以看到的是一阶导数已知的条件下插值函数曲线较为光滑，而二阶导数边界条件下的插值函数曲线稍有曲折，但不能确定是否是普遍现象。
总结
以上在本文中对三次样条插值插值两种形式的推导进行了简单的介绍，并提供了MATLAB的实现代码与测试用例。
参考文献
[1]李庆扬,王能超,易大义. 数值分析.第5版[M]. 清华大学出版社, 2008.
