紧压三次样条插值程序讲解.doc

文档介绍：求压紧三次样条函数的函数程序
function S=liti05_7(X,Y,dx0,dxn)
%Input -X is the 1xn abscissa vector
% -Y is the 1xn ordinate vector
% -dx0=S'(x0) first derivative boundary condition
% -dxn=S'(xn) first derivative boundary condition
%Output -S:rows of S are the confficients,indescending order,for the cubic interplants
N=length(X)-1; %求当前问题的规模数--的最大下标
H=diff(X); %求X的差分 h0 h1 …… hn-1
D=diff(Y)./H; %求 y 对 x 的一阶差商 d0 d1 …… dn-1
A=H(2:N-1); %为求线性方程组系数矩阵上次对角元做准备
B=2*(H(1:N-1)+H(2:N)); %求线性方程组系数矩阵主对角元做准备
C=H(2:N-1); %求线性方程组系数矩阵上次对角元
U=6*diff(D); %为求线性方程组常数列向量做准备
%压紧样条端点约束
B(1)=B(1)-H(1)/2; %3h0/2+2h1
U(1)=U(1)-3*(D(1)-dx0); %修改 u(1)
B(N-1)=B(N-1)-H(N)/2; % 2h(N-2)+3*h(N-1)/2
U(N-1)=U(N-1)-3*(dxn-D(N)); %修改 u(N-1)
A %输出线性方程组系数矩阵上次对角元向量
B %输出线性方程组系数矩阵主对角元向量
C %输出线性方程组系数矩阵下次对角元向量
U %输出线性方程组常数列向量
% 下面开始解三对角线行方程组
% 首先转化为上三角线性方程组
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
%计算端点x0 xn上的二阶导数
M(1)=3*(D(1)-dx0)/H(1)-M(2)/2;
M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;
M % 输出各节点上的二阶导数
for k=0:N-1
%计算第k个多项式的系数
S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1)); %算(x-X(k))三次项系数
S(k+1,2)=M(k+1)/2; %算(x-X(k))二次项系数
S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;%算(x-X(k))一次项系数
S(k+1,4)=Y(k+1); %算(x- X(k))零次项系数
end
hold on
plot(X,Y,'kO') % 画样点
s