MATLAB实现拉格朗日插值

MATLAB 实现拉格朗日 插值 Document serial number【NL89WT-NY98YT-NC8CB-NNUUT-NUT108】 数值分析上机报告 题目：插值法 学号：4 姓名：靳会有 一、调用 MATLAB 内带函数插值 1、MATLAB 内带插值函数列举如下： interp1 interpft interp2 interp3 interpn spline meshgrid ndgrid griddata 一维数据内插(查表法) 使用 FFT 方法的一维数据内插 二维数据内插(查表法) 三维数据内插(查表法) 多维数据内插(查表法) 三次样条内插 为三维绘图产生 X 和 Y 阵 为多维函数和内插产生阵列 数据网格 2、取其中的一维数据内插函数（interp1）为例，程序如下： 其调用格式为： yi=interp1(x, y, xi) yi=interp1(x, y, xi, ) 举例如下： x=0:10:100 y=[40 44 46 52 65 76 80 82 88 92 110]; xi=0:1:100 yi=interp1(x,y,xi, spline ) 3、其他内带函数调用格式为： Interpft 函数： y=interpft(x,n) y=interpft(x,n,dim) interp2 函数： ZI=interp2(X, Y, Z, XI, YI)， ZI=imerp2(Z, ntimes) ZI=interp2(Z, XI, YI) ，ZI=interp2(X, Y, Z, XI, YI, ) interp3 函数： VI=interp3(X,Y,Z,V,XI,YI,ZI)VI=interp3(V, ntimes) VI=interp3(V,XI,YI,ZI)VI=interp3(…, ) Interpn 函数： VI=interpn(X1, X2, X3, …, V, Y1, Y2, Y3, …) VI=interpn(V, ntimes) VI=interpn(V, Yl, Y2, Y3, …)VI=interpn(…, ) Spline 函数： yi=spline(x,y,xi) pp=spline(x,y) meshgrid 函数： [X,Y]=meshgrid(x,y) [X,Y]=meshgrid(x) [X,Y,Z]=meshgrid(x,y,z) Ndgrid 函数： [X1, X2, X3, …]=ndgrid(x1, x2, x3, …) [X1, X2, X3, …]=ndgrid(x) Griddata 函数： ZI=griddata(x, y, z, XI, YI) [XI, YI, ZI]=griddata(x, y, z, xi, yi) […]=griddata(… ) 二、自编函数插值 1、拉格朗日插值法： 建立 M 文件： function f = Language(x,y,x0) syms t l; if(length(x) == length(y)) n = length(x); else disp( x 和 y 的维数不相等！ ); return; %检错 end h=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; end simplify(h); if(nargin == 3) f = subs (h, t ,x0); %计算插值点的函数值 else f=collect(h); f = vpa(f,6); %将插值多项式的系数化成 6 位精度的小数 end 在 MATLAB 中输入： x=[18 31 66 68 70 72 70;] y=[23 33 52 51 43 40 46]; f=Language(x,y) plot(x,y) 结果为： f =Inf + (-t)*Inf - *t^2 + *t^3 - *t^4 + *t^5 - *t^6 图形如下：MATLAB 实现拉格朗日插值 建立如下拉格朗日插值函数： function y=lagrange(x0,y0,x); n=length(x0); m=length(x); for i=1:m z=x(i); s=; for k=1:n p=; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end 画图程序如下： x=[-5:1:5]; y=1./(1+x.^2); x0=[-5::5]; y0=lagrange(x,y,x0); y1=1./(1+x0.^2); plot(x0,y0, r ) hold on plot(x0,y1, g ) 注：画出的图形为 n =10 的图形 得到图形如下： n=10 的图 像 牛顿 K 次插值多项式 一、实验目的： 1、掌握牛顿插值法的基本思路和步骤。 2、培养编程与上机调试能力。 二、牛顿插值法基本思路与计算步骤： 给定插值点序列（x i , f (x i )),i0,1,,n,。构造牛顿插值多项式N n (u)。输 入要计算的函数点x,并计算N n (x)的值，利用牛顿插值公式，当增加一个节点 时，只需在后面多计算一项，而前面的计算仍有用；另一方面N n (x)的各项系数 恰好又是各阶均差，而各阶均差可用均差公式来计算。 为的 一阶均差。 为 均差表： 差 X0 X1 X2 X3 M f(X0) f(X1) f(X2) f(X3) M 零 阶 均一阶均差二阶均差 的 k 阶均差。 三阶均差 f[X0, X1] f[X1, X2] f[X2, X3] M f[X0,X1, X2] f[X1, X2,X3] M f[X0,X1, X2 X3] M 牛顿插值法计算步骤： x , f (x i )),i 0,1,,n, 1．输入n值及（ i；要计算的函数点 x 。 2．对给定的 x, 由 N (x) 计算 n 的值。 N (x) 3．输出 n 。 程序清单： function[c, d]=newpoly(x, y) %牛顿插值的 MATLAB 实现 %这里 x 为 n 个节点的横坐标所组成的向量，y 为纵坐标所组成的向量。 %c 为所求的牛顿插值多项式的系数构成的向量。 n=length(x);%取 x 的个数。 d=zeros(n, n);%构造 nXn 的空数组。 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)); end end c =d(n, n); for k=(n-1) : - 1 : 1 c =conv(c, poly(x(k)));% conv 求积，poly(x)将该多项式的系数 赋