(精選)三次樣條插值的MATLAB實現

1、MATLAB程序設計期中考查在許多問題中，通常根據實驗、觀測或經驗得到的函數表或離散點上的信息，去研究分析函數的有關特性。其中插值法是一種最基本的方法，以下給出最基本的插值問題三次樣條插值的基本提法：對插值區(qū)間Cz,b進行劃分：axxxb,函數y=f(x)在節(jié)點01nx上的值y=f(x)(=0,1,2,n)，并且如果函數S(jc)在每個小區(qū)間lx,x上iiiii+1是三次多項式，于L,b上有二階連續(xù)導數，則稱SG)是la,b上的三次樣條函數,如果S(x)在節(jié)點x上還滿足條件iSCx)=y(i-0,1,n)ii則稱S(x)為三次樣條插值函數。三次樣條插值問題提法：對L,b上給定的數表如下.x1求

2、一個分段三次多項式函數S6)滿足插值條件S(x)=y(i=0,1,n)式，并在ii插值區(qū)間L,b上有二階連續(xù)導數。這就需要推導三次樣條插值公式:記f(x)在節(jié)點x處的值為fCx)=m(i=0,1,n)(這不是給定插值問題數iii表中的已知值)。在每個小區(qū)間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).階導數連續(xù)有n-1個條件；(3).二階導TOC o 1-5 h z數連續(xù)有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)連續(xù)求兩次導，得:i+1i6x2x4x6x4x2xi+1m+i.h2ih2ii考慮SG)在節(jié)點x處的右極限值，得：iS()=-i+1mi+1yi+16(x+x+ii+1h3i-y)。于是i6(+(yh2iS心+0)=-mmihihii同理，在相鄰小區(qū)間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)在節(jié)點x處的左極限值為:i6h2i-1的連續(xù)性條件S心+0)=S心-0),這里i=1,2,.n-1，有下式成立:iis心-o)=2ihi

5、-14m+m-i-1hii-1(y-y)。利用SG)二階導數于節(jié)點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軟件來求區(qū)間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由所求出的表達式可知區(qū)間1,5可分為1,4ul4,5，對兩個區(qū)間分別用命令即可：針對第一個區(qū)間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.