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

求压紧三次样条函数的函数程序

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; %求当前问题的规模数--x的最大下标

H=diff(X); %求X的差分 h0 h1…… h n-1

D=diff(Y)./H; %求 y 对 x 的一阶差商 d0 d1…… d n-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 x n上的二阶导数

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') % 画样点

style=['r', 'g', 'b']; % 由于例题中N等于3，style只构造了三个元素

for k=1:N

%画第k个区间上的三次函数曲线段

x1=X(k):0.01:X(k+1);

y1=polyval(S(k,:),x1-X(k)); %算关于（x- X(k)）三次多项式的值

plot(x1,y1,style(k))

end

grid on

xlabel('x');

ylabel('y');

hold off

例求压紧三次样条曲线，经过点（0，0），（1，0.5），（2，2.0），（3，1.5），一阶导数的边界条件为S’(0)=0.2和S’(3)=-1

针对上例，程序运行结果如下

?x=[0 1 2 3];y=[0 0.5 2.0 1.5];

?dx0=0.2;dxn=-1;

?s=liti05_7(x,y,dx0,dxn)

A =

1

B =

3.5000 3.5000

C =

1

U =

5.1000 -10.5000

M =

-0.3600 2.5200 -3.7200 0.3600

s =

0.4800 -0.1800 0.2000 0

-1.0400 1.2600 1.2800 0.5000

0.6800 -1.8600 0.6800 2.0000

y

x