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

《紧压三次样条插值程序讲解》由会员上传分享，免费在线阅读，更多相关内容在应用文档-天天文库。

1、求压紧三次样条函数的函数程序functionS=liti05_7(X,Y,dx0,dxn)%Input-Xisthe1xnabscissavector%-Yisthe1xnordinatevector%-dx0=S'(x0)firstderivativeboundarycondition%-dxn=S'(xn)firstderivativeboundarycondition%Output-S:rowsofSaretheconfficients,indescendingorder,forthecubicinterplantsN=length(X)-1;%求当前问题的规模数--的最大下标H

2、=diff(X);%求X的差分h0h1……hn-1D=diff(Y)./H;%求y对x的一阶差商d0d1……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)

3、=U(N-1)-3*(dxn-D(N));%修改u(N-1)A%输出线性方程组系数矩阵上次对角元向量B%输出线性方程组系数矩阵主对角元向量C%输出线性方程组系数矩阵下次对角元向量U%输出线性方程组常数列向量%下面开始解三对角线行方程组%首先转化为上三角线性方程组fork=2:N-1temp=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);fork=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2))/B(k);end%计算端点x0xn上的二阶导数M(1

4、)=3*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;M%输出各节点上的二阶导数fork=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)）零次项系数endholdonplot(X,Y,'kO')%画样点st

5、yle=['r','g','b'];%由于例题中N等于3，style只构造了三个元素fork=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))endgridonxlabel('x');ylabel('y');holdoff例求压紧三次样条曲线，经过点（0，0），（1，0.5），（2，2.0），（3，1.5），一阶导数的边界条件为S’(0)=0.2和S’(3)=-1针对上例，程序运行结果如下?x=[0123];y=[00.52.

6、01.5];?dx0=0.2;dxn=-1;?s=liti05_7(x,y,dx0,dxn)A=1B=3.50003.5000C=1U=5.1000-10.5000M=-0.36002.5200-3.72000.3600s=0.4800-0.18000.20000-1.04001.26001.28000.50000.6800-1.86000.68002.0000