数值分析（宋）第2次大作业三次样条函数插值

1/5三次样条函数插值（数值分析第二次大作业）姓名：**学号：**班级：**1．给定插值条件如表1所示：表1.插值条件012345678.1258.49.09.4859.69.95910.16610.20.07740.0990.2800.600.7081.2001.8002.177作三次样条函数插值，取第一类边界条件：,。画插值曲线图像。解：三次样条插值实际上是由分段三次曲线并接而成，在连接点即样点上要求二阶导数连续。在本题中，共给出了个插值节点，故需要确定出7个小区间上的7条三次多项式曲线。设，，，根据三次样条插值理论，对于第一类边界条件，若令，，，，则三次样条表达式为(1)其中，且为下列矩阵方程的解,记为根据以上各式，编写Matlab程序（程序详见附录1）计算得到：其中是由组成的向量。将相应的和值代入三次样条表达式(1)即可得到所求三次样条插值曲线：S1(x)=0.175626*(8.4-x)^3+0.543630*(x-8.125)^3+0.050715*x-0.338313S2(x)=0.249164*(9-x)^3-0.143841*(x-8.4)^3+0.443148*x-3.67726S3(x)=-0.177947*(9.485-x)^3+1.93733*(x-9)^3+0.162228*x-1.15975S4(x)=8.17048*(9.6-x)^3-60.9703*(x-9.485)^3+1.85352*x-16.993S5(x)=-19.5309*(9.959-x)^3+54.0830*(x-9.6)^3-8.11695*x+79.5344S6(x)=93.7961*(10.166-x)^3-418.526*(x-9.959)^3+24.8511*x-247.124S7(x)=-2548.09*(10.2-x)^3+39730.7*(x-10.166)^3-37.786*x+386.033画出插值曲线图像如图1所示。图1三次样条插值曲线

2．逆时针旋转坐标轴45o保持（1．）中结点和边界条件的几何关系不变，再次作三次样条函数插值，画出插值曲线的图像。解：首先编写Matlab程序（详见程序附录2），实现逆时针旋转坐标轴45o并保持（1．）中结点和边界条件的几何关系不变，得到旋转后的节点坐标和边界条件，如表2所示。表2.逆时针旋转坐标轴后的节点坐标012345675.86.00976.5627.13127.28897.89068.46128.7519-5.6905-5.8697-6.166-6.2826-6.2876-6.1935-5.9157-5.6731而边界条件变为,。然后进行三线样条插值，画出插值曲线如图2所示。图2旋转坐标轴后的三次样条插值曲线3．比较（1．）、（2．）的结果，能得到什么结论？答：比较图1和图2可知，旋转前的三次样条插值曲线与旋转后的三次样条插值曲线在形状上已经完全不一样，说明三次样条插值不具有几何不变性。

程序附录附录1旋转前的三次样条插值clearall;X=[8.1258.49.09.4859.69.95910.16610.2];Y=[0.07740.0990.2800.600.7081.2001.8002.177];dy0=0.01087;dyn=100;[m,H,lambda,mu,D,A,M]=cubicspline(X,Y,dy0,dyn);%进行三次样条插值的M函数function[m,H,lambda,mu,D,A,M]=cubicspline(X,Y,dy0,dyn)m=length(X);A=zeros(m,m);n=m-1;H=zeros(1,n);H(1)=X(2)-X(1);lambda=zeros(1,n);lambda(1)=1;mu=zeros(1,n);mu(n)=1;D=zeros(1,m);D(1)=6*((Y(2)-Y(1))/H(1)-dy0)/H(1);fork=1:nhk=X(k+1)-X(k);H(k)=hk;endfork=1:n-1mu(k)=H(k)/(H(k)+H(k+1));endfork=2:nlambda(k)=1-mu(k-1);endD(m)=6*(dyn-(Y(m)-Y(n))/H(n))/H(n);fork=2:n;D(k)=6*((Y(k+1)-Y(k))/H(k)-(Y(k)-Y(k-1))/H(k-1))/(H(k-1)+H(k));endfori=1:m-1A(i,i)=2;A(i,i+1)=lambda(i);A(i+1,i)=mu(i);endA(m,m)=2;M=A\D';holdon;plot(X,Y,'bo');fork=1:nx=X(k):0.001:X(k+1);y=zeros(1,length(x));fort=1:length(x)y(t)=M(k)*(X(k+1)-x(t))^3/(6*H(k))+M(k+1)*(x(t)-X(k))^3/(6*H(k))+Y(k)*(X(k+1)-x(t))/H(k)-M(k)*H(k)*(X(k+1)-x(t))/6+Y(k+1)*(x(t)-X(k))/H(k)-M(k+1)*H(k)*(x(t)-X(k))/6;endplot(x,y);endholdoff;附录2旋转前的三次样条插值clearall;X=[8.1258.49.09.4859.69.95910.16610.2];Y=[0.07740.0990.2800.600.7081.2001.8002.177];dy0=0.01087;dyn=100;theta=45;L=length(X);u=zeros(2,1);fori=1:Lu=axisrotation([X(i),Y(i)],theta);X(i)=u(1);Y(i)=u(2);endu=qiexianxuanzhuan([dy0dyn],theta);XYdy0=u(1)dyn=u(2)[m,H,lambda,mu,D,A,M]=cubicspline(X,Y,dy0,dyn);%旋转节点的M函数functiony=axisrotation(x,theta)theta=theta*pi/180;A=[cos(theta)sin(theta);-sin(theta)cos(theta)];x=x';y=A*x;%旋转导数