MATLAB实现拉格朗日插值

数值分析上机报告

标准文案

题目：插值法学号：姓名：靳会有

201014924

实用文档

一、调用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, method)

举例如下：

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, method) interp3函数：

VI=interp3(X,Y,Z,V,XI,YI,ZI) VI=interp3(V, ntimes)

VI=interp3(V,XI,YI,ZI) VI=interp3(…, method)

Interpn函数：

VI=interpn(X1, X2, X3, …, V, Y1, Y2, Y3, …) VI=interpn(V, ntimes)

VI=interpn(V, Yl, Y2, Y3, …) VI=interpn(…, method)

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(… method) 二、自编函数插值 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 - 54329.8*t^2 + 1503.75*t^3 - 22.2065*t^4 + 0.16789*t^5 - 0.000512106*t^6

图形如下：

MATLAB实现拉格朗日插值

建立如下拉格朗日插值函数： 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(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:0.001: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的图像

牛顿一、实验目的：

1、掌握牛顿插值法的基本思路和步骤。 2、 培养编程与上机调试能力。

二、牛顿插值法基本思路与计算步骤：

给定插值点序列（xi,f(xi)),i0,1,,n,。构造牛顿插值多项式Nn(u)。输入要计算的函数点x,并计算Nn(x)的值，利用牛顿插值公式，当增加一个节点时，只需在后面多计算一项，而前面的计算仍有用；另一方面Nn(x)的各项系数恰好又是各阶均差，而各阶均差可用均差公式来计算。

K次插值多项式

为 的 一阶均差。

为

均差表：

的 k 阶均差。

x k标准文案

零阶均一阶均差 二阶均差 三阶均差 实用文档

差 X0 X1 X2 X3 M f(X0) f(X1) f(X2) f(X3) M 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 牛顿插值法计算步骤： 1． 输入n值及（i2． 对给定的x,由

x,f(xi)),i0,1,,n,；要计算的函数点x。

Nn(x)f(x0)(xx0)fx0,x1(xx0)(xx1)fx0,x1,x2L(xx0)

(xx1)L(xxn1)fx0,x1L,xn计算

Nn(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)将该多项式的系数赋给向量。 m=length(c);

c(m)=c(m)+d(k, k); end

五、测试数据与结果：

测试数据：（第三章习题第三题第2题）

f(x)=lnx的数值如表所示， 构造牛顿插值多项式并求ln0.53的值。 X 0.4 0.5 0.6 0.7 0.8 ln-0.91-0.69-0.5108-0.3577x 6291 3147 26 65 解： 由表可知x0=0.4, x1=0.5, x2=0.6, x3=0.7, x4=0.7，函数值： -0.223144 N(x)标准文案

实用文档

Y0=-0.916291, y1=-0.693147, y2=-0.510826, y3=-0.357765, y4=-0.223144 建立一个主程序np.m clc clear

newpoly([0.4,0.5,0.6,0.7,0.8],[ -0.916291, -0.693147, -0.510826, -0.357765, -0.223144])

计算结果如下： ans =

-0.3096 2.6083 -5.4861 5.6921 -2.4744 由此看出所求的牛顿多项式为：

432

P(x)= -0.3096x+2.6083x-5.4861x+5.6921x-2.4744 P(0.53)= -0.6347。

标准文案