(精选)三次样条插值的MATLAB实现

1、MATLAB程序设计期中考查在许多问题中，通常根据实验、观测或经验得到的函数表或离散点上的信息，去研究分析函数的有关特性。其中插值法是一种最基本的方法，以下给出最基本的插值问题三次样条插值的基本提法：对插值区间Cz,b进行划分：axxxb,函数y=f(x)在节点01nx上的值y=f(x)(=0,1,2,n)，并且如果函数S(jc)在每个小区间lx,x上iiiii+1是三次多项式，于L,b上有二阶连续导数，则称SG)是la,b上的三次样条函数,如果S(x)在节点x上还满足条件iSCx)=y(i-0,1,n)ii则称S(x)为三次样条插值函数。三次样条插值问题提法：对L,b上给定的数表如下.x1求

2、一个分段三次多项式函数S6)满足插值条件S(x)=y(i=0,1,n)式，并在ii插值区间L,b上有二阶连续导数。这就需要推导三次样条插值公式:记f(x)在节点x处的值为fCx)=m(i=0,1,n)(这不是给定插值问题数iii表中的已知值)。在每个小区间tx,x利用三次Hermite插值公式，得三次插值ii+1公式：SCx)=aCx)y+aCx)y+PCx)m+Pm，xgtx,x。为了得到这个公iii+1i+1iiii+1i+1ii+1式需要4n个条件：(1).非端点处的界点有2n个；(2).阶导数连续有n-1个条件；(3).二阶导TOC o 1-5 h z数连续有n-1个条件，其中边界条件

3、：QS心)=mS心)=m00nnQS(x)=aS(x)=a00nQ3aSCx)+aSCx)=PaSCx)+aSCx)=P304005n6n1Qy=ySf(x+0)=Sf(x-0)0n0S心+0)=S心-0)0n其中：aCLij0,i丰j1,i=jar(x)=0p(x)=0且(i,j=0,1)。ijij0,i丰j，i=j为对应变量的一阶导数。其推导过程如下：为了确定m的值,ic()(x-xSVx丿=i+i把S(x)展开为:这里hi(x+2(x-x)(xxCh+2(x-x)ih3i)2(x-x)(x-x)2(xiy+iii+1ih3i-x)i+mi+1-xi+1m+iih2ih2ii=x-x，对S

4、(x)连续求两次导，得:i+1i6x2x4x6x4x2xi+1m+i.h2ih2ii考虑SG)在节点x处的右极限值，得：iS()=-i+1mi+1yi+16(x+x+ii+1h3i-y)。于是i6(+(yh2iS心+0)=-mmihihii同理，在相邻小区间L,x上可得S心)的表达式为:i-1ii+1i+1-y)。i()6x-2x-4xSVx=i-1-m+h2i-1i-16x-4x-2x6(x+xi-1Lm+i-1ih2ih2i-1i-1iyi1)ii1及S)在节点x处的左极限值为:i6h2i-1的连续性条件S心+0)=S心-0),这里i=1,2,.n-1，有下式成立:iis心-o)=2ihi

5、-14m+m-i-1hii-1(y-y)。利用SG)二阶导数于节点x处ii-1i1丄2m+2hi1i1r11)+m+mIhh.)i-1i1=3ihi+1i、yy丄yyi+lL+ii-1h2h2丿ii111，用丄+丄除等式两hhi1i、丄l厂、亠边，并注意y二fiiy-yi+1ihi二ftx,x,上式可简记为:ii+1九m+2m+pm=gii-1iii+1i(i-1,2,n-1)hih+hi-1i17hP=1九=i-1iih+hi-1r2p0.00九mg1111九2卩.000mg2222=00072mgn-1n-1n-1n-100072mgnnn(*)g=3(7fix,x+pftx,xDiii-

6、1iiii+1最后求得mm的线性方程组为:1nn通过以上复杂的求解和迭代，就可以求解出插值函数的近似表达式。得出来的表达式就可以用MATLAB软件来求解。具体求解过程如下:已知n对数据点(x,y),(x,y)(x,y),(x,y),，假设函数关系为112233nny=f(x),但解析式不确定，数据插值就是构造函数关系式y=g(x),使x(=1,2,3,n)，满足关系g(x)=f(x)。iii例题:求满足下面函数表所给出的插值条件的三次自然样条函数。分析:表中所列出的是函数对点，首先要把对应的插值函数求出来，再用MATLAB软件来求区间1,5上间隔为0.5的各点的值。求解过程如下:因自然样条插值

7、函数的边界条件为Sn(x)=S心)=0,0n这里n=3，故确定m,m,m,m的方程组形式形如上面的(*)式，其中系数九，卩0123ii和g可按如下步骤进行：i0h:1,1;33g：ig1=3(pflx,x+九flx,x)=;1121012g2=3(gfL,x+九fL,xD=2232122g0=3fWxi-守y：=6;g3=3fL,x+23艺y；=-6.将上述参数带入(*)式，得到以下方程组：100丄0-612解得：1744由公式S(x)=-x4+1)2th+2(x-x)-ih3-x-i)2th+2(x-x)-i4+1h3yi+1-x-x4+1ih2h2i+1可知，S(x)=一一x3x-1,l1,4-x345103x-33,xel4,5MATLAB由所求出的表达式可知区间1,5可分为1,4ul4,5，对两个区间分别用命令即可：针对第一个区间127y一8x3+8x2+4x-1;其图像如下命令如下：x=l:4;y=(-l/8)*x3+(2/8)*x.2+(7/4)*x-l;xi=l:0.5:4;y1=interp1(x,y,xi,spline)其运行结果如下：y1=Columns1through60.87501.