东南大学_数值分析_第四章_三次样条插值.doc

东南大学《数值分析》上机练习——算法与程序设计实验报告
第四章 多项式插值与函数最佳逼近
——曲线拟合之3次样条插值
*****（学号） *****（姓名）
上机题目要求见教材P195,37题。
一、算法原理
题目要求编写第一边界条件的3次样条插值函数的通用程序，同时根据汽车门曲线值点构造三次紧压样条曲线函数。其基本原理如下
定义 设有N+1个点，其中。如果存在N个三次多项式，系数为满足如下性质：
(1)
(2)
则成为三次样条函数。
现证明其存在：
由于是分段三次多项式，其二阶导数是在区间内是分段线性的。根据线性拉格朗日插值可以表示为：
(3)
用代入上式，得
(4)
将上式积分两次，会引入两个积分常数，可得到如下形式：
(5)
将代入上式，并利用可得两个方程：
(6)
求解，并将所得的结果带入方程得
(7)
求式（7）的导数，并化简得
(8)
由上述方程可得如下方程
(9)
(10)
(11)
重组上述方程，得三角线性方程组，表示为
(12)
该式具有严格对角优势。算出系数后可由如下公式计算的样条系数。
(13)
二、流程图
开始
输入数据及边界条件
计算H，D，U
利用端点约束计算
求各段的样条系数
结束
用已知的N+1个点构造三次紧压样条曲线的问题。其通用程序流程图1所示。
图1关于求解三阶样条曲线算法流程图
具体步骤如下：
1) 计算，，
2) 联系端点约束条件求解样条函数的系数。
三、计算代码
核心代码
for k=2:N-1
temp=A(k-1)/B(k-1);
B(k)=B(k)-temp*C(k-1);
U(k)=U(k)-temp*U(k-1);
end
for k=N-2:-1:1
M(k+1)=(U(k)-C(k)*M(k+2))/B(k);
for k=0:N-1
S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));
S(k+1,2)=M(k+1)/2;
S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;
S(k+1,4)=Y(k+1);
end
完整代码
function S=Three_fit(X,Y,dx0,dxn)
%Input - X is a vector that contains a list of abscissas
% - Y is a vector that contains a list of ordinates
% - dx0=S'(x0) first derivative boundary condition
% - dxn=S'(xn) first derivative boundary condition
%Output - S rows of S are the coefficients , in descending order ,for
% the cubic interpolants
N=length(X)-1;
H=diff(X);
D=diff(Y)./H;
A=H(2:N-1);
B=2*(H(1:N-1)+H(2:N));
C=H(2:N);
U=6*diff(D);
%clamped spline endpoint constraints
B(1)=B(1)-H(1)/2;
U(1)=U(1)-3*(D(1)-dx0);
B(N-1)=B(N-1)-H(N)/2;
U(N-1)=U(N-1)-3*(dxn-D(N));
%Gauss elimination to solve mk
for k=2:N-1
temp=A(k-1)/B(k-1);
B(k)=B(k)-temp*C(k-1);
U(k)=U(k)-temp*U(k-1);
end
M(N)=U(N-1)/B(N-1);
for k=N-2:-1:1
M(k+1)=(U(k)-C(k)*M(k+2))/B(k);
end
M(1)=3*(D(1)-dx0)/H(1)-M(2)/2;
M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;
for k=0:N-1
S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));
S(k+1,2)=M(k+1)/2;
S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;
S(k+1,4)=Y(k+1);
end
四、计算结果及分析
根据题目数据调用函数如下
X=0:10;
Y=[2.51,3.30,4.04,4.70,5.22,5.54,5.78,5.40,5.57,5.70,5.80];
dx0=0.8;dxn=0.2;
S=Three_fit(X,Y,dx0,dxn)
plot(X,Y,'o'),hold on
yprint=[];
for ii=1:10
x=X(ii):0.1:X(ii+1);
y=polyval(S(ii,:),x-X(ii));%分段拟合
yprint=[yprint,polyval(S(ii,:),0.5)];plot(x,y)
end
xlabel('x','FontSize',20),ylabel('y','FontSize',20)
title('三次样条插值','FontSize',22)
hold off
disp('-------------------------------------------------------------------------------')
disp(' x(i) = 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5')
disp('-------------------------------------------------------------------------------')
disp(['S(i+0.5)= ' num2str(yprint(1)) ' ' num2str(yprint(2)) ' ' ...
num2str(yprint(3)) ' ' num2str(yprint(4)) ' ' num2str(yprint(5))...
' ' num2str(yprint(6)) ' ' num2str(yprint(7)) ' ' num2str(yprint(8))...
' ' num2str(yprint(9)) ' ' num2str(yprint(10))])
disp('-------------------------------------------------------------------------------')
输出三次样条函数的系数
S = -0.0085 -0.0015 0.8000 2.5100
-0.0045 -0.0270 0.7715 3.3000
-0.0037 -0.0404 0.7041 4.0400
-0.0409 -0.0514 0.6123 4.7000
0.1074 -0.1741 0.3868 5.2200
-0.2685 0.1479 0.3606 5.5400
0.4266 -0.6575 -0.1491 5.7800
-0.2679 0.6222 -0.1844 5.4000
0.0549 -0.1814 0.2565 5.5700
0.0584 -0.0168 0.0584 5.7000
打印出的值
-------------------------------------------------------------------------------
x(i) = 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5
-------------------------------------------------------------------------------
S(i+0.5)= 2.9086 3.6784 4.3815 4.9882 5.3833 5.7237 5.5944 5.4299 5.6598 5.7323
-------------------------------------------------------------------------------
输出拟合曲线:

Fig2 三次样条曲线
分析：由图形可以看出，三阶样条曲线非常光滑，用该法能获得很好的拟合效果。
五、结论
对数据进行多项式拟合在CAD，CAM和计算机图形处理软件中有许多应用。我们常常需要画出经过多个点的无误差的光滑的曲线。从数学的角度上看就是要求图像及其一阶，二阶导数在分段区间内连续。本次报告讨论了有关三阶样条函数的原理及应用，该法可使函数图像二阶连续可导，反映出曲线的光滑。该法具有重要的工程应用意义。