三次样条插值试验报告

1、三次样条插值实验报告专业班级学号姓名一、实验内容和要求1、阅读上面的文字和程序，试运行，检验程产和上面叙述的正确性。2、阅读上面的MATLAB程序：查资料，了解各MATLAB语句及命令。3、画程序流程图，理解并描述算法。4、修改上面的程序，能根据给定数据点，求(1)自然样条插值，边界S"(a)=O,S"(b)=O0回(2)第二类边界条件,S"(a)和S”(b)是确定的。5、使用上面的程序，根据数据点(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个三次样条插值和这些数据点。6、写实验报告(实验内容+算法描述+程序+写成分段函数的结果描述+截图)。二、算法说明定义：设有N+1个点，箕金;H如果存?8=%<%<,Y4=bN个三次多项式a(号数为，满是姆闵调心(x),8k3(x)(1)S(X)=Sk(x)=Sk0+8kl(X-5)+、式x-X1t尸+(X-%)3xw区,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则称函数为8次样条函数,因为是分段3；欠多项式，所以在区间上足分段线性的la,b(1)Pl+sd初队S”（）神户%-%(2)S；(x)=Mk+Mk+1<x<xk+pk=0,1,-sN-1.将（2）式积分两次，会引入两个积分常数，并得到（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可以简化成只包含未知系数的形题J为求这些值，必须使用（5）式的导数，即在处计算（6）,并简化结果可得到,其河新广孕L号4+44=巫p同理，在式（6）中用并计算消处的解战q（得（8）利用菰短鸣一导野诲绿监方程（6）、（7）,可得到包含Mi尚甯喙联系式（9）14TMl+2（%t+hk）Mk+hkMk+1=4其中74=6（dk-dk_i）,k=l,2,、N-L方程组（9）中的未知数是要求的值，馔篇的项可通过数据点集进行简里队冲硼到的常量。因此方程（9）是包含N+1

5、个未知数，具有N-1个线性方程的不定方程组。所以需要另外两个方程组才能求解，即边界条件。如果已知，财与）$（、）M。=（-SXq）-%2（10）hN-l2根据（耳£,410）求出后Mk时利用下面的公式计算的样条条幽鼠,。=yi，%=4-%Qm1Me）,产学“，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!实验结果第一类边界条件：»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