计算方法大作业(第二次)

《计算方法大作业(第二次)》由会员分享，可在线阅读，更多相关《计算方法大作业(第二次)（8页珍藏版）》请在装配图网上搜索。

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个插值节点，任意给定一组常数，要求构造一个插值三次样条函数，使得如下插值条件得以满足：，i=0，1，n经过插值点的三次样条函数是一组三次多项式，即有：

2、由节点处的连续性可知：由节点处的一阶与二阶光滑性可知：又设，记，则。再根据边界条件，从而可相继解出用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=0.01087;un=100; xx1=x1(1):0.001:x1(end); yy1 b1 c1 d1=spline3(x1,y1,xx1,1,u1,un); fprintf(ttb1ttc1ttd1

3、n);b1c1d1 disp(b1 c1(1:end-1,1) d1); 0.01087 0.14489 0.368 0.17405 0.4485 -0.393 0.2878 -0.25891 2.1153 1.5294 2.8188 -69.141 -0.56548 -21.035 73.614 12.794 58.247 -512.32 -28.949 -259.9 42279 plot(x1,y1,bo,xx1,yy1,r-); grid on画出插值曲线的图像。图1 三次样条曲线2逆时针旋转座标轴45o 保持（1）中结点和边界条件的几何关系不变，再次作三次样条函数插值，画出插值曲线的图

4、像。坐标轴逆时针旋转45，相当于节点顺时针旋转45。设为旋转前的坐标，为旋转后的坐标，则可以得到如下关系：故旋转后的节点坐标为： 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);endfprintf(tttx2ttty2n); disp(x2 y2); 5.8 -5.6905 6.0097 -5.8697 6.562 -6.166 7.1312 -6.2826 7.2889 -6.2876 7.8906 -6.1935 8

5、.4612 -5.9157 8.7519 -5.6731端点处的一阶导数为： v1=(u1+tan(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

6、,1) d2); -0.97849 0.67221 -0.38277 -0.74704 0.43138 -0.090754 -0.35362 0.28102 -0.034909 -0.067629 0.22141 0.053338 0.0061747 0.24664 0.0046897 0.3081 0.2551 0.10233 0.6992 0.43028 0.12195 plot(x2,y2,b+,xx2,yy2,m-.); grid on;图2 旋转后的三次样条曲线3比较（1）、（2）的结果，能得到什么结论？将（1）中所得的三次样条曲线整体顺时针旋转45，并与二题（2）中的三次样条曲线画

7、在同一幅图中比较，得 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,c,d=spline3(x,y,xx,flag,vl,vr)%三次样条插值函数%

8、（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(输入数据应成对！);endn=length(x);a=zeros(n-1,1);b=a;d=a;dx=a;dy=a;A=zeros

9、(n);B=zeros(n,1);for i=1:n-1 a(i)=y(i); dx(i)=x(i+1)-x(i); dy(i)=y(i+1)-y(i);endfor i=2:n-1 A(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=0 A(1,1)=1; A(n,n)=1;end%端点一阶导数条件if flag=1 A(1,1)=2*dx(1);A(1,2)=dx(1); A(n,n-1)=

10、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=2 A(1,1)=2;A(n,n)=2; B(1,1)=vl;B(n,1)=vr;endc=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;endmm,nn=size(xx);yy=zeros(mm,nn);for i=1:mm*nn for ii=1:n-1 if xx(i)=x(ii)&xx(i)x(ii+1) j=ii; break; elseif xx(i)=x(n) j=n-1; end end yy(i)=a(j)+b(j)*(xx(i)-x(j)+c(j)*(xx(i)-x(j)2+d(j)*(xx(i)-x(j)3;endend第 8 页 共 8 页免责声明：图文来源于网络搜集，版权归原作者所以若侵犯了您的合法权益，请作者与本上传人联系，我们将及时更正删除。