计算方法大作业第二次

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

3、,y1,xx1,1,u1,un);>> fprintf('ttb1ttc1ttd1n');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,

4、9;r-');>> grid on画出插值曲线的图像。图1 三次样条曲线2逆时针旋转座标轴45o 保持（1）中结点和边界条件的几何关系不变，再次作三次样条函数插值，画出插值曲线的图像。坐标轴逆时针旋转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

5、(i);end>>fprintf('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.4612 -5.9157 8.7519 -5.6731端点处的一阶导数为：>> v1=(u1+tan(theta)/(1-u1*tan(theta);>> vn=(un+tan(theta)/(1-un*tan(theta);>> f

6、printf('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

7、138 -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）中的三次样条曲线画在同一幅图中

8、比较，得>> 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旋转前后三次样条曲线几何比较比较图中两条曲线可知，曲线

9、不重合，故三次样条插值不具备几何不变性。附录：三次样条插值函数程序functionyy,b,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 lengt

10、h(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);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%自然样条端点条件（端点

11、二阶导数为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)=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=