《计算方法大作业(第二次)》由会员分享，可在线阅读，更多相关《计算方法大作业(第二次)（7页珍藏版）》请在金锄头文库上搜索。
1、数值计算第二次大作业1给定插值条件如下：给定插值条件如下：i 0 1 2 3 4 5 6 7Xi 8.125 8.4 9.0 9.485 9.6 9.959 10.166 10.2Yi 0.0774 0.099 0.280 0.60 0.708 1.200 1.800 2.177作三次样条函数插值，取第一类边界条件作三次样条函数插值，取第一类边界条件Y0=0.01087 Y7=100根据题目要求，首先要构造三次样条函数，三次样条函数的构造过程如下：设有共n个插值节点，任意给定一组常数，要求构造一0x1xnx0y1yny个插值三次样条函数，使得如下插值条件得以满足：)(xS，i=0，1，n( )
2、iiS xy经过插值点的三次样条函数是一组三次多项式，即有：23 1111111112 23 222222222323 111111111( )()()() , , ( )()()() ,( )()()() ,nnnnnnnnnnS xab xxc xxd xxxx x Sxab xxc xxdxxxx xSxabxxcxxdxxxxx M由节点处的连续性可知：11( ),(),1,2,1.iiiiiiS xy S xyinL23 2112112112123 1111111,1,2,1, ()()()()()()iinnnnnnnnnnnay in yyb xxc xxd xxyybxxcxx
3、dxx LM由节点处的一阶与二阶光滑性可知： 11( )( ),( )( ),1,2,iiiiiiiiSxS xSxS xinL又设，记，则 1()/ 2nnncSx11,1,2,1iiiiiixxyy in L。再根据边界条件，从而可相继解出1,1,2,13ii i iccdinL,iib c用 matlab 编程，编写三次样条函数（见附录） ，对第一题求解： format short g; x1=8.125,8.4,9.0,9.485,9.6,9.959,10.166,10.2; y1=0.0774,0.099,0.280,0.60,0.708,1.200,1.800,2.177; u1=
4、0.01087;un=100; xx1=x1(1):0.001:x1(end); yy1 b1 c1 d1=spline3(x1,y1,xx1,1,u1,un); fprintf(ttb1ttc1ttd1n); b1c1d1 disp(b1 c1(1:end-1,1) d1);0.01087 0.14489 0.3680.17405 0.4485 -0.3930.2878 -0.25891 2.11531.5294 2.8188 -69.141-0.56548 -21.035 73.61412.794 58.247 -512.32-28.949 -259.9 42279 plot(x1,y1,
5、bo,xx1,yy1,r-); grid on 画出插值曲线的图像。画出插值曲线的图像。图 1 三次样条曲线2逆时针旋转座标轴逆时针旋转座标轴 45o 保持（保持（1 ）中结点和）中结点和边界条件的几何关系不变，再次作三边界条件的几何关系不变，再次作三 次样条函数插值，画出插值曲线的图像。次样条函数插值，画出插值曲线的图像。坐标轴逆时针旋转坐标轴逆时针旋转 45，相当于节点顺时针旋转，相当于节点顺时针旋转 45。设。设为旋转前的坐标，为旋转前的坐标，( , )Tx y为旋转后的坐标，则可以得到如下关系：为旋转后的坐标，则可以得到如下关系：( ,)Tx ycossinsincosxxyy 故旋转
6、后的节点坐标为：故旋转后的节点坐标为： theta=-pi/4; for i=1:length(x1) x2(i)=cos(theta)*x1(i)-sin(theta)*y1(i); y2(i)=sin(theta)*x1(i)+cos(theta)*y1(i); end fprintf(tttx2ttty2n); disp(x2 y2);5.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+ta
7、n(theta)/(1-u1*tan(theta); vn=(un+tan(theta)/(1-un*tan(theta); fprintf(tttv1tttvnn); v1vn disp(v1 vn);-0.97849 0.9802 则旋转后的三次样条的系数及图像为： xx2=x2(1):0.001:x2(end); yy2 b2 c2 d2=spline3(x2,y2,xx2,1,v1,vn); fprintf(tttb2tttc2tttd2n); b2c2d2 disp(b2 c2(1:end-1,1) d2);-0.97849 0.67221 -0.38277-0.74704 0.43
8、138 -0.090754-0.35362 0.28102 -0.034909-0.067629 0.22141 0.0533380.0061747 0.24664 0.00468970.3081 0.2551 0.102330.6992 0.43028 0.12195 plot(x2,y2,b+,xx2,yy2,m-.); grid on;图 2 旋转后的三次样条曲线 3比较（比较（1 ） 、 （2 ）的结果，能得到什么结论？）的结果，能得到什么结论？ 将（将（1 ）中所得的三次样条曲线整体顺时针旋转）中所得的三次样条曲线整体顺时针旋转 45，并与二题，并与二题（2 ）中的三次样）中的三次样
9、 条曲线画在同一幅图中比较，得条曲线画在同一幅图中比较，得 for i=1:length(xx1) xx3(i)=cos(theta)*xx1(i)-sin(theta)*yy1(i); yy3(i)=sin(theta)*xx1(i)+cos(theta)*yy1(i); end plot(x2,y2,bo,xx3,yy3,r-,xx2,yy2,m); grid on; legend(节点,旋转前,旋转后);图 3 旋转前后三次样条曲线几何比较 比较图中两条曲线可知，曲线不重合，故三次样条插值不具备几何不变性。附录：附录： 三次样条插值函数程序三次样条插值函数程序 functionyy,b,
10、c,d=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 表示插值点处的函数值 if length(x)=length(y)error(输入数据应成
11、对！); end n=length(x); a=zeros(n-1,1); b=a;d=a;dx=a;dy=a;A=zeros(n);B=zeros(n,1); for i=1:n-1a(i)=y(i);dx(i)=x(i+1)-x(i);dy(i)=y(i+1)-y(i); end for i=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 %自然样条端点条件（端点二阶导数为 0） if flag=0A(1,1)=1;A(n,n)=1
12、; end %端点一阶导数条件 if flag=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 %端点二阶导数条件 if flag=2A(1,1)=2;A(n,n)=2;B(1,1)=vl;B(n,1)=vr; end c=AB; for(i=1:n-1)d(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 mm,nn=size(xx); yy=zeros(mm,nn); for i=1:mm*nnfor ii=1:n-1if xx(i)=x(ii)break;elseif xx(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; end end
