數(shù)值分析作業(yè)-三次樣條插值

1、數(shù)值計算方法作業(yè)實驗名稱實驗4.3三次樣條插值函數(shù)（P126）4.5三次樣條插值函數(shù)的收斂性（P127）實驗時間姓名班級學號成績實驗4.3 三次樣條差值函數(shù)實驗目的：掌握三次樣條插值函數(shù)的三彎矩方法。實驗函數(shù):x0.00.10.20.30.4F(x)0.50000.53980.57930.61790.7554求f(0.13)和f(0.36)的近似值實驗內(nèi)容：（1） 編程實現(xiàn)求三次樣條插值函數(shù)的算法，分別考慮不同的邊界條件；（2） 計算各插值節(jié)點的彎矩值；（3） 在同一坐標系中繪制函數(shù)f(x)，插值多項式，三次樣條插值多項式的曲線比較插值結(jié)果。實驗4.5 三次樣條差值函數(shù)的收斂性實驗目的：多項式

2、插值不一定是收斂的，即插值的節(jié)點多，效果不一定好。對三次樣條插值函數(shù)如何呢？理論上證明三次樣條插值函數(shù)的收斂性是比較困難的，通過本實驗可以證明這一理論結(jié)果。實驗內(nèi)容：按照一定的規(guī)則分別選擇等距或非等距的插值節(jié)點，并不斷增加插值節(jié)點的個數(shù)。實驗要求：（1） 隨著節(jié)點個數(shù)的增加，比較被逼近函數(shù)和三樣條插值函數(shù)的誤差變化情況，分析所得結(jié)果并與拉格朗日插值多項式比較；（2） 三次樣條插值函數(shù)的思想最早產(chǎn)生于工業(yè)部門。作為工業(yè)應用的例子，考慮如下例子：某汽車制造商根據(jù)三次樣條插值函數(shù)設計車門曲線，其中一段數(shù)據(jù)如下：0123456789100.00.791.532.192.713.033.272.893.

3、063.193.290.80.2算法描述：拉格朗日插值： 其中是拉格朗日基函數(shù)，其表達式為：牛頓插值：其中三樣條插值：所謂三次樣條插值多項式Sn(x)是一種分段函數(shù)，它在節(jié)點Xi(a<X0<X1<Xn<b)分成的每個小區(qū)間xi-1,xi上是三次多項式，其在此區(qū)間上的表達式如下：式中Mi=.因此，只要確定了Mi的值，就確定了整個表達式，Mi的計算方法如下：令則Mi滿足如下n-1個方程：常用的邊界條件有如下幾類：（1） 給定區(qū)間兩端點的斜率m0,mn,即（2） 給定區(qū)間兩端點的二階導數(shù)M0，Mn,即（3） 假設y=f(x)是以b-a為周期的周期函數(shù)，則要求三次樣條插值函數(shù)S

4、（x）也為周期函數(shù)，對S（x）加上周期條件對于第一類邊界條件有對于第二類邊界條件有其中那么解就可以為 對于第三類邊界條件，由此推得，其中，那么解就可以為：程序代碼：1拉格朗日插值函數(shù)Lang.mfunction f=lang(X,Y,xi)%X為已知數(shù)據(jù)的橫坐標%Y為已知數(shù)據(jù)的縱坐標%xi插值點處的橫坐標%f求得的拉格朗日插值多項式的值 n=length(X); f=0; for i=1:n l=1; for j=1:i-1 l=l.*(xi-X(j)/(X(i)-X(j); end; for j=i+1:n l=l.*(xi-X(j)/(X(i)-X(j); end;%拉格朗日基函數(shù) f=f

5、+l*Y(i); end fprintf('%dn',f)return2 牛頓插值函數(shù)newton.mfunction f=newton(X,Y,xi)%X為已知數(shù)據(jù)的橫坐標%Y為已知數(shù)據(jù)的縱坐標%xi插值點處的橫坐標%f求得的拉格朗日插值多項式的值n=length(X);newt=X',Y'%計算差商表for j=2:n for i=n:-1:1 if i>=j Y(i)=(Y(i)-Y(i-1)/(X(i)-X(i-j+1); else Y(i)=0; end end newt=newt,Y'end %計算牛頓插值 f=newt(1,2); f

6、or i=2:n z=1; for k=1:i-1 z=(xi-X(k)*z; end f=f+newt(i-1,i)*z; end fprintf('%dn',f) return3三次樣條插值第一類邊界條件Threch.mfunction S=Threch1(X,Y,dy0,dyn,xi) % X為已知數(shù)據(jù)的橫坐標%Y為已知數(shù)據(jù)的縱坐標%xi插值點處的橫坐標%S求得的三次樣條插值函數(shù)的值 %dy0左端點處的一階導數(shù)% dyn右端點處的一階導數(shù)n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros

7、(1,n-2);for i=1:n%求函數(shù)的一階差商 h(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i)/h(i);end for i=2:n%求函數(shù)的二階差商 f2(i)=(f1(i)-f1(i-1)/(X(i+1)-X(i-1);d(i)=6*f2(i);end d(1)=6*(f1(1)-dy0)/h(1); d(n+1)=6*(dyn-f1(n-1)/h(n-1);%¸賦初值A=zeros(n+1,n+1); B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1 B(i)=h(i)/(h(i)+h(i+1);C(i)=1-B

8、(i);end A(1,2)=1;A(n+1,n)=1;for i=1:n+1 A(i,i)=2;end for i=2:n A(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=Ad;syms x;for i=1:n Sx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(x-X(i). +M(i)/2*(x-X(i)2+(M(i+1)-M(i)/(6*h(i)*(x-X(i)3);digits(4); Sx(i)=vpa(Sx(i);%三樣條插值函數(shù)表達式end for i=1:n disp('S(x)=');

9、fprintf('%s (%d,%d)n',char(Sx(i),X(i),X(i+1);end for i=1:n if xi>=X(i)&&xi<=X(i+1) S=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(xi-X(i)+M(i)/2*(xi-X(i)2+(M(i+1)-M(i)/(6*h(i)*(xi-X(i)3; end enddisp('xi S');fprintf('%d,%dn',xi,S);return 4 三次樣條插值第二類邊界條件Threch2.mfunction S

10、x=Threch2(X,Y,d2y0,d2yn,xi) X為已知數(shù)據(jù)的橫坐標%Y為已知數(shù)據(jù)的縱坐標%xi插值點處的橫坐標%S求得的三次樣條插值函數(shù)的值 %d2y0左端點處的二階導數(shù)% d2yn右端點處的二階導數(shù)n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求一階差商h(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i)/h(i);end for i=2:n%求二階差商f2(i)=(f1(i)-f1(i-1)/(X(i+1)-X(i-1);d(i)=6*

11、f2(i);end d(1)=2*d2y0; d(n+1)=2*d2yn;%賦初值 A=zeros(n+1,n+1); B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1 B(i)=h(i)/(h(i)+h(i+1);C(i)=1-B(i);end A(1,2)=0;A(n+1,n)=0;for i=1:n+1 A(i,i)=2;end for i=2:n A(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=Ad;syms x;for i=1:n Sx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i

12、)*(x-X(i). +M(i)/2*(x-X(i)2+(M(i+1)-M(i)/(6*h(i)*(x-X(i)3);digits(4); Sx(i)=vpa(Sx(i);end for i=1:n disp('S(x)=');fprintf('%s (%d,%d)n',char(Sx(i),X(i),X(i+1);end for i=1:n if xi>=X(i)&&xi<=X(i+1) S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(xi-X(i)+M(i)/2*(xi-X(i)2+(M(i+1)

13、-M(i)/(6*h(i)*(xi-X(i)3; end enddisp('xi S');fprintf('%d,%dn',xi,S);return5插值節(jié)點處的插值結(jié)果main3.mclearclcX=0.0,0.1,0.2,0.3,0.4;Y=0.5000,0.5398,0.5793,0.6179,0.7554;xi=0.13;%xi=0.36;disp('xi=0.13');%disp('xi=0.36');disp('拉格朗日插值結(jié)果');lang(X,Y,xi);disp('牛頓插值結(jié)果'

14、);newton(X,Y,xi);disp('三次樣條第一類邊界條件插值結(jié)果');Threch1(X,Y,0.40,0.36,xi);%0.4,0.36分別為兩端點處的一階導數(shù)disp('三次樣條第二類邊界條件插值結(jié)果');Threch2(X,Y,0,-0.136,xi);%0,-0.136分別為兩端點處的二階導數(shù)6將多種插值函數(shù)即原函數(shù)圖像畫在同一張圖上main2.mclearclcX=0.0,0.1,0.2,0.3,0.4;Y=0.5000,0.5398,0.5793,0.6179,0.7554;a=linspace(0,0.4,21);NUM=21;L=z

15、eros(1,NUM);N=zeros(1,NUM);S=zeros(1,NUM);B=zeros(1,NUM);for i=1:NUM xi=a(i); L(i)=lang(X,Y,xi);% 拉格朗日插值 N(i)=newton(X,Y,xi);% 牛頓插值 B(i)=normcdf(xi,0,1);%原函數(shù) S(i)=Threch1(X,Y,0.4,0.36,xi);%三次樣條函數(shù)第一類邊界條件endplot(a,B,'-r');hold on;plot(a,L,'b');hold on;plot(a,N,'r');hold on;plot

16、(a,S,'r+');hold on;legend('原函數(shù)','拉格朗日插值','牛頓插值','三次樣條插值',2);hold off7增加插值節(jié)點觀察誤差變化main4.mclear;clc;N=5;%4.5第一問Ini=zeros(1,1001);a=linspace(-1,1,1001);Ini=1./(1+25*a.2);for i=1:3 %節(jié)點數(shù)量變化次數(shù) N=2*N; t=linspace(-1,1,N+1);%插值節(jié)點 ft=1./(1+25*t.2);%插值節(jié)點函數(shù)值 val=linspace(

17、-1,1,101); for j=1:101 L(j)=lang(t,ft,val(j); S(j)=Threch1(t,ft,0.074,-0.074,val(j);%三樣條第一類邊界條件插值 endplot(a,Ini,'k')%原函數(shù)圖象hold onplot(val,L,'r')%拉格朗日插值函數(shù)圖像hold onplot(val,S,'b')%三次樣條插值函數(shù)圖像str=sprintf('插值節(jié)點為%d時的插值效果',N);title(str); legend('原函數(shù)','拉格朗日插值'

18、,'三次樣條插值');%顯示圖例hold off figureend8車門曲線main5.mclearclcX=0,1,2,3,4,5,6,7,8,9,10;Y=0.0,0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29;dy0=0.8;dyn=0.2;n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:nh(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i)/h(i);end for i=2:nf2(i)=(f1(i)-f1(i-1)/(X(i+1)-X(i-1);d(i)=6*f2(i);end d(1)=6*(f1(1)-dy0)/h(1); d(n+1)=6*(dyn-f1