三次样条插值实验报告专业班级学号姓名一、 实验内容和要求1、阅读上面的文字和程序，试运行，检验程序和上面叙述的正确性。

2、阅读上面的 MATLAB 程序；查资料，了解各 MATLAB 语句及命令。

3、画程序流程图，理解并描述算法。

4、修改上面的程序，能根据给定数据点，求（1）自然样条插值，边界 S ′′(a )=0, S ′′(b )=0 （2）第二类边界条件，S ′′(a )和S ′′(b )是确定的。

5、使用上面的程序，根据数据点（0，1），（1，0），（2，0），（3，1），（4，2）， （5，2），和（6，1），求三种不同的三次样条插值，其中S ′(0)=−0.6，S ′(6)=-1.8，S ′′(0)=1，S ′′(6)=−1；S ′′(0)，S ′′(6)=0.在同一坐标系中，画出这 3 个三次样条插值和这些数据点。

6、写实验报告（实验内容+算法描述+程序+写成分段函数的结果描述+截图）。

二、 算法说明定义：设有N+1个点 ，其中 。

如果存在N 个三次多项式 ，系数为 ，满足如下性质： （1）（2）01n a x x x b =<<<=()k S x 1[,],0,1,,1k k x x x k N +∈=-0[,]Nk k k x y =,0,1,2,3(),(),(),()k k k k S x S x S x S x 23,0,1,2,3()S ()()()()k k k k k k k k S x x s s x x s x x s x x ==+-+-+-()k k S x y =0,1,,k N =（3）则称函数 为三次样条函数。

因为是分段三次多项式，所以 在区间 上是分段线性的。

（1）用 代入上式，可得（2）将（2）式积分两次，会引入两个积分常数，并得到（3） 将 代入上（3），并使用 ，可分别得到包含的方程： （4）求解这两个方程，求出 ，而且将这些值代入方程（3）中，可得到三次多项式方程：（5）表达式5可以简化成只包含未知系数的形式。

为求这些值，必须使用（5）式的导数，即（6） 1lim ()lim ()kkk k x x x x S x S x +--→→=1,2,,1k n =-1lim ()lim ()kkkk x x x x S x S x +--→→''=1,2,,1k n =-1lim ()lim ()kkkk x x x x S x S x +--→→''''=1,2,,1k n =-()S x 1111()S (x )S (x )k kkk k k k k k x x x x S x x x x x ++++--''''''=+--()S x ()S x ''[,]a b 111(),()kkk k kk kM S x M S x h x x +++''''===-和11()M k kkk k k k x x x x S x M h h ++--''=+1,0,1,, 1.k k x x x k N +≤≤=-33111()()()()()66k k k k k k k k k kkMM S x x x x x p x x q x x h h +++=-+-+-+-1,k k x x +和+1+1()()k k k k k k S x y S x y ==和k k p q 和2211,66k k k k k k k k k k M M y h p h y h q h ++=+=+k k p q 和3311111()()()()()6666k k k k k k k kk k k k k k k kkM M y M hy M h S x x x x x x x x x h h h h +++++⎛⎫⎛⎫=--+-+--+-- ⎪⎪⎝⎭⎝⎭{}k M 221111()()()226k k k k k k k kkk k k k k k k M M y M h y m hS x x x x x h h h h h ++++⎛⎫'=--+---+- ⎪⎝⎭在处计算（6），并简化结果可得到 ，其中 （7） 同理，在式（6）中用 并计算在 处的解，可得 （8）利用节点处一阶导函数连续及方程（6）、（7），可得到包含的重要关系式（9）其中方程组（9）中的未知数是要求的值 ，其他的项可通过数据点集 进行简单的计算得到的常量。

因此方程（9）是包含N+1个未知数，具有N-1个线性方程的不定方程组。

所以需要另外两个方程组才能求解，即边界条件。

如果已知，则 100001113(())23(())2N N N N N MM d S x h M M S x d h ---'=--'=--（10）（11） 根据（9）（10）求出 后，可利用下面的公式计算 的样条系数 。

k x 1()36k k k k k k k M m S x h h d +'=--+1k k k k y y d h +-=11,(),k k k S x -'-代替可得到k x 11111()36k k k k k k k Mm S x h h d -----'=++11,,k k k M M M -+11112()k k k k k k k k h M h h M h M μ---++++=16(),1,2,, 1.k k k d d k N μ-=-=-{}0Nk M {}0(,)Nk k k x y =0(),()n S x S x ''{}0Nk M ()k S x 111112222322222111N N N N N N N N N b c M a b c M a b C M a b M μμμμ---------⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦,0k k s y =，1,1(2)6k k k k k h M M s d ++=-，,22k k Ms =，1,36k k k k M M s h +-=，1k k k k y y d h +-=，1k k k h x x +=-三、 源程序csfit1:第一类边界条件function S=csfit1(X,Y,dx0,dxn)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);B(1)=B(1)-H(1)/2;U(1)=U(1)-3*(D(1)-dx0); B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(dxn-D(N)); for k=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); EndM(N)=U(N-1)/B(N-1); for k=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2))/B(k); EndM(1)=3*(D(1)-dx0)/H(1)-M(2)/2; M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;{},k jsfor k=0:N-1S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));S(k+1,2)=M(k+1)/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+1);endcsfit2:第二类边界条件function S=csfit2(X,Y,dx0,dxn)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(1)=U(1)-dx0;U(N-1)=U(N-1)-dxn;for k=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);EndM(N)=U(N-1)/B(N-1);for k=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2))/B(k);EndM(1)=dx0;M(N+1)=dxn;for k=0:N-1S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));S(k+1,2)=M(k+1)/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+1);end画图：x1=0:.01:1;y1=polyval(S1(1,:),x1-X(1));x2=1:.01:2;y2=polyval(S1(2,:),x2-X(2));x3=2:.01:3;y3=polyval(S1(3,:),x3-X(3));x4=3:.01:4;y4=polyval(S1(4,:),x4-X(4));x5=4:.01:5;y5=polyval(S1(5,:),x5-X(5));x6=5:.01:6;y6=polyval(S1(6,:),x6-X(6));>> plot(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,X,Y,'.') >> hold on>> x1=0:.01:1;y1=polyval(S2(1,:),x1-X(1));x2=1:.01:2;y2=polyval(S2(2,:),x2-X(2));x3=2:.01:3;y3=polyval(S2(3,:),x3-X(3));x4=3:.01:4;y4=polyval(S2(4,:),x4-X(4));x5=4:.01:5;y5=polyval(S2(5,:),x5-X(5));x6=5:.01:6;y6=polyval(S2(6,:),x6-X(6));>> plot(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,X,Y,'.') >> hold on>> x1=0:.01:1;y1=polyval(S3(1,:),x1-X(1));x2=1:.01:2;y2=polyval(S3(2,:),x2-X(2));x3=2:.01:3;y3=polyval(S3(3,:),x3-X(3));x4=3:.01:4;y4=polyval(S3(4,:),x4-X(4));x5=4:.01:5;y5=polyval(S3(5,:),x5-X(5));x6=5:.01:6;y6=polyval(S3(6,:),x6-X(6));>> plot(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,X,Y,'.')四、实验结果第一类边界条件：第二类边界条件：自然边界条件：图像：五、说明与分析实验结果显示，对同一组数据，不同边界条件，会导致形成曲线有一定差距，但总体趋势不变。