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

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

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

3、063.193.290.80.2算法描述：拉格朗日插值： 其中是拉格朗日基函數(shù)，其表達(dá)式為：牛頓插值：其中三樣條插值：所謂三次樣條插值多項(xiàng)式Sn(x)是一種分段函數(shù)，它在節(jié)點(diǎn)Xi(aX0X1Xn=j Y(i)=(Y(i)-Y(i-1)/(X(i)-X(i-j+1); else Y(i)=0; end end newt=newt,Y;end %計(jì)算牛頓插值 f=newt(1,2); for 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三次樣條插值第一類邊界條件Threc

4、h.mfunction S=Threch1(X,Y,dy0,dyn,xi) % X為已知數(shù)據(jù)的橫坐標(biāo)%Y為已知數(shù)據(jù)的縱坐標(biāo)%xi插值點(diǎn)處的橫坐標(biāo)%S求得的三次樣條插值函數(shù)的值 %dy0左端點(diǎn)處的一階導(dǎo)數(shù)% dyn右端點(diǎn)處的一階導(dǎo)數(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%求函數(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

5、+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(chǔ)=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)=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(

6、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ù)表達(dá)式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)&xi=X(i+1) S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(xi-X(i

7、)+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);return5插值節(jié)點(diǎn)處的插值結(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é)果);newton(X,Y,xi);disp(三次樣條第一類邊界條件插

8、值結(jié)果);Threch1(X,Y,0.40,0.36,xi);%0.4,0.36分別為兩端點(diǎn)處的一階導(dǎo)數(shù)disp(三次樣條第二類邊界條件插值結(jié)果);Threch2(X,Y,0,-0.136,xi);%0,-0.136分別為兩端點(diǎn)處的二階導(dǎo)數(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=zeros(1,NUM);N=zeros(1,NUM);S=zeros(1,NUM);B=zeros(1,NU

9、M);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(a,S,r+);hold on;legend(原函數(shù),拉格朗日插值,牛頓插值,三次樣條插值,2);hold off7增加插值節(jié)點(diǎn)觀察誤差變化main4.mclear;

10、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é)點(diǎn)數(shù)量變化次數(shù) N=2*N; t=linspace(-1,1,N+1);%插值節(jié)點(diǎn) ft=1./(1+25*t.2);%插值節(jié)點(diǎn)函數(shù)值 val=linspace(-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(va

11、l,L,r)%拉格朗日插值函數(shù)圖像hold onplot(val,S,b)%三次樣條插值函數(shù)圖像str=sprintf(插值節(jié)點(diǎn)為%d時(shí)的插值效果,N);title(str); legend(原函數(shù),拉格朗日插值,三次樣條插值);%顯示圖例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);f

12、1=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(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(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;x=zeros(1,n);S=zeros(1,n);for i=1:n x(i)=X(i)+0.5; S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(x(i)-X(i)+M(i)/2*(x(i)-X(i)2+(M(i+1)-M(i)/(6*h(i)*(x(i)-X