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

1、數(shù)值分析上機(jī)報(bào)告題目：插值法學(xué)號(hào)：201014924姓名：靳會(huì)有一、調(diào)用MATLAB內(nèi)帶函數(shù)插值1、MATLAB內(nèi)帶插值函數(shù)列舉如下：interp1interpftinterp2interp3interpnsplinemeshgridndgridgriddata一維數(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ù)（interp1）為例，程序如下：其調(diào)用格式為： yi=interp1(x, y, xi) yi=interp1(x, y, x

2、i, method) 舉例如下： x=0:10:100y=40 44 46 52 65 76 80 82 88 92 110;xi=0:1:100yi=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, Z, XI, YI)， ZI=imerp2(Z, ntimes)ZI=interp2(Z, XI, YI) ，ZI=interp2(X, Y, Z, XI, YI, method) interp3函數(shù)：VI=in

3、terp3(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, Y2, Y3, ) VI=interpn(, method)Spline函數(shù)：yi=spline(x,y,xi)pp=spline(x,y)meshgrid函數(shù)： X,Y=meshgrid(x,y) X,Y=meshgrid(x) X,Y

4、,Z=meshgrid(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, ZI=griddata(x, y, z, xi, yi) =griddata( method)二、自編函數(shù)插值1、拉格朗日插值法：建立M 文件：function f = Language(x,y,x0)syms t l;if(length(x) = length(y) n = length(x);else disp('x和y的維

5、數(shù)不相等！'); return; %檢錯(cuò)endh=sym(0);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(nargin = 3) f = subs (h,'t',x0); %計(jì)算插值點(diǎn)的函數(shù)值else f=collect(h); f = vpa(f,6); %將插值多項(xiàng)式的系數(shù)化成6位精度的小數(shù)end在MATLAB中輸入：x=18 31 66 68 7

6、0 72 70; y=23 33 52 51 43 40 46;f=Language(x,y)plot(x,y)結(jié)果為：f =Inf + (-t)*Inf - 54329.8*t2 + 1503.75*t3 - 22.2065*t4 + 0.16789*t5 - 0.000512106*t6圖形如下：MATLAB實(shí)現(xiàn)拉格朗日插值建立如下拉格朗日插值函數(shù)：function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0

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

8、輸入要計(jì)算的函數(shù)點(diǎn)并計(jì)算的值，利用牛頓插值公式，當(dāng)增加一個(gè)節(jié)點(diǎn)時(shí)，只需在后面多計(jì)算一項(xiàng)，而前面的計(jì)算仍有用；另一方面的各項(xiàng)系數(shù)恰好又是各階均差，而各階均差可用均差公式來(lái)計(jì)算。為 的 一階均差。為 的 k 階均差。均差表：零階均差一階均差二階均差三階均差X0f(X0)X1f(X1)fX0, X1X2f(X2)fX1, X2fX0,X1, X2X3f(X3)fX2, X3fX1, X2,X3fX0,X1, X2 X3MMMMM 牛頓插值法計(jì)算步驟：1 輸入值及（；要計(jì)算的函數(shù)點(diǎn)。2 對(duì)給定的由 計(jì)算的值。3輸出。 程序清單：functionc, d=newpoly(x, y) %牛頓插值的MATL

9、AB實(shí)現(xiàn) %這里 x為n個(gè)節(jié)點(diǎn)的橫坐標(biāo)所組成的向量，y為縱坐標(biāo)所組成的向量。 %c為所求的牛頓插值多項(xiàng)式的系數(shù)構(gòu)成的向量。n=length(x);%取x的個(gè)數(shù)。d=zeros(n, n);%構(gòu)造nXn的空數(shù)組。d(: , 1)=y'for j=2 : n for k=j : n d(k, j)=(d(k, j-1) - d(k-1, j-1) / (x(k)-x(k-j+1); endend c =d(n, n); for k=(n-1) : - 1 : 1 c =conv(c, poly(x(k);% conv求積，poly(x)將該多項(xiàng)式的系數(shù)賦給向量。 m=length(c);

10、c(m)=c(m)+d(k, k); end五、測(cè)試數(shù)據(jù)與結(jié)果：測(cè)試數(shù)據(jù)：（第三章習(xí)題第三題第2題）f(x)=lnx的數(shù)值如表所示， 構(gòu)造牛頓插值多項(xiàng)式并求ln0.53的值。X0.40.50.60.70.8lnx-0.916291-0.693147-0.510826-0.357765-0.223144解： 由表可知x0=0.4, x1=0.5, x2=0.6, x3=0.7, x4=0.7，函數(shù)值：Y0=-0.916291, y1=-0.693147, y2=-0.510826, y3=-0.357765, y4=-0.223144建立一個(gè)主程序np.mclcclearnewpoly(0.4,0.5,0.6,0.7,0.8, -0.916291, -0.693147, -0