計算方法大作業(yè)第二次

1、數(shù)值計算第二次大作業(yè)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作三次樣條函數(shù)插值，取第一類邊界條件Y0=0.01087 Y7=100根據(jù)題目要求，首先要構(gòu)造三次樣條函數(shù)，三次樣條函數(shù)的構(gòu)造過程如下：設(shè)有共n個插值節(jié)點，任意給定一組常數(shù)，要求構(gòu)造一個插值三次樣條函數(shù)，使得如下插值條件得以滿足：，i=0，1，n經(jīng)過插值點的三次樣條函數(shù)是一組三次多項式，即有：由節(jié)點處的連續(xù)性可知：由節(jié)點處的一階與二階光

2、滑性可知：又設(shè)，記，則。再根據(jù)邊界條件，從而可相繼解出用matlab編程，編寫三次樣條函數(shù)（見附錄），對第一題求解：>> 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逆時針旋轉(zhuǎn)座標(biāo)軸45o 保持（1）中結(jié)點和邊界條件的幾何關(guān)系不變，再次作三次樣條函數(shù)插值，畫出插值曲線的圖像。坐標(biāo)軸逆時針旋轉(zhuǎn)45°，相當(dāng)于節(jié)點順時針旋轉(zhuǎn)45°。設(shè)為旋轉(zhuǎn)前的坐標(biāo)，為旋轉(zhuǎn)后的坐標(biāo)，則可以得到如下關(guān)系：故旋轉(zhuǎn)后的節(jié)點坐標(biāo)為：>> 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端點處的一階導(dǎo)數(shù)為：>> 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則旋轉(zhuǎn)后的三次樣條的系數(shù)及圖像為：>> 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 旋轉(zhuǎn)后的三次樣條曲線3比較（1）、（2）的結(jié)果，能得到什么結(jié)論？將（1）中所得的三次樣條曲線整體順時針旋轉(zhuǎn)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('節(jié)點','旋轉(zhuǎn)前','旋轉(zhuǎn)后');圖3旋轉(zhuǎn)前后三次樣條曲線幾何比較比較圖中兩條曲線可知，曲線

9、不重合，故三次樣條插值不具備幾何不變性。附錄：三次樣條插值函數(shù)程序functionyy,b,c,d=spline3(x,y,xx,flag,vl,vr)%三次樣條插值函數(shù)%（x,y)為插值節(jié)點，xx為插值點%flag表端點邊界條件類型；%flag=0：自然樣條（端點二階導(dǎo)數(shù)為0）；%flag=1：第一類邊界條件（端點一階導(dǎo)數(shù)給定）；%flag=2：第二類邊界條件（端點二階導(dǎo)數(shù)給定）；%vl,vr表示左右端點處的在邊界條件值；%樣條函數(shù)為：Si(x)=yi+bi*(x-xi)+ci*(x-xi)2+di*(x-xi)3%b,c,d分別為各子區(qū)間上的系數(shù)值%yy表示插值點處的函數(shù)值if lengt

10、h(x)=length(y) error('輸入數(shù)據(jù)應(yīng)成對！');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、二階導(dǎo)數(shù)為0）if flag=0 A(1,1)=1; A(n,n)=1;end%端點一階導(dǎo)數(shù)條件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%端點二階導(dǎo)數(shù)條件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=