东南大学_数值分析_第四章_三次样条插值

东南大学《数值分析》上机练习——算法与程序设计实验报告第四章多项式插值与函数最佳逼近——曲线拟合之3次样条插值*****（学号）*****（姓名）上机题目要求见教材P195,37题。一、算法原理题目要求编写第一边界条件的3次样条插值函数的通用程序，同时根据汽车门曲线值点构造三次紧压样条曲线函数()Sx。其基本原理如下定义设0{(,)}Nkkkxy有N+1个点，其中01Naxxxb。如果存在N个三次多项式()kSx，系数为,0,1,2,3,,kkkkSSSS和满足如下性质：23,0,1,2,3()()()()kkkkkkkkSxssxxsxxsxx(1)111''111''''111(),0,1,...,()(),0,1,...,2()(),0,1,...,2()(),0,1,...,2kkkkkkkkkkkkkkSxyforkNSxSxforkNSxSxforkNSxSxforkN(2)则成()Sx为三次样条函数。现证明其存在：由于()Sx是分段三次多项式，其二阶导数是在区间0[,]Nxx内是分段线性的。根据线性拉格朗日插值()()kSxSx可以表示为：11111()()()kkkkkkkkkxxxxSxSxSxxxxx(3)用111(),()kkkkkkkmSxmSxxx和h代入上式，得11()kkkkkkkmmSxxxxxhh(4)将上式积分两次，会引入两个积分常数，可得到如下形式：33111()()()66kkkkkkkkkkkmmSxxxxxpxxqxxhh(5)第二章非线性方程的解法将1kkxx和代入上式，并利用11(),()kkkkkkySxySx可得两个方程：2211;66kkkkkkkkkkmmyhphyhqh(6)求解,kkpq，并将所得的结果带入方程(5)得3311111()66()()66kkkkkkkkkkkkkkkkkmmSxxxxxhhymhymhxxxxhh(7)求式（7）的导数，并化简得1111112()6(),1,...,1,()/kkkkkkkkkkkkkkkhmhhmhmuuddkNdyyh其中(8)由上述方程可得如下方程011121002()hhmhmuhm(9)11112()kkkkkkkkhmhhmhmu(10)12211112()NNNNNNNNhmhhmuhm(11)重组上述方程，得三角线性方程组HMV，表示为11111222223223222111NNNNNNNNNmvbcabcmvabmvabcabmv(12)该式具有严格对角优势。算出系数{}km后可由如下公式计算()kSx的样条系数,{}kjs。1,0,11,2,3(2),,6,26kkkkkkkkkkkkkhmmsysdmmmssh(13)东南大学《数值分析》上机练习——算法与程序设计实验报告二、流程图用已知的N+1个点0{(,)}Nkkkxy构造三次紧压样条曲线的问题。其通用程序流程图1所示。图1关于求解三阶样条曲线算法流程图具体步骤如下：1)计算1kkkxxh，1()/kkkkkdyyh，16()kkkudd2)联系端点约束条件求解样条函数的系数。三、计算代码核心代码fork=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);endfork=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2))/B(k);fork=0:N-1S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));开始输入数据及边界条件计算H，D，U利用端点约束计算求各段的样条系数结束第二章非线性方程的解法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完整代码functionS=Three_fit(X,Y,dx0,dxn)%Input-Xisavectorthatcontainsalistofabscissas%-Yisavectorthatcontainsalistofordinates%-dx0=S'(x0)firstderivativeboundarycondition%-dxn=S'(xn)firstderivativeboundarycondition%Output-SrowsofSarethecoefficients,indescendingorder,for%thecubicinterpolantsN=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);U=6*diff(D);%clampedsplineendpointconstraintsB(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));%Gausseliminationtosolvemkfork=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);fork=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;fork=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东南大学《数值分析》上机练习——算法与程序设计实验报告四、计算结果及分析根据题目数据调用函数如下X=0:10;Y=[2.51,3.30,4.04,4.70,5.22,5.54,5.78,5.40,5.57,5.70,5.80];dx0=0.8;dxn=0.2;S=Three_fit(X,Y,dx0,dxn)plot(X,Y,'o'),holdonyprint=[];forii=1:10x=X(ii):0.1:X(ii+1);y=polyval(S(ii,:),x-X(ii));%分段拟合yprint=[yprint,polyval(S(ii,:),0.5)];plot(x,y)endxlabel('x','FontSize',20),ylabel('y','FontSize',20)title('三次样条插值','FontSize',22)holdoffdisp('-------------------------------------------------------------------------------')disp('x(i)=0.51.52.53.54.55.56.57.58.59.5')disp('-------------------------------------------------------------------------------')disp(['S(i+0.5)='num2str(yprint(1))''num2str(yprint(2))''...num2str(yprint(3))''num2str(yprint(4))''num2str(yprint(5))...''num2str(yprint(6))''num2str(yprint(7))''num2str(yprint(8))...''num2str(yprint(9))''num2str(yprint(10))])disp('-------------------------------------------------------------------------------')输出三次样条函数的系数S=-0.0085-0.00150.80002.5100-0.0045-0.02700.77153.3000-0.0037-0.04040.70414.0400-0.0409-0.05140.61234.70000.1074-0.17410.38685.2200-0.26850.14790.36065.54000.4266-0.6575-0.14915.7800-0.26790.6222-0.18445.40000.0549-0.18140.25655.57000.0584-0.01680.05845.7000第二章非线性方程的解法打印出(+0.5)0,1,,9Sii，的值-------------------------------------------------------------------------------x(i)=0.51.52.53.54.55.56.57.58.59.5-------------------------------------------------------------------------------S(i+0.5)=2.90863.67844.38154.98825.38335.72375.59445.42995.65985.7323-------------------------------------------------------------------------------输出拟合曲线:0123456789102.533.544.555.56xy三次样条插值X-YS1(x)S3(x)S4(x)S5(x)S6(x)S7(x)S8(x)S9(x)S10(x)S11(x)Fig2三次样条曲线分析：由图形可以看出，三阶样条曲线非常光滑，用该法能获得很好的拟合效果。五、结论对数据进行多项式拟合在CAD，CAM和计算机图形处理软件中有许多应用。我们常常需要画出经过多个点的无误差的光滑的曲线。从数学的角度上看就是要求图像及其一阶，二阶导数在分段区间内连续。本次报告讨论了有关三阶样条函数的原理及应用，该法可使函数图像二阶连续可导，反映出曲线的光滑。该法具有重要的工程应用意义。