关于三次样条插值,计算方法比较复杂,但是静下心来仔细研究也是可以理解的。
本文借鉴文章来源:http://www.cnki.com.cn/Article/CJFDTotal-BGZD200611035.htm
定义:
简单来说就是给定了一些在区间[a,b]的数据点{x1,x2,x3.....xn},对应函数值{y1,y2,y3.....yn},函数在[xj,xj+1] (j=1,2,...n-1此处根据你的编译器所定,matlab数组下标从1开始的)上有表达式S(x),且满足下面条件:
1. S(x)是一个三次多项式,在这里设为
2. S(xj)=yj (2-1)
3. S(xj-0)=S(xj+0) (j=2,3....,n-1)这就是保证了在非端点处的其它点连续 (2-2)
4. S'(xj-0)=S'(xj+0) (j=2,3....,n-1)一阶导数连续 (3-1)
5. S''(xj-0)=S''(xj+0) (j=2,3....,n-1)二阶导数连续 (3-2)
【注】3 4 5的区间从2开始到n-1是因为两个端点不需要判断是否连续,端点处没连续这一说。
有一个说法“n个未知数需要n个方程才能确定唯一解”,具体对不对,可以参考线性代数的知识。我们的最终目标是求出每个区间的式(1)或者函数值。 共有n-1个区间,每个区间四个参数aj,bj,cj,dj,那么就总共需要4(n-1)个求未知数。在(2-1)中给出了n个方程,(2-2)(3-1)(3-2)总共给出了3(n-2)个方程,所以依据唯一解方法可知还需要4(n-1)-3(n-2)-n=2个方程。
对于剩下的两个方程,三次样条插值给出了两个边界约束方程,刚好凑齐两个,并且有三种,可以依照自己的兴趣选择一种便于实现的。
1. 给定了端点处的一阶导数值:S'(x1)=y1' S'(xn)=yn'
2.给定了端点处的二阶导数值:S''(x1)=y1'' S''(xn)=yn''。自然边界条件:y1''=yn''=0
3.给定了一个周期性条件:如果f(x)是以b-a为周期的函数,在端点处便满足:S'(x1+0)=S'(xn-0),S''(x1+0)=S''(xn-0)
下面的推导是以边界方程1为例的:
推导过程:
(一) 利用二阶导得到一些式子:
![](https://i-blog.csdnimg.cn/blog_migrate/add977802a3837c05eea7383e0fe2465.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/c51a698e7880b21584675492ff0c9192.jpeg)
(二)依据S(x)得到一些式子
![](https://i-blog.csdnimg.cn/blog_migrate/18614b245314f5638e6bf81e3fd6f140.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/782c6d6e85b5ff1a670970d41055c16f.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/ae1b1696b9163eb007fee2ecef930480.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/9a8d5eca351bd6dd7fc179e209b45ed6.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/e466834cd47031d91ac5f089cf6d16fe.jpeg)
(三)利用一阶导数得到一些式子
![](https://i-blog.csdnimg.cn/blog_migrate/3ea73ab104d1975985d52bb944901d82.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/c45c0aeb380f6681536f468ef5f66dda.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/a81da3bdc7dec6948ddabe476028fccb.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/5f52f37aed7e950988aa4394009728be.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/62fc6d1ac28cb36c074ae98af98d8777.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/7b5f6a87da1f5d00dca30517a734566f.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/ca789d8e598363e5f98702624d2e7f21.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/c50ff92f9d68d1983b76ad19e1c0073b.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/78ee9ff27dc916c99ee93610ca7837b1.jpeg)
(四)带入边界条件
![](https://i-blog.csdnimg.cn/blog_migrate/6e4a16769150413cccafa788014fcf59.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/5d98a0417bdc0cedd3c28c5a1d614c73.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/1d07323db9046f0c499e4c9a8fe339db.jpeg)
(五)结果
![](https://i-blog.csdnimg.cn/blog_migrate/ab1583a09ea44742680705650bf80c48.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/2496e74a369b7283be9cb924fb2f12cd.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/ff042d4bf52317f39cdb5eab223a25dd.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/983416561cf25c3fcf7237f60a8deb4d.jpeg)
![](https://i-blog.csdnimg.cn/blog_migrate/691ac406320adafbf8968771688feccd.jpeg)
<pre name="code" class="plain">function f = SplineThree(x,y,dy1,dyn,x0)
%x,y为输入的已知点
%x0为待求插值点的横坐标
%dy1,dyn为约束条件,是端点处的一阶导数值
n=length(x);
for j=1:n-1
h(j)=x(j+1)-x(j);
end
%得到式子右边的结果矩阵
d(1,1)=6*((y(2)-y(1))/h(1)-dy1)/h(1); %等式(11)的结果矩阵的上端点值
d(n,1)=6*(dyn-(y(n)-y(n-1))/h(n-1))/h(n-1); %等式(11)的结果矩阵的下端点值
for i=2:n-1
u(i)=h(i-1)/(h(i-1)+h(i));
d(i,1)=6*((y(i+1)-y(i))/h(i)-(y(i)-y(i-1))/h(i-1))/(h(i-1)+h(i));
end
%得到系数矩阵
A(1,1)=2;
A(1,2)=1;
A(n,n-1)=1;
A(n,n)=2;
for i=2:n-1
A(i,i-1)=u(i);
A(i,i)=2;
A(i,i+1)=1-u(i);
end
%解方程
M=A\d;
format long
syms t;
%得到每个区间的式子
for i=1:n-1
a(i)=y(i+1)-M(i+1)*h(i)^2/6-((y(i+1)-y(i))/h(i)-(M(i+1)-M(i))*h(i)/6)*x(i+1);
b(i)=((y(i+1)-y(i))/h(i)-(M(i+1)-M(i))*h(i)/6)*t;
c(i)=(t-x(i))^3*M(i+1)/(6*h(i));
e(i)=(x(i+1)-t)^3*M(i)/(6*h(i));
f(i)=a(i)+b(i)+c(i)+e(i);
%f(i)=M(i)*(x(i+1)-t)^3/(6*h(i))+M(i+1)*(t-x(i))^3/(6*h(i))+(y(i)-M(i)*h(i)^2/6)*(x(i+1)-t)/h(i)+(y(i+1)-x(i+1)*h(i)^2/6)*(t-x(i))/h(i);
% f(i)=((x(j+1)-x)^3)*M(i)/(6*h(i))+((x-x(i))^3)*M(i+1)/(6*h(i))+(y(i)-M(i)*(h(i)^2)/6)*((x(i+1)-x)/h(i))+(y(i+1)-(M(i+1)*(h(i)^2)/6))*((x-x(i))/h(i));
end
%化简
f=collect(f);
f=vpa(f,6);
if(nargin==5)
nn=length(x0);
for i=1:nn
for j=1:n-1
if(x0(i)>=x(j)&x0(i)<=x(j+1))
yynum(i)=subs(f(j),'t',x0(i)); %计算插值点的函数值.subs是替换函数,把x0用t替换
end
end
end
f=yynum;
else
f=collect(f); %将插值多项式展开
f=vpa(f,6); %将插值多项式的系数化成6位精度的小数
end
end
SplineThree.m
<pre name="code" class="plain">% x=[-1.5 0 1 2];
% y=[0.125 -1 1 9];
% dy1=0.75;
% dyn=14;
% x0=1.5;
% answer=SplineThree(x,y,dy1,dyn);
%
% X=[-1.5 0 1 2];
% Y=[0.125 -1 1 9];
% dY=[0.75 14]
% m=5;
% spline3(X,Y,dY,x0,m)
X=0:2*pi;
Y=sin(X);
dY=[1,1];
dy1=1;
dyn=1;
xx=0:0.5:6;
m=5;
% for i=1:length(xx)
% spline3(X,Y,dY,xx,m);
% end
yy=SplineThree(X,Y,dy1,dyn,xx);
plot(xx,yy,'r')
![](https://i-blog.csdnimg.cn/blog_migrate/f70d3b12f5b2929242e235d672f43318.jpeg)
本文已经同步到微信公众号中,公众号与本博客将持续同步更新运动捕捉、机器学习、深度学习、计算机视觉算法,敬请关注