自编的三次样条插值matlab程序(含多种边界条件)

数值计算第二次大作业——验证三次样条函数插值是否有几何不变性（1）给定的插值条件如下：i01234567Xi8.1258.49.09.4859.69.95910.1710.2Yi0.07740.0990.280.600.7081.2001.8002.177端点边界条件为第一类边界条件（给定一阶导数）：'0'7Y0.01087Y100三次样条函数的构造过程如下：设1231nnxxxxx共n个插值节点，则经过数据点1122,,,,,,nnxyxyxy的三次样条Sx是一组三次多项式：11231111111112232222222223321111111,,,,,,,,nnnnnnnnnnSxabxxcxxdxxxxxSxabxxcxxdxxxxxSxabxxcxxdxxxxx(1.1)由节点处的连续性可知：112321121121121231111111,,1,2,1.,1,2,1.,iiiiiiiinnnnnnnnnnnSxySxyinayinyybxxcxxdxxyybxxcxxdxx(1.2)由节点处的一阶与二阶光滑性可知：''112''1222112112122''21112212212112221121221112,,1,2,,.0230230262026iiiiiiiinnnnnnnnnnnnnnnnnSxSxSxSxinSxSxbcxxdxxbSxSxbcxxdxxbSxSxcdxxcSxSxc21212nnnndxxc(1.3)又设12nnncSx，记11,,1,2,,1iiiiiixxyyin，则由(1.3)可得：1,1,2,,1.3iiiiccdin(1.4)从（1.2）解得：212,1,2,,1.3iiiiiiiiiiiibcdccin(1.5)将(1.4)与(1.5)代入(1.3)得：21111222321122221111223,23nnnnnnnnnnncccccc(1.6)增加两个端点边界条件，因为11112,2nnncSxcSx，故有：1.第零类边界条件：自然样条，10,0.ncc2.第一类边界条件：给定端点一阶导数，设''1111,nnnSxvSxv，则有：1112111111112323nnnnnnnccvccv3.第二类边界条件：给定端点二阶导数，设1111,nnnSxvSxv，则有：1122nncvcv结合(1.6)及所给的边界条件即可解出ic,而11,233iiiiiiiiiiccdbcc，故可得到最终各个子区间上的三次样条函数。根据以上过程进行matlab编程，编写三次样条spline3函数，具体见附录。因此所编函数可对第一题求解：clear;clc;formatshortg;x1=[8.1258.499.4859.69.95910.16610.2]';y1=[0.07740.0990.280.60.7081.21.82.177]';u1=0.01087;un=100;xx1=[x1(1):0.001:x1(end)]';[yy1b1c1d1]=spline3(x1,y1,xx1,1,u1,un);fprintf('\t\tb1\t\tc1\t\td1\n');disp([b1c1(1:end-1,1)d1]);plot(x1,y1,'bo',xx1,yy1,'--k');gridon;b1c1d10.010870.144890.3680.174050.4485-0.3930.2878-0.258912.11531.52942.8188-69.141-0.56548-21.03573.61412.79458.247-512.32-28.949-259.94227988.599.51010.500.511.522.5388.599.51010.500.511.522.53图表一、三次样条曲线（2）坐标轴逆时针旋转45度，相当于节点顺时针旋转45度。设,Txy为旋转前的坐标，'',Txy为旋转后的坐标，则有：'cossin'sincosxxyy故旋转后的节点坐标为：theta=-pi/4;fori=1:length(x1)x2(i)=cos(theta)*x1(i)-sin(theta)*y1(i);y2(i)=sin(theta)*x1(i)+cos(theta)*y1(i);endfprintf('\t\t\tx2\t\t\ty2\n');disp([x2'y2']);x2y25.8-5.69056.0097-5.86976.562-6.1667.1312-6.28267.2889-6.28767.8906-6.19358.4612-5.91578.7519-5.6731端点处的一阶导数为：v1=(u1+tan(theta))/(1-u1*tan(theta));vn=(un+tan(theta))/(1-un*tan(theta));fprintf('\t\t\tv1\t\t\tvn\n');disp([v1vn]);v1vn-0.978490.9802则旋转后的三次样条的系数及图像为：xx2=[x2(1):0.001:x2(end)]';[yy2b2c2d2]=spline3(x2,y2,xx2,1,v1,vn);fprintf('\t\t\tb2\t\t\tc2\t\t\td2\n');disp([b2c2(1:end-1,1)d2]);plot(x2,y2,'*b',xx2,yy2,'-.k');gridon;b2c2d2-0.978490.67221-0.38277-0.747040.43138-0.090754-0.353620.28102-0.034909-0.0676290.221410.0533380.00617470.246640.00468970.30810.25510.102330.69920.430280.121955.566.577.588.59-6.4-6.3-6.2-6.1-6-5.9-5.8-5.7-5.65.566.577.588.59-6.4-6.3-6.2-6.1-6-5.9-5.8-5.7-5.6图表二、旋转后的三次样条曲线（3）将第（1）题中所得的样条曲线整体旋转-45度（即顺时针旋转45度），并与第（2）题的曲线画在同一幅图上，得：fori=1:length(xx1)xx3(i)=cos(theta)*xx1(i)-sin(theta)*yy1(i);yy3(i)=sin(theta)*xx1(i)+cos(theta)*yy1(i);endplot(xx3,yy3,'--',xx2,yy2,'-',x2,y2,'ok');gridon;legend('旋转前','旋转后');5.566.577.588.599.5-6.8-6.6-6.4-6.2-6-5.8-5.6-5.4-5.2旋转前旋转后5.566.577.588.599.5-6.8-6.6-6.4-6.2-6-5.8-5.6-5.4-5.2旋转前旋转后图表三、旋转前后样条曲线几何比较比较上图中的两条曲线，易知曲线不重合，故有以下结论：三三三次次次样样样条条条函函函数数数插插插值值值不不不具具具备备备几几几何何何不不不变变变性性性！！！附附附录录录三次样条插值函数matlab的m文件程序：function[yybcd]=spline3(x,y,xx,flag,vl,vr)%三次样条插值函数%(x,y)为插值节点，xx为插值点；%flag表端点边界条件类型：%flag=0:自然样条(端点二阶导数为0)；%flag=1:第一类边界条件(端点一阶导数给定)；%flag=2:第二类边界条件(端点二阶导数给定)；%vl,vr表左右端点处的在边界条件值。%样条函数为：Si(x)=yi+bi*(x-xi)+ci*(x-xi)^2+di*(x-xi)^3%b,c,d分别为各子区间上的系数值%yy表插值点处的函数值.iflength(x)~=length(y)error('输入数据应成对!');endn=length(x);a=zeros(n-1,1);b=a;d=a;dx=a;dy=a;A=zeros(n);B=zeros(n,1);fori=1:n-1a(i)=y(i);dx(i)=x(i+1)-x(i);dy(i)=y(i+1)-y(i);endfori=2:n-1A(i,i-1)=dx(i-1);A(i,i)=2*(dx(i-1)+dx(i));A(i,i+1)=dx(i);B(i,1)=3*(dy(i)/dx(i)-dy(i-1)/dx(i-1));end%自然样条端点条件(端点二阶导数为零)%ifflag==0A(1,1)=1;A(n,n)=1;end%---------------------------------%%端点一阶导数条件%ifflag==1A(1,1)=2*dx(1);A(1,2)=dx(1);A(n,n-1)=dx(n-1);A(n,n)=2*dx(n-1);B(1,1)=3*(dy(1)/dx(1)-vl);B(n,1)=3*(vr-dy(n-1)/dx(n-1));end%---------------%%端点二阶导数条件%ifflag==2A(1,1)=2;A(n,n)=2;B(1,1)=vl;B(n,1)=vr;end%---------------%c=A\B;fori=1:n-1d(i)=(c(i+1)-c(i))/(3*dx(i));b(i)=dy(i)/dx(i)-dx(i)*(2*c(i)+c(i+1))/3;end[mmnn]=size(xx);yy=zeros(mm,nn);fori=1:mm*nnforii=1:n-1ifxx(i)=x(ii)&xx(i)x(ii+1)j=ii;break;elseifxx(i)==x(n)j=n-1;endendyy(i)=a(j)+b(j)*(xx(i)-x(j))+c(j)*(xx(i)-x(j))^2+d(j)*(xx(i)-x(j))^3;endend