數(shù)值分析（宋）第2次大作業(yè)三次樣條函數(shù)插值

1/5三次樣條函數(shù)插值（數(shù)值分析第二次大作業(yè)）姓名：**學(xué)號：**班級：**1．給定插值條件如表1所示：表1.插值條件012345678.1258.49.09.4859.69.95910.16610.20.07740.0990.2800.600.7081.2001.8002.177作三次樣條函數(shù)插值，取第一類邊界條件：,。畫插值曲線圖像。解：三次樣條插值實際上是由分段三次曲線并接而成，在連接點即樣點上要求二階導(dǎo)數(shù)連續(xù)。在本題中，共給出了個插值節(jié)點，故需要確定出7個小區(qū)間上的7條三次多項式曲線。設(shè)，，，根據(jù)三次樣條插值理論，對于第一類邊界條件，若令，，，，則三次樣條表達(dá)式為(1)其中，且為下列矩陣方程的解,記為根據(jù)以上各式，編寫Matlab程序（程序詳見附錄1）計算得到：其中是由組成的向量。將相應(yīng)的和值代入三次樣條表達(dá)式(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．逆時針旋轉(zhuǎn)坐標(biāo)軸45o保持（1．）中結(jié)點和邊界條件的幾何關(guān)系不變，再次作三次樣條函數(shù)插值，畫出插值曲線的圖像。解：首先編寫Matlab程序（詳見程序附錄2），實現(xiàn)逆時針旋轉(zhuǎn)坐標(biāo)軸45o并保持（1．）中結(jié)點和邊界條件的幾何關(guān)系不變，得到旋轉(zhuǎn)后的節(jié)點坐標(biāo)和邊界條件，如表2所示。表2.逆時針旋轉(zhuǎn)坐標(biāo)軸后的節(jié)點坐標(biāo)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而邊界條件變?yōu)?。然后進(jìn)行三線樣條插值，畫出插值曲線如圖2所示。圖2旋轉(zhuǎn)坐標(biāo)軸后的三次樣條插值曲線3．比較（1．）、（2．）的結(jié)果，能得到什么結(jié)論？答：比較圖1和圖2可知，旋轉(zhuǎn)前的三次樣條插值曲線與旋轉(zhuǎn)后的三次樣條插值曲線在形狀上已經(jīng)完全不一樣，說明三次樣條插值不具有幾何不變性。

程序附錄附錄1旋轉(zhuǎn)前的三次樣條插值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);%進(jìn)行三次樣條插值的M函數(shù)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旋轉(zhuǎn)前的三次樣條插值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);%旋轉(zhuǎn)節(jié)點的M函數(shù)functiony=axisrotation(x,theta)theta=theta*pi/180;A=[cos(theta)sin(theta);-sin(theta)cos(theta)];x=x';y=A*x;%旋轉(zhuǎn)導(dǎo)數(shù)