紧压三次样条插值程序讲解.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 interplantsN=length(X)-1; %求当前问题的规模数-的最大下标 H=diff(X); %求X的差分 h0 h1 hn-1D=diff(Y)./H; %求 y 对 x 的一阶差商 d0 d1 dn-1A=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+2h1U(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)/2U(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)）零次项系数endhold onplot(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)endgrid onxlabel(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 = 1B = 3.5000 3.5000C = 1U = 5.1000 -10.5000M = -0.3600 2.5200 -3.7200 0.3