数值计算第二次大作业
1．给定插值条件如下：
i 0 1 2 3 4 5 6 7 Xi 8.125 8.4 9.0 9.485 9.6 9.959 10.166 10.2 Yi 0.0774 0.099 0.280 0.60 0.708 1.200 1.800 2.177 作三次样条函数插值，取第一类边界条件 Y0’=0.01087 Y7’=100
根据题目要求，首先要构造三次样条函数，三次样条函数的构造过程如下：
设有x0＜x1＜…＜xn共n个插值节点，任意给定一组常数y0，y1，…，yn，要求构造一个插值三次样条函数S(x)，使得如下插值条件得以满足：
S(xi)?yi，i=0，1，…，n
经过插值点的三次样条函数是一组三次多项式，即有：
23?S1(x)?a1?b1(x?x1)?c1(x?x1)?d1(x?x1),x?[x1,x2],?23S2(x)?a2?b2(x?x2)?c2(x?x2)?d2(x?x2),x?[x2,x3],? ????S(x)?a?b(x?x)?c(x?x)2?d(x?x)3,x?[x,x]n?1n?1n?1n?1n?1n?1n?1n?1n?n?1由节点处的连续性可知：
Si(xi)?yi,Si(xi?1)?yi?1,i?1,2,?n?1.
ai?yi,i?1,2,?n?1,??23y2?y1?b1(x2?x1)?c1(x2?x1)?d1(x2?x1)? ????y?y?b(x?x)?c(x?x)2?d(x?x)3n?1n?1nn?1n?1nn?1n?1nn?1?n由节点处的一阶与二阶光滑性可知：
Si?1(xi)?Si(xi),Si?1(xi)?Si(xi),i?1,2,?,n
''''''又设cn?Sn?1(xn)/2，记?i?xi?1?xi,?i?yi?1?yi,i?1,2,?,n?1，则
di?ci?1?ci3?i,i?1,2,?,n?1。再根据边界条件，从而可相继解出bi,ci
''用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('\\t\\tb1\\t\\tc1\\t\\td1\\n'); b1 c1 d1 >> 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．）中结点和边界条件的几何关系不变，再次作三次样条函数插值，画出插值曲线的图像。
坐标轴逆时针旋转45°，相当于节点顺时针旋转45°。设(x,y)T为旋转前的坐标，
(x',y')T为旋转后的坐标，则可以得到如下关系：
?sin???x????cos???y??x'??cos??????y'??sin?
故旋转后的节点坐标为：
>> 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('\\t\\t\\tx2\\t\\t\\ty2\\n'); >> 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)); >> fprintf('\\t\\t\\tv1\\t\\t\\tvn\\n'); v1 vn >> 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('\\t\\t\\tb2\\t\\t\\tc2\\t\\t\\td2\\n'); b2 c2 d2 >> disp([b2 c2(1:end-1,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．）中的三次样条曲线画在同一幅图中比较，得
>> 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旋转前后三次样条曲线几何比较
比较图中两条曲线可知，曲线不重合，故三次样条插值不具备几何不变性。
附录：
三次样条插值函数程序
function[yy,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 length(x)~=length(y)
error('输入数据应成对！'); end
n=length(x); a=zeros(n-1,1); b=a;d=a;dx=a;dy=a;