自编的三次样条插值matlab程序(含多种边界条件)
数值计算第二次大作业
——验证三次样条函数插值是否有几何不变性
（1）给定的插值条件如下：
i 0 1 2 3 4 5 6 7
Xi 8.125 8.4 9.0 9.485 9.6 9.959 10.17 10.2
Yi 0.0774 0.099 0.28 0.60 0.708 1.200 1.800 2.177
端点边界条件为第一类边界条件（给定一阶导数）：
.Y' .
0.01087
.0.
'
Y .100
.7
.
三次样条函数的构造过程如下：
设
x1 .
x2 .
x3 .
xn.1 .xn 共
n个插值节点，则经过数据点
.
x1, y1 .
,.
x2, y2 .
,.,.xn , yn .
的三次样条
S .
x.
是一组三次多项式：
.S .........231111111112,,xabxxcxxdxxxxx.........,
.
.
23
.
.
.
a ..
x .
x .
c ..
x .
d .x .
x , x ..x , x ,
.S2 x 2 b22 .
2 x 2 .
22 .
23 .
.
(1.1)
.
.
.
S .x .
x .3, x ..x , x .
.
n.
......21111111nnnnnnxabxxcxxd.............
n .
n.1 n.1
.
由节点处的连续性可知：
S .
x ..
y , S .
x ..
y ,i .1, 2, .n .1.
ii iii.1 i.1
.a .
yi , .1, 2, .n .1.
ii
.
.y y ......2321121121121bxxcxxdxx......., (1.2)
..
.
.
.y .
y .
b .
x .
x ..
c .x .
x .2 .
d .x .
x .3
.
nn.1 n.1 nn.1 n.1 nn.1 n.1 nn.1
由节点处的一阶与二阶光滑性可知：
Si
''""
11,iiiiiixSxSxS.....
..
.
.
..xi .
,i .1, 2, ., n.
.
2
''
0 .
S .x ..
S .x ..
b .
2c .x .
x ..
3d .x .
x ..
b
12 221121 121 2
.
..
.
2
''
(1.3)
.0 .
Sn.2 .
xn.1 ..
Sn.1 .
xn.1 ..
bn.2 .
2cn.2 .xn.1 .
xn.2 ..
3dn.2 .xn.1 .
xn.2 ..
bn.1
..
""
.0 .
S .
x ..
S .
x ..
2c .
6d .x .
x ..
2c
1222 1121 2
.
.
.
.
.0 .
Sn
""
211122nnnnxSxc........
..
..
6dn.2 .
xn.1 .
xn.2 ..
2cn.1
又设
cn .Sn
"
1..
xn .
2 ，记
..
x .
x , ..
y .
yi , .1, 2, ., n .1 ，则由(1.3)可得：
ii.1 ii i.1 i
c .
c
i.1 i
di .
, i .
1, 2, ., n .1. (1.4)
3.i
从（1.2）解得：
bi .
.i .
ci.i .
di.i
212iiiicc....
.
.
.
., i .1, 2, ., n .1. (1.5)
.i .i 3
将(1.4)与(1.5)代入(1.3)得：
.
..2 .1 .
..c .
2..1 ..
.
2 ..2c3 .
3
11 2 c ...,
...2 .1 .
.
..
(1.6)
.
..n.1 .n.2 .
..n.2cn.2 .
2..n.2 ..n.1 .cn.1 ..n.1cn .
3...
.
.
..n.1 .n.2 .
增加两个端点边界条件，因为
2c ""
1111,2nnSxcS....
.
.
xn .
，故有：
1.第零类边界条件：自然样条，
c1 .
0, cn .
0.
2.第一类边界条件：给定端点一阶导数，设
S ''
1111,nxvS...
..
xn ..vn ，则有：
.2.c ..c .
3..
..
v .
.
.
..
.
.n .
11121111111123nnnnnnccv..........
.
3.第二类边界条件：给定端点二阶导数，设
S""
1111,nxvS...
..
xn ..vn ，则有：
.2c1 .
v1.
2c .
v
.nn
ci .
c ..
结合(1.6)及
所给的边界条件即可解出
.ci .,而
di ..
311,2iiiiiibcc....
.i.i3
.
.
，故
可得到最终各个子区间上的三次样条函数。
根据以上过程进行
matlab编程，编写三次样条
spline3函数，具体见附录。
因此所编函数可对第一题求解：
clear;clc;
format short g;
x1=[8.125 8.4 9 9.485 9.6 9.959 10.166 10.2]';
y1=[0.0774 0.099 0.28 0.6 0.708 1.2 1.8 2.177]';
u1=0.01087;un=100;
xx1=[x1(1):0.001:x1(end)]';
[yy1 b1 c1 d1]=spline3(x1,y1,xx1,1,u1,un);
fprintf('\t\tb1\t\tc1\t\td1\n');
disp([b1 c1(1:end-1,1) d1]);
plot(x1,y1,'bo',xx1,yy1,'--k');
grid on;
b1 c1
0.01087
0.17405
0.2878
1.5294
-0.56548
12.794
-28.949
d1
0.14489
0.4485
-0.25891
2.8188
-21.035
58.247
-259.9
0.368
-0.393
2.1153
-69.141
73.614
-512.32
42279
8 8.5 9 9.5 10 10.5
0
0.5
1
1.5
2
2.5
3
图表一、三次样条曲线
（2）坐标轴逆时针旋转45度，相当于节点顺时针旋转45度。设
.
, .T 为旋转前的
xy
''
坐标，
.
x , y.T 为旋转后的坐标，则有：
x ' cos ..sin .
x
...
...
.
...
...
y ' sin .
cos .
y
...
...
故旋转后的节点坐标为：
theta=-pi/4;
for i=1:length(x1)
x2(i)=cos(theta)*x1(i)-sin(theta)*y1(i);
y2(i)=sin(theta)*x1(i)+cos(theta)*y1(i);
end
fprintf('\t\t\tx2\t\t\ty2\n');
disp([x2' y2']);
x2 y2
5.8 -5.6905
6.0097 -5.8697
6.562 -6.166
7.1312 -6.2826
7.2889 -6.2876
7.8906 -6.1935
8.4612 -5.9157
8.7519 -5.6731
端点处的一阶导数为：
v1=(u1+tan(theta))/(1-u1*tan(theta));
vn=(un+tan(theta))/(1-un*tan(theta));
fprintf('\t\t\tv1\t\t\tvn\n');
disp([v1 vn]);
v1 vn
-0.97849 0.9802
则旋转后的三次样条的系数及图像为：
xx2=[x2(1):0.001:x2(end)]';
[yy2 b2 c2 d2]=spline3(x2,y2,xx2,1,v1,vn);
fprintf('\t\t\tb2\t\t\tc2\t\t\td2\n');
disp([b2 c2(1:end-1,1) d2]);
plot(x2,y2,'*b',xx2,yy2,'-.k');
grid on;
b2 c2 d2
-0.97849 0.67221 -0.38277
-0.74704 0.43138 -0.090754
-0.35362 0.28102 -0.034909
-0.067629 0.22141 0.053338
0.0061747 0.24664 0.0046897
0.3081 0.2551 0.10233
0.6992 0.43028 0.12195
-
-
-
-
-
-
-
-
-
5.5 6 6.5 7 7.5 8 8.5 9
-6.4
-6.3
-6.2
-6.1
-6
-5.9
-5.8
-5.7
-5.6
图表二、旋转后的三次样条曲线
（3）将第（1）题中所得的样条曲线整体旋转-45度（即顺时针旋转45度），并与
第（2）题的曲线画在同一幅图上，得：
for i=1:length(xx1)
xx3(i)=cos(theta)*xx1(i)-sin(theta)*yy1(i);
yy3(i)=sin(theta)*xx1(i)+cos(theta)*yy1(i);
end
plot(xx3,yy3,'--',xx2,yy2,'-',x2,y2,'ok');grid on;
legend('旋转前','旋转后');
-
-
-
-
-
-
-
-
-
5.5 6 6.5 7 7.5 8 8.5 9 9.5
-6.8
-6.6
-6.4
-6.2
-6
-5.8
-5.6
-5.4
-5.2
旋转前
旋转后
图表三、旋转前后样条曲线几何比较
比较上图中的两条曲
线，易知曲线不重合，故有以下结论：
三
三三三次
次次次样
样样样条
条条条函
函函函数
数数数插
插插插值
值值值不
不不不具
具具具备
备备备几
几几几何
何何何不
不不不变
变变变性
性性性！
！！！
附
附附附录
录录录
三次样条插值函数matlab的m文件程序：
function [yy b c d]=spline3(x,y,xx,flag,vl,vr)
%三次样条插值函数
% (x,y)为插值节点，
xx为插值点；
% flag表端点边界条件类型：
% flag=0 :自然样条
(端点二阶导数为
0)；
% flag=1 :第一类边界条件
(端点一阶导数给定
)；
% flag=2 :第二类边界条件
(端点二阶导数给定
)；
% vl,vr表左右端点处的在边界条件值。
%样条函数为：
Si(x)=yi+bi*(x-xi)+ci*(x-xi)^2+di*(x-xi)^3
% b,c,d分别为各子区间上的系数值
% yy表插值点处的函数值
.
if length(x)~=length(y)
error('输入数据应成对
!');
end
n=length(x);
a=zeros(n-1,1);
b=a;d=a;dx=a;dy=a;
A=zeros(n);B=zeros(n,1);
for i=1:n-1
a(i)=y(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)=dx(i-1);
A(i,i)=2*(dx(i-1)+dx(i));
A(i,i+1)=dx(i);
B(i,1)=3*(dy(i)/dx(i)-dy(i-1)/dx(i-1));
end
%自然样条端点条件
(端点二阶导数为零
)%
if flag==0
A(1,1)=1;
A(n,n)=1;
end
% ---------------------------------%
%端点一阶导数条件
%
if flag==1
A(1,1)=2*dx(1);A(1,2)=dx(1);
A(n,n-1)=dx(n-1);A(n,n)=2*dx(n-1);
B(1,1)=3*(dy(1)/dx(1)-vl);
B(n,1)=3*(vr-dy(n-1)/dx(n-1));
end
% ---------------%
%端点二阶导数条件
%
if flag==2
A(1,1)=2;A(n,n)=2;
B(1,1)=vl;B(n,1)=vr;
end
% ---------------%
c=A\B;
for i=1:n-1
d(i)=(c(i+1)-c(i))/(3*dx(i));
b(i)=dy(i)/dx(i)-dx(i)*(2*c(i)+c(i+1))/3;
end
[mm nn]=size(xx);
yy=zeros(mm,nn);
for i=1:mm*nn
for ii=1:n-1
if xx(i)>=x(ii) & xx(i)
break;
elseif xx(i)==x(n)
j=n-1;
end
end
yy(i)=a(j)+b(j)*(xx(i)-x(j))+c(j)*(xx(i)-x(j))^2+d(j)*(xx(i)-x(j))^3;
end
end