三次樣條插值試驗報告

1、三次樣條插值實驗報告專業(yè)班級學號姓名一、實驗內(nèi)容和要求1、閱讀上面的文字和程序，試運行，檢驗程產(chǎn)和上面敘述的正確性。2、閱讀上面的MATLAB程序：查資料，了解各MATLAB語句及命令。3、畫程序流程圖，理解并描述算法。4、修改上面的程序，能根據(jù)給定數(shù)據(jù)點，求(1)自然樣條插值，邊界S"(a)=O,S"(b)=O0回(2)第二類邊界條件,S"(a)和S”(b)是確定的。5、使用上面的程序，根據(jù)數(shù)據(jù)點(0,1),(1,0),(2,0),(3,1),(4,2),(5.2),和(6,1),求三種不同的三次樣條插值，其中S'(O)=-0.6,5*(6)=-1.8,

2、<(0)=1,S"(6)=-l：S"(0),S"(6)=0.在同一坐標系中，畫出這3個三次樣條插值和這些數(shù)據(jù)點。6、寫實驗報告(實驗內(nèi)容+算法描述+程序+寫成分段函數(shù)的結(jié)果描述+截圖)。二、算法說明定義：設(shè)有N+1個點，箕金;H如果存?8=%<%<,Y4=bN個三次多項式a(號數(shù)為，滿是姆閔調(diào)心(x),8k3(x)(1)S(X)=Sk(x)=Sk0+8kl(X-5)+、式x-X1t尸+(X-%)3xw區(qū),4j,k=0,1,N-1(2)S(%)=ykk=0,l,，N(3)lim8k(x)=lim(x)xf埼xtWlimS；(x)=lim(x)10X

3、T、lim鼠(x)=lim"(x)XT&k=l,2,*,n-lk=l,2,n-lk=l,2/«,n-l則稱函數(shù)為8次樣條函數(shù),因為是分段3；欠多項式，所以在區(qū)間上足分段線性的la,b(1)Pl+sd初隊S”（）神戶%-%(2)S；(x)=Mk+Mk+1<x<xk+pk=0,1,-sN-1.將（2）式積分兩次，會引入兩個積分常數(shù)，并得到（3）s1t（X）=（xk+1-x）3+乎（X）3+Pk（k一X）+Q（X-%）6hk6hk將代入加，并使用，可分別得到賴自班腥（看”）=丫皿（4）yk=?供pA，yk4i=乎h；+qAoo求解這兩個方程，求出，而以糊這些值

4、代入方程（3）中，可得到三次多項式方程：Sk（x）=-,（4*-x）3+等（X-4尸+（興-x）+（管-（Xf）（5）"表達式5可以簡化成只包含未知系數(shù)的形題J為求這些值，必須使用（5）式的導數(shù)，即在處計算（6）,并簡化結(jié)果可得到,其河新廣孕L號4+44=巫p同理，在式（6）中用并計算消處的解戰(zhàn)q（得（8）利用菰短鳴一導野誨綠監(jiān)方程（6）、（7）,可得到包含Mi尚甯喙聯(lián)系式（9）14TMl+2（%t+hk）Mk+hkMk+1=4其中74=6（dk-dk_i）,k=l,2,、N-L方程組（9）中的未知數(shù)是要求的值，饌篇的項可通過數(shù)據(jù)點集進行簡里隊沖硼到的常量。因此方程（9）是包含N+1

5、個未知數(shù)，具有N-1個線性方程的不定方程組。所以需要另外兩個方程組才能求解，即邊界條件。如果已知，財與）$（、）M。=（-SXq）-%2（10）hN-l2根據(jù)（耳£,410）求出后Mk時利用下面的公式計算的樣條條幽鼠,。=yi，%=4-%Qm1Me）,產(chǎn)學“，dk=必一%三、源程序csfitl:第一類邊界條件functionS=csfitl(X,Y,dxO,dxn)N=length(X)-l;H=diff(X);D=diff(Y)./H;A=H(2:N-1);B=2*(H(1;N-1)+H(2:N);C=H(2:N);C=H(2:N);U=6-diff(D);U(l)=U(l)-3*

6、(D(l)-dxO);B(N-l)=B(N-l)-H(N)/2;U(N-l)=U(N-l)-3*(dxn-D(N);fork=2:N4temp=A(k-l)/B(k-l);B(k)=B(k)-temp*C(k-l);U(k)=U(k)-temp*U(k4);EndM(N|=U(N-1)/B(N-1);fork=N-2:4:1M(k+l)=(U(k)-C(k)*M(k+2)/B(k);EndM(l)=3*(D(l)-dxO)/H(l)-M(2)/2;M(N+l)=3*(dxn-D(N)/H(N)-M(N)/2;fork=O:N4S(k+l,l)=(M(k+2)-M(k+l)/(6*H(k+l);

7、S(k+l/2)=M(k+l)/2;S(k+13)=D(k+l)-H(k+l)*(2*M(k+l)+M(k+2)/6;S(k+l/4)=Y(k+l);endcsfit2:第二類邊界條件functionS=csfit2(X,Y,dxOrdxn)N=length(X)-1;H=diff(X);D=diff(Y)./H;A=H(2:N-1);B=2*(H(1:N-1)+H(2:N);C=H(2:N);C=H(2:N)；U=6*diff(D);U(l)=U(l)-dxO;U(N-1)=U(N-l)-dxn;fork=2:N-ltemp=A(k-1)/B(k-1);B(k)=B(k)-temp*C(k-

8、1);U(k)=U(k)-temp*U(k-1);EndM(N)=U(N-D/B(N-l);fork=N-2:-1:1M(kfl)=(U(k)-C(k)*M(k+2)/B(k);EndM(l)=dxO;M(N+l)=dxn;fork=0:N-lS(k+1,1)=(M(k+2)-M(k+1)/(6*H(k+1);S(kflr2)=M(k+l)/2;S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2)/6;S(k+1,4)=Y(k+l);end畫圖：xl=0:.01:l;yl=polyval(Sl(l,:),xl-X(l);x2=l：,01:2;y2=polyval(Sl

9、(2,:),x2-X(2);x3=2:.01:33=polyval(Sl(3,:)/x3-X(3);x4=3:.01:44=polyval(Sl(4,:),x4-X(4);x5=4:.01:55=polyval(Sl(5,:)/x5-X(5);x6=5:.01:66=polyval(Sl(6,:),x6-X(6);»plot(xLyLx2,y2,x3,y3,x4,y4,x5,y5,x6,y6,XY'.')»holdon»xl=0:.01:171=polyval(S2(l/:)/xl-X(l);x2=l:.01:22=polyval(S2(2,:),

10、x2-X(2);x3=2:.01:33=polyval(S2(3/:)/x3-X(3);x4=3:.01:44=polyval(S2(4,:),x4-X(4);x5=4:.01:55=polyval(S2(5,:),x5-X(5);x6=5:.01:66=polyval(S2(6/:),x6-X(6);»plot(xLyLx2,y2,x3,y3,x4,y4,x5,y5,x6,y6,XY'.')»holdon»xl=0:.01:171=polyval(S3(l/:)/xl-X(l);x2=l:.01:2;y2=polyval(S3(2,:),x2-X

11、(2);x3=2:.01:3;y3=polyval(S3(3,:),x3-X)；x4=3:.01:44=polyval(S3(4,:)/x4-X(4);x5=4:.01:5,y5=polyval(S3(5,:),x5-X(5);x6=5:.01:66=polyval(S3(6,:)/x6-X(6);»plot(xLyLx2,y2,x3,y3,x4,y4,x5,y5,x6,y6,XY'.')!1!實驗結(jié)果第一類邊界條件：»K=0123456：Y=1001221：»dx0=-0.6:dxn=-1.8;»Sl=csfit1(X,Y,dxO,dx

12、n)SI=0.6400-1.0400-0.60001.0000-0.12000.8800-0.76000-0.16000.52000.64000-0.24000.04001.20001.0000Q.1200-0.6UUUU.bbUO2.UU0U-0.2400一0.3200-0.44002.0000第二類邊界條件:>>ditO=l;dxn二一1;>>S2=csfit2(IjY?<1x0,±cn)S2=-0.01110.5000-1.48891.00000.05560.4667-0.52220-0.21110.63330.57780-0.2111Q1.2111l.OQOQ0.0556-0.63330.57782.0000-0.0111-0.4667-0.52222.0000自然邊界條件：»dx0=0:dxn=0:»S3