MATLAB實現(xiàn)拉格朗日插值

1、數(shù)值分析上機(jī)報告題目：插值法學(xué)號：201014924姓名：靳會有一、調(diào)用MATLAB內(nèi)帶函數(shù)插值1、MATLAB內(nèi)帶插值函數(shù)列舉如下:interpl interpft interp2 interp3 interpn spline meshgrid ndgrid griddata一維數(shù)據(jù)內(nèi)插（查表法） 使用FFT方法的一維數(shù)據(jù)內(nèi)插 二維數(shù)據(jù)內(nèi)插（查表法） 三維數(shù)據(jù)內(nèi)插（查表法） 多維數(shù)據(jù)內(nèi)插（查表法） 三次樣條內(nèi)插為三維繪圖產(chǎn)生X和Y陣 為多維函數(shù)和內(nèi)插產(chǎn)生陣列 數(shù)據(jù)網(wǎng)格2、取其中的一維數(shù)據(jù)內(nèi)插函數(shù)（interpl ）為例，程序如下:其調(diào)用格式為：yi=interp1(x, y, xi)yi=i

2、nterp1(x, y, xi, method)舉例如下：x=0:10:100y=40 44 46 52 65 76 80 82 88 92 110;xi=0:1:100 yi=interp1(x,y,xi,'spline')3、其他內(nèi)帶函數(shù)調(diào)用格式為：Interpft 函數(shù)：y=interpft(x,n)y=interpft(x,n,dim)interp2 函數(shù)：ZI=interp2(X, Y,乙 XI, YI) ， ZI=imerp2(Z, ntimes)Zl=interp2(Z,XI, YI) ，ZI=interp2(X, Y, Z, XI, YI,method) int

3、erp3 函數(shù)：VI=interp3(X,Y,Z,V,XI,YI,ZI) VI=interp3(V, ntimes)VI=interp3(V,XI,YI,ZI)VI=interp3(,method)Interpn 函數(shù)：VI=interpn(X1, X2, X3,V, Y1, Y2, Y3,)VI=interpn(V, ntimes)VI=interpn(V, Yl, 丫2, 丫3,) method)Spline 函數(shù)：yi=spline(x,y,xi)pp=spli ne(x,y)meshgrid 函數(shù)：X,Y=meshgrid(x,y)X,Y=meshgrid(x)X,Y,Z=meshgr

4、id(x,y,z)Ndgrid 函數(shù)：X1, X2, X3,=ndgrid(x1, x2, x3,X1, X2, X3,=ndgrid(x)Griddata 函數(shù)：ZI=griddata(x, y, z, XI, YI)XI, YI, Zl=griddata(x, y, z, xi, yi)=griddata( method)二、自編函數(shù)插值1、拉格朗日插值法：建立M文件：fun cti on f = Lan guage(x,y,x0)syms t l;if(len gth(x) = len gth(y)VI=interpn()n = len gth(x);elsedisp('x和y的

5、維數(shù)不相等！');return;% 檢錯endh=sym(O);for (i=1: n)l=sym(y(i);for(j=1:i-1)l=l*(t-x(j)/(x(i)-x(j);end;for(j=i+1: n)l=l*(t-x(j)/(x(i)-x(j);end;h=h+l;endsimplify(h);if(n argin = 3)f = subs (h,'t',xO); %計算插值點的函數(shù)值elsef=collect(h);f = vpa(f,6); % 將插值多項式的系數(shù)化成6位精度的小數(shù)end在MATLAB中輸入：x=18 31 66 68 70 72 70

6、;y=23 33 52 51 43 40 46;f=La nguage(x,y)Plot(x,y)結(jié)果為：f =I nf + (-t)*I nf - 54329.8*2 + 1503.75*tA3 - 22.2065P4 + 0.16789*tA5-0.000512106*tA6圖形如下：MATLAB實現(xiàn)拉格朗日插值建立如下拉格朗日插值函數(shù)：fun cti on y=lagra nge(x0,y0,x);n=len gth(x0);m=le ngth(x);for i=1:mz=x(i);s=0.0;for k=1: np=1.0;for j=1: nif j二k p=p*(z-xO(j)/(

7、xO(k)-xO(j);endends=p*yO(k)+s;endy(i)=s;end畫圖程序如下：x=-5:1:5; y=1./(1+x.A2);xO二卜5:0.001:5; y0=lagra nge(x,y,x0); y1=1./(1+xO.A2);plot(x0,y0,'r')hold onplot(xO,y1,'g')注：畫出的圖形為n =10的圖形得到圖形如下：牛頓K次插值多項式一、實驗?zāi)康模?、掌握牛頓插值法的基本思路和步驟。2、培養(yǎng)編程與上機(jī)調(diào)試能力。二、牛頓插值法基本思路與計算步驟:給定插值點序列(Xi, f(Xi) ,i 0,1,n,。構(gòu)造牛頓

8、插值多項式 Nn(u)。輸入要計算的函數(shù)點X,并計算Nn(x)的值，利用牛頓插值公式，當(dāng)增加一個節(jié)點時，只需在后面多計算 一項，而前面的計算仍有用；另一方面Nn(X)的各項系數(shù)恰好又是各階均差，而各階均差可用均差公式來計算。/(X)一階均差。/坯X.L斗二丿氐丄 兀牛耳一/1%斗丄兀討/(X)的k階均差。均差表:Xk零階均差一階均差二階均差三階均差Xof(Xo)X1f(X1)fXo, X1X2f(X2)fX 1, X2fXo,X1, X2X3f(X3)fX2, X3fX 1, X2 ,X3fX o,X1, X2X3MMMMM牛頓插值法計算步驟:1輸入n值及(Xi, f (Xi) ,i0,1,m

9、 ；要計算的函數(shù)點x。2.對給定的x，由Nn(X) f (Xo) (X Xo) f Xo,X1(X Xo)(X xjf Xo,X1, X2L(X Xo)(X XjL (x Xn 1)fXo,X1L , Xn計算Nn(X)的值。3 輸出 Nn(x)。程序清單: functionc, d=n ewpoly(x, y)%牛頓插值的MATLAB實現(xiàn)%這里x為n個節(jié)點的橫坐標(biāo)所組成的向量，y為縱坐標(biāo)所組成的向量。%c為所求的牛頓插值多項式的系數(shù)構(gòu)成的向量。n=length(x);%取 x 的個數(shù)。d=zeros(n, n);% 構(gòu)造nXn的空數(shù)組。d(: , 1)=y'for j=2 : nfo

10、r k=j : nd(k, j)=(d(k, j-1) - d(k-1, j-1) / (x(k)-x(k-j+1);endendc =d( n, n);for k=(n-1): - 1 : 1c =conv(c, poly(x(k);% conv 求積，poly(x)將該多項式的系數(shù)賦給向量。 m=le ngth(c);c(m)=c(m)+d(k, k);end五、測試數(shù)據(jù)與結(jié)果：測試數(shù)據(jù)：(第三章習(xí)題第三題第 2題)f(x)=lnx的數(shù)值如表所示，構(gòu)造牛頓插值多項式并求ln0.53的值。X0.40.50.60.70.8In-0.91-0.69-0.5108-0.3577-0.22314x6291314726654解：由表可知 X0=0.4, x 1=0.5, x 2=0.6, x 3=0.7, x 4=0.7，函數(shù)值：Y0=-0.916291,y1=-0.693147,y2=-0.510826,y3=-0.357765,y4=-0.223144建立一個主程序 n p.m clcclearn ewpoly(0.4,0.5,0.6,0.7,0.