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

1、数值计算第二次大作业验证三次样条函数插值是否有几何不变性（1）给定的插值条件如下:101234567Xt8.1258.49.09.4859.69.95910.1710.2£0.07740.0990.280.600.7081.2001.8002.177端点边界条件为第一类边界条件（给定一阶导教）:Yo = 0.01087y7 = 100三次样条函数的构造过程如下：设竝共"个插值节点，则经过数据点（心）,（忑，'2），（耳，儿）的三次样条S（x）是一组三次多项式：S （x） = q + 勺（x-xj + q （x-+ d （x-xj5,x g xnxj, （X） = 2

2、 +fe2（X"X2）+C2（X-J2）2 + 2（X" J2）5>XeX2>Jj（11）S”T （x） = % +（X -V1）+ c“T （X - .-1 ）2 + <-1 （X - x J'，x w 卜 I,兀由节点处的连续性可知：£（兀）=y,S：（和）=y屮，心 1,2,“-1.q = y,i = l,2, I.一儿=必+人（ 一再）+ 9（农一兀）'+ 4（农一呂）'，（1.2）由节点处的一阶与二阶光滑性可知:S,t (兀)=£ (xJ,S；t (兀)=5；(X.),/ = 1,2,/.0=Sl(x2

3、)-S2(x2) = bl + 2cl(x2-xl)+3dl(x2-xl)2-b2f o = S；T (儿.J 一 S；T ( El J =+ 2c“ (兀-1 一爲L J + 迅 T (Xi - x“2)' - g (I 3)0 = S；(X?) - S；(X?) = 2q + 6£ (x2 -xj- 2c20 = S；L2 () - S；I (x”. J = 2cn_2 + 6d”T (兀T -兀L2 ) - K-l(A ASt)/訂犬n-ln-2又设 C" = S；T (x” )/2 ,记 <5,=和-禺4 =畑->；，i = 1,2,上一1 ,则

4、由(1.3)可得:di =:c , / = 1,2,从(12)解得：bf = -(2 + c/41), i = 1,2,-将(1.4)与(1®代入(1.3)得:恥 + 2(4 + 6匕 + 6£ = 3 宁增加两个端点边界条件，因为2C = S；（xJ,2c; = S：t（£）,故有:1. 第零类边界条件：自然样条，G = 0,cfl = 0.2. 第一类边界条件：给定端点一阶导数，设St(兀)=*,则有：2恥+ <5心=3(织/4-叫)i 氏+ 2 玄品=3( V,- m3. 第二类边界条件：给定端点二阶导数，设S；(xJ =片,SJe)"”，则

5、有：c c5结合(1.6)及所给的边界条件即可解出q，而 < =弋0 =寸(2q + %),故5oioi 5可得到最终各个子区间上的三次样条函数.根携以上过程进行matlab编写三次样条spline3函数，具体见附录.因此所编函数可对第一题求解： clear;clc;format short g;xl=8.125 8.4 9 9.485 9.6 9.959 10.166 10.21;yl二0.0774 0.099 0.28 0.6 0.708 1.2 1.8 2.177;ul = 0 01087;un=l00;xxl=xl(1):0001:xl(end);yyl bl cl dl=spl

6、ine3(xl,yl,xxl,l,ul,un);fprintf(ttblttclttdln1);disp(bl cl(1:end-1/1) dl);plot(xl ,yl, bo,xxlyylJ k1)；grid on;bl cl dl0.010870.144890.3680174050.4485-0.3930.2878-0.258912.11531.52942.818869.141-0.5654821.03573.61412.79458.247-51232-28.949-259.942279坐标轴逆时针旋转45度，相当于节点顺时针旋转45度设为旋转蔚的坐标，为旋转后的坐标，则有：卜甘。sO故

7、旋持后的节点坐标为：一 sin 0cose;theta=-pi/4;for i=l:length(xl)x2(i)=cos(theta)*xl(i)-sin(theta)*yl (i); y2(i)-sin(theta)*xl(i)+cos(theta)*yl (i); end fprintf(tttx2ttty2n)； disp(x2- y2 M)；x2y25.86.0097-56905-5.86976.562-6.1667.13127.28897.89068.4612-62826 -6.2876 一61935一591578.75195.6731端点处的一阶导数为：vl=(ul+tan(th

8、eta)/(l-ul*tan(theta); vn=(un+tan(theta)/(l-un*tan(theta); fprintf(tttvltttvnn1);disp(vl vn)；vlvn-0.978490.9802则旋转后的三次样条的系数及图像为：xx2=x2(1) :000l:x2(©nd)yy2 b2 c2 d2=spline3(x2z y2,xx2z1,vlz vn);fprintf(tttb2tttc2tttd2n1); disp(b2 c2(l:end-lrl) d2);plot(x2ry2rb'/xx2/yy2/-.k-);grid on;b2c2d2-0

9、.978490.67221-0.38277-0.747040431380.090754-0.353620.28102-0.0349090.0676290.221410.0533380.00617470.246640.004689703081025510102330.69920.430280.12195(3)将第(1)题中所得的样条曲线搭体旋转45度(即顺时针旋转45度)，并与 第(2)題的曲线画在同一幅图上，得：for i=l:length(xxl)xx3 (i)=cos(theta)*xxl(i)-sin(theta)*yyl(i);yy3(i)=sin(theta)*xxl(i)+cos(

10、theta)*yyl(i);endplot(xx3,yy3,.xx2,yy2，',x2,y2,ok');grid on; legend (1旋转前1 J旋转后1 )；图裘三、旋转萌后样条曲线几何比校 比校上图中的两条曲线，易知曲线不重合，故有以下结论：三次样条爲裁插值不具备几何不变性|三次样条插值函Mmatlab的m文件程序： function yy b c d=spline3(x/y,xx/flag,vlzvr) %三次样条插值因数吻区y)为插值节点，xx为插值点；% fl昭表端点边界条件类型：呀 flag=O :自然样条(端点二阶导数为0);吩 flag=l:第一类边界条件

11、(端点一阶导数给定)；吩 Oag=2 :第二类边界条件(端点二阶导数给定)； 呀vl,vr左右端点处的在边界条件值.% 样条函数为：Si(x)二yi+bir*xi)+ci*(xxi)A2+diqxxi)八3吩b,c,d分别为各子区间上的系数值吩yy裘插值点处的函数值if length(x)=length(y)error('输入数据应成对!');endn=length(x);a=zeros(n-l,l);b=a;d=a;dx=a;dy=a;A=zeros(n);B=zeros(n,l)；for i=l:n-la(i)=y(i)；dx(i)=x(i+l)-x(i);dy(i)=y(

12、i+l)-y(i);endfor i=2:n-lA(i,i-l)=dx(i.l);A(i/i)=2*(dx(M)+dx(i);A(izi+l)=dx(i);B(i4)=3*(dy(i)/dx(i)-dy(i-l)/dx(i-l);end%自然样条端点条件(端点、二阶导数为辈)if flag=0A(l,l)=l;A(n,n)=l;end吻吻呀端点一阶导数条伴if flag=lA(lzl)=2*dx(l);A(l/2)=dx(l);A(n,n-l)=dx(n-l);A(nzn)=2*dx(n-l);B(l,l)=3*(dy(l)/dx(l).vl);B(n,l)=3*(vr-dy( n-1 )/dx( n-1);end%呀%端点二阶导数条件if flag=2A(l,l)=2p(n,n)=2;B(lzl)=vl;B(nzl)=vr;end96% c=AB;for i=l:n-ld(i)=(c(i+l)c(i)/(3*dx(i);b(i)=dy(i)/dx(i).dx(i)*(2*c(i)+c(i+l)Aendmm nn=size