第二章 多项式与插值
? 来源于实际、又广泛用于实际。
? 多项式插值的主要目的是用一个多项式
拟合离散点上的函数值，使得可以用该
多项式估计数据点之间的函数值。
? 可导出数值积分方法，有限差分近似
? 关注插值多项式的表达式、精度、选点
效果。
2.1 关于多项式 MATLAB命令
? 一个多项式的幂级数形式可表示为,
? 也可表为嵌套形式
? 或因子形式
N阶多项式 n个根，其中包含重根和复根。若多
项式所有系数均为实数，则全部复根都将以共轭对
的形式出现
1
1 2 1
nn
nny c x c x c x c
?
?? ? ? ? ?L
1 2 3 1( ( ( ) ) )nny c x c x c x c x c ?? ? ? ? ?LL
1 1 2( ) ( ) ( )ny c x r x r x r? ? ? ?L
? 幂系数：在 MATLAB里，多项式用行向量表示，其
元素为多项式的系数，并从左至右按降幂排列。
例,
被表示为 >> p=[2 1 4 5]
>> poly2sym(p)
ans =
2*x^3+x^2+4*x+5
? Roots,多项式的零点可用命令 roots求的。
例,>> r=roots(p) 得到
r =
0.2500 + 1.5612i
0.2500 - 1.5612i
-1.0000
所有零点由一个列向量给出。
322 4 5y x x x? ? ? ?
? Poly,由零点可得原始多项式的各系数，但可能相差
一个常数倍。
例,>> poly(r)
ans =
1.0000 0.5000 2.0000 2.5000
注意：若存在重根，这种转换可能会降低精度。
例,
>> r=roots([1 -6 15 -20 15 -6 1])
r =
1.0042 + 0.0025i
1.0042 - 0.0025i
1.0000 + 0.0049i
1.0000 - 0.0049i
0.9958 + 0.0024i
0.9958 - 0.0024i
舍入误差的影响，与计算精度有关。
6 6 5 4 3 2( 1 ) 6 1 5 2 0 1 5 6 1y x x x x x x x? ? ? ? ? ? ? ? ?
? polyval,可用命令 polyval计算多项式的值。
例,计算 y(2.5)
>> c=[3,-7,2,1,1]； xi=2.5; yi=polyval(c,xi)
yi =
23.8125
如果 xi是含有多个横坐标值的数组，则 yi也
为与 xi长度相同的向量。
>> c=[3,-7,2,1,1]; xi=[2.5,3];
>> yi=polyval(c,xi)
yi =
23.8125 76.0000
4 3 23 7 2 1y x x x x? ? ? ? ?
? polyfit:给定 n+1个点将可以唯一确定一个 n阶多项式。利
用命令 polyfit可容易确定多项式的系数。
例,
>> x=[1.1,2.3,3.9,5.1];
>> y=[3.887,4.276,4.651,2.117];
>> a=polyfit(x,y,length(x)-1)
a =
-0.2015 1.4385 -2.7477 5.4370
>> poly2sym(a)
ans =
-403/2000*x^3+2877/2000*x^2-27477/10000*x+5437/1000
多项式为
Polyfit的第三个参数是多项式的阶数。
320, 2 0 1 5 1, 4 3 8 5 2, 7 4 7 7 5, 4 3 7 0y x x x? ? ? ? ?
多项式积分,
功能：求多项式积分
调用格式,py=poly_itg(p)
p:被积多项式的系数
py：求积后多项式的系数
poly_itg.m
function py=poly_itg(p)
n=length(p);
py=[p.*[n:-1:1].^(-1),0]
不包括最后一项积分常数
1212
1212
nn n
nn
cccY y d x x x x c x c
nn
?
??? ? ? ? ? ? ??? L
11 2 1nn nny c x c x c x c? ?? ? ? ? ?L
多项式微分,
? Polyder,求多项式一阶导数的系数。
调用格式为, b=polyder(c )
c为多项式 y的系数,b是微分后的系数,
其值为,
1
1 2 1
nn
nny c x c x c x c
?
?? ? ? ? ?L
' 1 2
12 ( 1 )
nn
ny n c x n c x c
??? ? ? ? ?L
12[,( 1 ),,]nn c n c c? L
两个多项式的和与差,
命令 poly_add:求两个多项式的和,其调用格式为,
c= poly_add(a,b)
多项式 a减去 b，可表示为,
c= poly_add(a,-b)
1
1 2 1
nn
b n ny b x b x b x b
?
?? ? ? ? ?L
1
1 2 1
mm
a m my a x a x a x a
?
?? ? ? ? ?L
功能：两个多项式相加
调用格式,b=poly_add(p1,p2)
b:求和后的系数数组
? poly_add.m
function p3=poly_add(p1,p2)
n1=length(p1);
n2=length(p2);
if n1==n2 p3=p1+p2;end
if n1>n2 p3=p1+[zeros(1,n1-n2),p2];end
if n1<n2 p3=[zeros(1,n2-n1),p1]+p2;end
? m阶多项式与 n阶多项式的乘积是 d=m+n阶的多项式,
计算 系数的 MATLAB命令是,c=conv(a,b)
? 多项式 除多项式 的除法满足,
其中 是商,是除法的余数。多项式 和
可由命令 deconv算出。
例,[q,r]=deconv(a,b)
1
1 2 1
dd
c a b d dy y y c x c x c x c
?
?? ? ? ? ? ?L
cy
by ay
a q b ry y y y??
qy ry qy
ry
11 2 1mma m my a x a x a x a? ?? ? ? ? ?L
11 2 1nnb n ny b x b x b x b? ?? ? ? ? ?L
? 例
>> a=[2,-5,6,-1,9]; b=[3,-90,-18];
>> c=conv(a,b)
c =
6 -195 432 -453 9 -792 -162
>> [q,r]=deconv(c,b)
q =
2 -5 6 -1 9
r =
0 0 0 0 0 0 0
>> poly2sym(c)
ans =
6*x^6-195*x^5+432*x^4-453*x^3+9*x^2-792*x-162
2.2 Lagrange插值
? 方法介绍
对给定的 n个插值点 及对应的函
数值,利用 n次 Lagrange插值多项式，
则对插值区间内任意 x的函数值 y可通过下式
求的,
? MATLAB实现
12,,,nx x xL
12,,,ny y yL
1 1
( ) ( )
nn
j
k
k j kj
jk
xx
y x y
xx? ?
?
?
?
?? ?
function y=lagrange(x0,y0,x)
ii=1:length(x0); y=zeros(size(x));
for i=ii
ij=find(ii~=i); y1=1;
for j=1:length(ij),y1=y1.*(x-x0(ij(j))); end
y=y+y1*y0(i)/prod(x0(i)-x0(ij));
end
? 算例：给出 f(x)=ln(x)的数值表，用 Lagrange计算
ln(0.54)的近似值。
>> x=[0.4:0.1:0.8]；
>> y=[-0.916291,-0.693147,-0.510826,-0.356675,-0.223144];
>> lagrange(x,y,0.54)
ans =
-0.6161 （精确解 -0.616143)
1 1
( ) ( )
nn
j
k
k j kj
jk
xx
y x y
xx? ?
?
?
?
?? ?
2.3 Hermite插值
? 方法介绍
不少实际问题不但要求在节点上函数值相等，而且
要求导数值也相等，甚至要求高阶导数值也相等，满足
这一要求的插值多项式就是 Hermite插值多项式。下面
只讨论函数值与一阶导数值个数相等且已知的情况。
已知 n个插值点 及对应的函数值
和一阶导数值 。则对插值区间
内任意 x的函数值 y的 Hermite插值公式,
12,,,nx x xL
12,,,ny y yL '' '12,,,ny y yL
'
1
2
11
( ) [ ( ) ( 2 ) ]
1
( ) ;
n
i i i i i i
i
n n
j
ii
jj i j i j
jiji
y x h x x a y y y
xx
ha
x x x x
?
??
??
? ? ? ?
?
??
??
?
??其中
? MATLAB实现
% hermite.m
function y=hermite(x0,y0,y1,x)
n=length(x0); m=length(x);
for k=1:m yy=0.0;
for i=1:n h=1.0; a=0.0;
for j=1:n
if j~=i
h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;
a=1/(x0(i)-x0(j))+a;
end
end
yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));
end
y(k)=yy;
end
'
1
2
11
( ) [ ( ) ( 2 ) ]
1
( ) ;
n
i i i i i i
i
n n
j
ii
jj i j i j
jiji
y x h x x a y y y
xx
ha
x x x x
?
??
??
? ? ? ?
?
??
??
?
??其中
? 算例,对给定数据,试构造 Hermite多项式求出
sin0.34的近似值。
>> x0=[0.3,0.32,0.35];
>> y0=[0.29552,0.31457,0.34290];
>> y1=[0.95534,0.94924,0.93937];
>> y=hermite(x0,y0,y1,0.34)
y =
0.3335
>> sin(0.34) ％与精确值比较
ans =
0.3335
>> x=[0.3:0.005:0.35];y=hermite(x0,y0,y1,x);
>> plot(x,y)
>> y2=sin(x); hold on
>> plot(x,y2,'--r')
2.4 Runge现象和分段插值
? 问题的提出：根据区间 [a,b]上给出的节点做
插值多项式 p(x)的近似值，一般总认为 p(x)的
次数越高则逼近 f(x)的精度就越好，但事实并
非如此。
? 反例,
在区间 [-5,5]上的各阶导数存在，但在此
区间上取 n个节点所构成的 Lagrange插值多项
式在全区间内并非都收敛。
? 取 n=10,用 Lagrange插值法进行插值计算。
2
1()
1fx x? ?
>> x=[-5:1:5]; y=1./(1+x.^2); x0=[-5:0.1:5];
>> y0=lagrange(x,y,x0);
>> y1=1./(1+x0.^2);
％绘制图形
>> plot(x0,y0,'--r')
％插值曲线
>> hold on
>> plot(x0,y1,‘-b')
％原曲线
? 为解决 Rung问题，引入分段插值。
? 算法分析：所谓分段插值就是通过插值点用折
线或低次曲线连接起来逼近原曲线。
? MATLAB实现 可调用内部函数。
– 命令 1 interp1
? 功能, 一维数据插值（表格查找）。该命令对数据点之
间计算内插值。它找出一元函数 f(x)在中间点的数值。其
中函数 f(x)由所给数据决定。
? 格式 1 yi = interp1(x,Y,xi)
%返回插值向量 yi，每一元素对应于参量 xi，同时由向
量 x与 Y的内插值决定。参量 x指定数据 Y的点。若 Y为一矩阵，
则按 Y的每列计算。
? 算例
对于 t,beta, alpha分别有两组数据与之对应，用分段线
性插值法计算当 t=321,440,571时 beta, alpha的值。
2.5 分段插值
>> temp=[300,400,500,600]';
>> beta=1000*[3.33,2.50,2.00,1.67]';
>> alpha=10000*[0.2128,0.3605,0.5324,0.7190]';
>> ti=[321,400,571]';
>> propty=interp1(temp,[beta,alpha],ti);
％ propty=interp1(temp,[beta,alpha],ti,’linear’);
>> [ti,propty]
ans =
1.0e+003 *
0.3210 3.1557 2.4382
0.4000 2.5000 3.6050
0.5710 1.7657 6.6489
? 格式 2 yi = interp1(Y,xi)
%假定 x=1:N,其中 N为向量 Y的长度, 或者为矩阵
Y的行数 。
? 格式 3 yi = interp1(x,Y,xi,method) %用指定的算法计算
插值,
’ nearest’：最近邻点插值, 直接完成计算；
’ linear’：线性插值 （ 缺省方式 ）, 直接完成计算；
’ spline’：三次样条函数插值 。
’ cubic’,分段三次 Hermite插值 。
其它, 如 ’ v5cubic’ 。
对于超出 x 范围的 xi 的分量, 使用方
法 ’ nearest’,’ linear’,’ v5cubic’的插值算法, 相应
地将返回 NaN。 对其他的方法, interp1将对超出的分
量执行外插值算法 。
? yi = interp1(x,Y,xi,method,'extrap')
? yi = interp1(x,Y,xi,method,extrapval) %确定超出 x范围
的 xi中的分量的外插值 extrapval,其值通常取 NaN或 0。
? 算例
>> year=1900:10:2010;
>> product=[75.995,91.972,105.711,123.203,131.669,..,
150.697,179.323,203.212,226.505,249.633,256.344,267.893];
>> p1995 = interp1(year,product,1995)
p1995 =
252.9885
>> x = 1900:1:2010;
>> y = interp1(year,…
product,x,'cubic');
>> plot(year,…
product,'o',x,y)
? 例,已知的数据点来自函数
根据生成的数据进行插值处理，得出较平滑的曲线
直接生成数据。
>> x=0:.12:1;
>> y=(x.^2-3*x…
+5).*exp(-5*x…
).*sin(x);
>> plot(x,y,x,y,'o')
? 调用 interp1( )函数,
>> x1=0:.02:1; y0=(x1.^2-3*x1+5).*exp(-5*x1).*sin(x1);
>> y1=interp1(x,y,x1); y2=interp1(x,y,x1,'cubic');
>> y3=interp1(x,y,x1,'spline'); y4=interp1(x,y,x1,'nearest');
>> plot(x1,[y1',y2',y3',y4'],':',x,y,'o',x1,y0)
? 误差分析
>> [max(abs(y0(1:49) …
-y2(1:49))),
max(abs(y0-y3)),
max(abs(y0-y4))]
ans =
0.0177 0.0086 0.1598
>> x0=-1+2*[0:10]/10;
>> y0=1./(1+25*x0.^2);
>> x=-1:.01:1;
>> y=lagrange(x0,y0,x);
% Lagrange 插值
>> ya=1./(1+25*x.^2);
>> plot(x,ya,x,y,':')
例
>> y1=interp1(x0,y0,x,'cubic'); y2=interp1(x0,y0,x,'spline');
>> plot(x,ya,x,y1,':',x,y2,'--')
– 命令 2 interp2
? 功能 二维数据内插值（表格查找）
? 格式 1 ZI = interp2(X,Y,Z,XI,YI)
%返回矩阵 ZI，其元素包含对应于参量 XI与
YI（可以是向量、或同型矩阵）的元素。参量 X
与 Y必须是单调的，且相同的划分格式，就像由
命令 meshgrid生成的一样。若 Xi与 Yi中有在 X与
Y范围之外的点，则相应地返回 NaN。
? 格式 2 ZI = interp2(Z,XI,YI) %缺省地,X=1:n、
Y=1:m，其中 [m,n]=size(Z)。再按第一种情形进
行计算。
? 格式 3 ZI = interp2(X,Y,Z,XI,YI,method) %用指
定的算法 method计算二维插值,
’ linear’：双线性插值算法 （ 缺省算法 ） ；
’ nearest’：最临近插值；
’ spline’：三次样条插值；
’ cubic’：双三次插值 。
? 算例,
>> years=1950:10:1990;
>> service=10:10:30;
>> wage = [150.697 199.592 187.625
179.323 195.072 250.287
203.212 179.092 322.767
226.505 153.706 426.730
249.633 120.281 598.243];
>> w = interp2(service,years,wage,15,1975)
w =
190.6288
? 例
>> [x,y]=meshgrid(-3:.6:3,-2:.4:2);
>> z=(x.^2-2*x).*exp…
(-x.^2-y.^2-x.*y);
>> surf(x,y,z),
axis([-3,3,-2,2,-0.7,1.5])
? 选较密的插值点，用默认的线性插值算
法进行插值
>> [x1,y1]=meshgrid(-3:.2:3,-2:.2:2);
>> z1=interp2(x,y,z,x1,y1);
>> surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])
? 立方和样条插值,
>> z1=interp2(x,y,z,x1,y1,'cubic');
>> z2=interp2(x,y,z,x1,y1,'spline');
>> surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])
>> figure;surf(x1,y1,z2),axis([-3,3,-2,2,-0.7,1.5])
? 算法误差的比较
> z=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
>> surf(x1,y1,abs(z-z1)),axis([-3,3,-2,2,0,0.08])
>> figure;surf(x1,y1,abs(z-z2)),axis([-3,3,-2,2,0,0.025])
二维一般分布数据的插值
? 功能,可对非网格数据进行插值
? 格式,z=griddata(x0,y0,z0,x,y,’method’)
’ v4 ’:MATLAB4.0提供的插值算法，公认效果较好；
’ linear’：双线性插值算法（缺省算法）；
’ nearest’：最临近插值；
’ spline’：三次样条插值；
’ cubic’：双三次插值。
? 例,
在 x为［ 3,3］,y为［－ 2,2］矩形区域随机
选择一组坐标，用 ’ v4 ’与’ cubic’插值法进
行处理，并对误差进行比较。
222(,) ( 2 ) x y x yz f x y x x e ? ? ?? ? ?
>> x=-3+6*rand(200,1);y=-2+4*rand(200,1);
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> [x1,y1]=meshgrid(-3:.2:3,-2:.2:2);
>> z1=griddata(x,y,z,x1,y1,'cubic');
>> surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])
>> z2=griddata(x,y,z,x1,y1,'v4');
>> figure;surf(x1,y1,z2),axis([-3,3,-2,2,-0.7,1.5])
误差分析
>> z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
>> surf(x1,y1,abs(z0-z1)),axis([-3,3,-2,2,0,0.15])
>> figure;surf(x1,y1,abs(z0-z2)),axis([-3,3,-2,2,0,0.15])
? 例,
在 x为［ 3,3］,y为［－ 2,2］矩形区域随机选择一组坐标中，
对分布不均匀数据，进行插值分析。
>> x=-3+6*rand(200,1); y=-2+4*rand(200,1);
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); % 生成已知数据
>> plot(x,y,'x') % 样本点的二维分布
>> figure,plot3(x,y,z,'x'),axis([-3,3,-2,2,-0.7,1.5]),grid
222(,) ( 2 ) x y x yz f x y x x e ? ? ?? ? ?
? 去除在 (-1,-1/2)点为圆心，以 0.5为半径的圆内的点。
>> x=-3+6*rand(200,1); y=-2+4*rand(200,1); % 重新生成样本点
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> ii=find((x+1).^2+(y+0.5).^2>0.5^2); % 找出不满足条件的点坐标
>> x=x(ii); y=y(ii); z=z(ii); plot(x,y,'x')
>> t=[0:.1:2*pi,2*pi]; x0=-1+0.5*cos(t); y0=-0.5+0.5*sin(t);
>> line(x0,y0) % 在图形上叠印该圆，可见，圆内样本点均已剔除
? 用新样本点拟合出曲面
>> [x1,y1]=meshgrid(-3:.2:3,-2:.2:2);
>> z1=griddata(x,y,z,x1,y1,'v4');
>> surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])
? 误差分析
>> z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
>> surf(x1,y1,abs(z0-z1)),axis([-3,3,-2,2,0,0.1])
>> contour(x1,y1,abs(z0-z1),30); hold on,plot(x,y,‘x’);
line(x0,y0) ％误差的二维等高线图
– 命令 3 interp3
三维网格生成用 meshgrid( )函数,调用格式,
[x,y,z]=meshgrid(x1,y1,z1)
其中 x1,y1,z1为这三维所需要的分割形式，应以向
量形式给出，返回 x,y,z为网格的数据生成，均为
三维数组。
griddata3( ) 三维非网格形式的插值拟合
– 命令 4 interpn
n维网格生成用 ndgrid( )函数,调用格式,
[x1,x2,…,xn]=ndgrid[v1,v2,…,vn]
griddatan( ) n维非网格形式的插值拟合
interp3 ( ),interpn( )调用格式同 interp2()函数一致；
griddata3( ),griddatan( )调用格式同 griddata( )函数
一致。
? 例,
通过函数生成一些网格型样本点，试根据样
本点进行拟合，并给出拟合误差。
>> [x,y,z]=meshgrid(-1:0.2:1); [x0,y0,z0]=meshgrid(-1:0.05:1);
>> V=exp(x.^2.*z+y.^2.*x+z.^2.*y…
).*cos(x.^2.*y.*z+z.^2.*y.*x);
>> V0=exp(x0.^2.*z0+y0.^2.*x0…
+z0.^2.*y0).*cos(x0.^2.*y0.*z0+z0.^2.*y0.*x0);
>> V1=interp3(x,y,z,V,x0,y0,z0,'spline'); err=V1-V0; max(err(:))
ans =
0.0419
2 2 2 22(,,) c o s ( )x z y x z yV x y z e x y z z y x????
2.6样条插值的 MATLAB表示
定义一个三次样条函数类,
S=csapi(x,y)
其中 x=[x1,x2,….,xn],y=[y1,y2,…,yn] 为样本
点。 S返回样条函数对象的插值结果，包括子
区间点、各区间点三次多项式系数等。
? 可用 fnplt()绘制出插值结果，其调用格式,
fnplt(S)
? 对给定的向量 xp，可用 fnval()函数计算，其调
用格式,
yp=fnval(S,xp)
其中得出的 yp是 xp上各点的插值结果。
? 例,
>> x0=[0,0.4,1,2,pi]; y0=sin(x0);
>> sp=csapi(x0,y0),fnplt(sp,':'); hold on,
sp =
form,'pp'
breaks,[0 0.4000 1 2 3.1416]
coefs,[4x4 double]
pieces,4
order,4
dim,1
>> ezplot('sin(t)',[0,pi]); plot(x0,y0,'o')
>> sp.coefs
ans =
-0.1627 0.0076 0.9965 0
-0.1627 -0.1876 0.9245 0.3894
0.0244 -0.4804 0.5238 0.8415
0.0244 -0.4071 -0.3637 0.9093
在 (0.4000,1)区间内，插值多项式可以表示为,
32
2 ( ) 0, 1 6 2 7 ( 0, 4 ) 0, 1 8 7 6 ( 0, 4 )
0, 9 2 4 5 ( 0, 4 ) 0, 3 8 9 4
S x x x
x
? ? ? ? ?
? ? ?
例
点，用三次样条插值的方法对这些数据进行拟合
>> x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x);
>> sp=csapi(x,y); fnplt(sp)
>> c=[sp.breaks(1:4)' sp.breaks(2:5)' sp.coefs(1:4,:),..,
sp.breaks(5:8)' sp.breaks(6:9)' sp.coefs(5:8,:) ]
c =
Columns 1 through 7
0 0.1200 24.7396 -19.3588 4.5151 0 0.4800
0.1200 0.2400 24.7396 -10.4526 0.9377 0.3058 0.6000
0.2400 0.3600 4.5071 -1.5463 -0.5022 0.3105 0.7200
0.3600 0.4800 1.9139 0.0762 -0.6786 0.2358 0.8400
Columns 8 through 12
0.6000 -0.2404 0.7652 -0.5776 0.1588
0.7200 -0.4774 0.6787 -0.4043 0.1001
0.8400 -0.4559 0.5068 -0.2621 0.0605
0.9600 -0.4559 0.3427 -0.1601 0.0356
? 格式
S=csapi({x1,x2,…,x n},z)
处理多个自变量的网格数据三次样条插值类,
>> x0=-3:.6:3; y0=-2:.4:2; [x,y]=ndgrid(x0,y0); % 注意这里
只能用 ndgrid,
否则生成的 z 矩阵
顺序有问题
>> z=(x.^2-2*x).*…
exp(-x.^2-y.^2-x.*y);
>> sp=csapi({x0,y0},z);
>> fnplt(sp);
例
函数 spline
– 功能 三次样条数据插值
– 格式 yy = spline(x,y,xx)
– 例,对离散分布在 y=exp(x)sin(x)函数曲线上的数
据点进行样条插值计算,
>>x = [0 2 4 5 8 12 12.8 17.2 19.9 20];
>> y = exp(x).*sin(x);
>>xx = 0:.25:20;
>>yy = spline(x,y,xx);
>>plot(x,y,'o',xx,yy)
? 来源于实际、又广泛用于实际。
? 多项式插值的主要目的是用一个多项式
拟合离散点上的函数值，使得可以用该
多项式估计数据点之间的函数值。
? 可导出数值积分方法，有限差分近似
? 关注插值多项式的表达式、精度、选点
效果。
2.1 关于多项式 MATLAB命令
? 一个多项式的幂级数形式可表示为,
? 也可表为嵌套形式
? 或因子形式
N阶多项式 n个根，其中包含重根和复根。若多
项式所有系数均为实数，则全部复根都将以共轭对
的形式出现
1
1 2 1
nn
nny c x c x c x c
?
?? ? ? ? ?L
1 2 3 1( ( ( ) ) )nny c x c x c x c x c ?? ? ? ? ?LL
1 1 2( ) ( ) ( )ny c x r x r x r? ? ? ?L
? 幂系数：在 MATLAB里，多项式用行向量表示，其
元素为多项式的系数，并从左至右按降幂排列。
例,
被表示为 >> p=[2 1 4 5]
>> poly2sym(p)
ans =
2*x^3+x^2+4*x+5
? Roots,多项式的零点可用命令 roots求的。
例,>> r=roots(p) 得到
r =
0.2500 + 1.5612i
0.2500 - 1.5612i
-1.0000
所有零点由一个列向量给出。
322 4 5y x x x? ? ? ?
? Poly,由零点可得原始多项式的各系数，但可能相差
一个常数倍。
例,>> poly(r)
ans =
1.0000 0.5000 2.0000 2.5000
注意：若存在重根，这种转换可能会降低精度。
例,
>> r=roots([1 -6 15 -20 15 -6 1])
r =
1.0042 + 0.0025i
1.0042 - 0.0025i
1.0000 + 0.0049i
1.0000 - 0.0049i
0.9958 + 0.0024i
0.9958 - 0.0024i
舍入误差的影响，与计算精度有关。
6 6 5 4 3 2( 1 ) 6 1 5 2 0 1 5 6 1y x x x x x x x? ? ? ? ? ? ? ? ?
? polyval,可用命令 polyval计算多项式的值。
例,计算 y(2.5)
>> c=[3,-7,2,1,1]； xi=2.5; yi=polyval(c,xi)
yi =
23.8125
如果 xi是含有多个横坐标值的数组，则 yi也
为与 xi长度相同的向量。
>> c=[3,-7,2,1,1]; xi=[2.5,3];
>> yi=polyval(c,xi)
yi =
23.8125 76.0000
4 3 23 7 2 1y x x x x? ? ? ? ?
? polyfit:给定 n+1个点将可以唯一确定一个 n阶多项式。利
用命令 polyfit可容易确定多项式的系数。
例,
>> x=[1.1,2.3,3.9,5.1];
>> y=[3.887,4.276,4.651,2.117];
>> a=polyfit(x,y,length(x)-1)
a =
-0.2015 1.4385 -2.7477 5.4370
>> poly2sym(a)
ans =
-403/2000*x^3+2877/2000*x^2-27477/10000*x+5437/1000
多项式为
Polyfit的第三个参数是多项式的阶数。
320, 2 0 1 5 1, 4 3 8 5 2, 7 4 7 7 5, 4 3 7 0y x x x? ? ? ? ?
多项式积分,
功能：求多项式积分
调用格式,py=poly_itg(p)
p:被积多项式的系数
py：求积后多项式的系数
poly_itg.m
function py=poly_itg(p)
n=length(p);
py=[p.*[n:-1:1].^(-1),0]
不包括最后一项积分常数
1212
1212
nn n
nn
cccY y d x x x x c x c
nn
?
??? ? ? ? ? ? ??? L
11 2 1nn nny c x c x c x c? ?? ? ? ? ?L
多项式微分,
? Polyder,求多项式一阶导数的系数。
调用格式为, b=polyder(c )
c为多项式 y的系数,b是微分后的系数,
其值为,
1
1 2 1
nn
nny c x c x c x c
?
?? ? ? ? ?L
' 1 2
12 ( 1 )
nn
ny n c x n c x c
??? ? ? ? ?L
12[,( 1 ),,]nn c n c c? L
两个多项式的和与差,
命令 poly_add:求两个多项式的和,其调用格式为,
c= poly_add(a,b)
多项式 a减去 b，可表示为,
c= poly_add(a,-b)
1
1 2 1
nn
b n ny b x b x b x b
?
?? ? ? ? ?L
1
1 2 1
mm
a m my a x a x a x a
?
?? ? ? ? ?L
功能：两个多项式相加
调用格式,b=poly_add(p1,p2)
b:求和后的系数数组
? poly_add.m
function p3=poly_add(p1,p2)
n1=length(p1);
n2=length(p2);
if n1==n2 p3=p1+p2;end
if n1>n2 p3=p1+[zeros(1,n1-n2),p2];end
if n1<n2 p3=[zeros(1,n2-n1),p1]+p2;end
? m阶多项式与 n阶多项式的乘积是 d=m+n阶的多项式,
计算 系数的 MATLAB命令是,c=conv(a,b)
? 多项式 除多项式 的除法满足,
其中 是商,是除法的余数。多项式 和
可由命令 deconv算出。
例,[q,r]=deconv(a,b)
1
1 2 1
dd
c a b d dy y y c x c x c x c
?
?? ? ? ? ? ?L
cy
by ay
a q b ry y y y??
qy ry qy
ry
11 2 1mma m my a x a x a x a? ?? ? ? ? ?L
11 2 1nnb n ny b x b x b x b? ?? ? ? ? ?L
? 例
>> a=[2,-5,6,-1,9]; b=[3,-90,-18];
>> c=conv(a,b)
c =
6 -195 432 -453 9 -792 -162
>> [q,r]=deconv(c,b)
q =
2 -5 6 -1 9
r =
0 0 0 0 0 0 0
>> poly2sym(c)
ans =
6*x^6-195*x^5+432*x^4-453*x^3+9*x^2-792*x-162
2.2 Lagrange插值
? 方法介绍
对给定的 n个插值点 及对应的函
数值,利用 n次 Lagrange插值多项式，
则对插值区间内任意 x的函数值 y可通过下式
求的,
? MATLAB实现
12,,,nx x xL
12,,,ny y yL
1 1
( ) ( )
nn
j
k
k j kj
jk
xx
y x y
xx? ?
?
?
?
?? ?
function y=lagrange(x0,y0,x)
ii=1:length(x0); y=zeros(size(x));
for i=ii
ij=find(ii~=i); y1=1;
for j=1:length(ij),y1=y1.*(x-x0(ij(j))); end
y=y+y1*y0(i)/prod(x0(i)-x0(ij));
end
? 算例：给出 f(x)=ln(x)的数值表，用 Lagrange计算
ln(0.54)的近似值。
>> x=[0.4:0.1:0.8]；
>> y=[-0.916291,-0.693147,-0.510826,-0.356675,-0.223144];
>> lagrange(x,y,0.54)
ans =
-0.6161 （精确解 -0.616143)
1 1
( ) ( )
nn
j
k
k j kj
jk
xx
y x y
xx? ?
?
?
?
?? ?
2.3 Hermite插值
? 方法介绍
不少实际问题不但要求在节点上函数值相等，而且
要求导数值也相等，甚至要求高阶导数值也相等，满足
这一要求的插值多项式就是 Hermite插值多项式。下面
只讨论函数值与一阶导数值个数相等且已知的情况。
已知 n个插值点 及对应的函数值
和一阶导数值 。则对插值区间
内任意 x的函数值 y的 Hermite插值公式,
12,,,nx x xL
12,,,ny y yL '' '12,,,ny y yL
'
1
2
11
( ) [ ( ) ( 2 ) ]
1
( ) ;
n
i i i i i i
i
n n
j
ii
jj i j i j
jiji
y x h x x a y y y
xx
ha
x x x x
?
??
??
? ? ? ?
?
??
??
?
??其中
? MATLAB实现
% hermite.m
function y=hermite(x0,y0,y1,x)
n=length(x0); m=length(x);
for k=1:m yy=0.0;
for i=1:n h=1.0; a=0.0;
for j=1:n
if j~=i
h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;
a=1/(x0(i)-x0(j))+a;
end
end
yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));
end
y(k)=yy;
end
'
1
2
11
( ) [ ( ) ( 2 ) ]
1
( ) ;
n
i i i i i i
i
n n
j
ii
jj i j i j
jiji
y x h x x a y y y
xx
ha
x x x x
?
??
??
? ? ? ?
?
??
??
?
??其中
? 算例,对给定数据,试构造 Hermite多项式求出
sin0.34的近似值。
>> x0=[0.3,0.32,0.35];
>> y0=[0.29552,0.31457,0.34290];
>> y1=[0.95534,0.94924,0.93937];
>> y=hermite(x0,y0,y1,0.34)
y =
0.3335
>> sin(0.34) ％与精确值比较
ans =
0.3335
>> x=[0.3:0.005:0.35];y=hermite(x0,y0,y1,x);
>> plot(x,y)
>> y2=sin(x); hold on
>> plot(x,y2,'--r')
2.4 Runge现象和分段插值
? 问题的提出：根据区间 [a,b]上给出的节点做
插值多项式 p(x)的近似值，一般总认为 p(x)的
次数越高则逼近 f(x)的精度就越好，但事实并
非如此。
? 反例,
在区间 [-5,5]上的各阶导数存在，但在此
区间上取 n个节点所构成的 Lagrange插值多项
式在全区间内并非都收敛。
? 取 n=10,用 Lagrange插值法进行插值计算。
2
1()
1fx x? ?
>> x=[-5:1:5]; y=1./(1+x.^2); x0=[-5:0.1:5];
>> y0=lagrange(x,y,x0);
>> y1=1./(1+x0.^2);
％绘制图形
>> plot(x0,y0,'--r')
％插值曲线
>> hold on
>> plot(x0,y1,‘-b')
％原曲线
? 为解决 Rung问题，引入分段插值。
? 算法分析：所谓分段插值就是通过插值点用折
线或低次曲线连接起来逼近原曲线。
? MATLAB实现 可调用内部函数。
– 命令 1 interp1
? 功能, 一维数据插值（表格查找）。该命令对数据点之
间计算内插值。它找出一元函数 f(x)在中间点的数值。其
中函数 f(x)由所给数据决定。
? 格式 1 yi = interp1(x,Y,xi)
%返回插值向量 yi，每一元素对应于参量 xi，同时由向
量 x与 Y的内插值决定。参量 x指定数据 Y的点。若 Y为一矩阵，
则按 Y的每列计算。
? 算例
对于 t,beta, alpha分别有两组数据与之对应，用分段线
性插值法计算当 t=321,440,571时 beta, alpha的值。
2.5 分段插值
>> temp=[300,400,500,600]';
>> beta=1000*[3.33,2.50,2.00,1.67]';
>> alpha=10000*[0.2128,0.3605,0.5324,0.7190]';
>> ti=[321,400,571]';
>> propty=interp1(temp,[beta,alpha],ti);
％ propty=interp1(temp,[beta,alpha],ti,’linear’);
>> [ti,propty]
ans =
1.0e+003 *
0.3210 3.1557 2.4382
0.4000 2.5000 3.6050
0.5710 1.7657 6.6489
? 格式 2 yi = interp1(Y,xi)
%假定 x=1:N,其中 N为向量 Y的长度, 或者为矩阵
Y的行数 。
? 格式 3 yi = interp1(x,Y,xi,method) %用指定的算法计算
插值,
’ nearest’：最近邻点插值, 直接完成计算；
’ linear’：线性插值 （ 缺省方式 ）, 直接完成计算；
’ spline’：三次样条函数插值 。
’ cubic’,分段三次 Hermite插值 。
其它, 如 ’ v5cubic’ 。
对于超出 x 范围的 xi 的分量, 使用方
法 ’ nearest’,’ linear’,’ v5cubic’的插值算法, 相应
地将返回 NaN。 对其他的方法, interp1将对超出的分
量执行外插值算法 。
? yi = interp1(x,Y,xi,method,'extrap')
? yi = interp1(x,Y,xi,method,extrapval) %确定超出 x范围
的 xi中的分量的外插值 extrapval,其值通常取 NaN或 0。
? 算例
>> year=1900:10:2010;
>> product=[75.995,91.972,105.711,123.203,131.669,..,
150.697,179.323,203.212,226.505,249.633,256.344,267.893];
>> p1995 = interp1(year,product,1995)
p1995 =
252.9885
>> x = 1900:1:2010;
>> y = interp1(year,…
product,x,'cubic');
>> plot(year,…
product,'o',x,y)
? 例,已知的数据点来自函数
根据生成的数据进行插值处理，得出较平滑的曲线
直接生成数据。
>> x=0:.12:1;
>> y=(x.^2-3*x…
+5).*exp(-5*x…
).*sin(x);
>> plot(x,y,x,y,'o')
? 调用 interp1( )函数,
>> x1=0:.02:1; y0=(x1.^2-3*x1+5).*exp(-5*x1).*sin(x1);
>> y1=interp1(x,y,x1); y2=interp1(x,y,x1,'cubic');
>> y3=interp1(x,y,x1,'spline'); y4=interp1(x,y,x1,'nearest');
>> plot(x1,[y1',y2',y3',y4'],':',x,y,'o',x1,y0)
? 误差分析
>> [max(abs(y0(1:49) …
-y2(1:49))),
max(abs(y0-y3)),
max(abs(y0-y4))]
ans =
0.0177 0.0086 0.1598
>> x0=-1+2*[0:10]/10;
>> y0=1./(1+25*x0.^2);
>> x=-1:.01:1;
>> y=lagrange(x0,y0,x);
% Lagrange 插值
>> ya=1./(1+25*x.^2);
>> plot(x,ya,x,y,':')
例
>> y1=interp1(x0,y0,x,'cubic'); y2=interp1(x0,y0,x,'spline');
>> plot(x,ya,x,y1,':',x,y2,'--')
– 命令 2 interp2
? 功能 二维数据内插值（表格查找）
? 格式 1 ZI = interp2(X,Y,Z,XI,YI)
%返回矩阵 ZI，其元素包含对应于参量 XI与
YI（可以是向量、或同型矩阵）的元素。参量 X
与 Y必须是单调的，且相同的划分格式，就像由
命令 meshgrid生成的一样。若 Xi与 Yi中有在 X与
Y范围之外的点，则相应地返回 NaN。
? 格式 2 ZI = interp2(Z,XI,YI) %缺省地,X=1:n、
Y=1:m，其中 [m,n]=size(Z)。再按第一种情形进
行计算。
? 格式 3 ZI = interp2(X,Y,Z,XI,YI,method) %用指
定的算法 method计算二维插值,
’ linear’：双线性插值算法 （ 缺省算法 ） ；
’ nearest’：最临近插值；
’ spline’：三次样条插值；
’ cubic’：双三次插值 。
? 算例,
>> years=1950:10:1990;
>> service=10:10:30;
>> wage = [150.697 199.592 187.625
179.323 195.072 250.287
203.212 179.092 322.767
226.505 153.706 426.730
249.633 120.281 598.243];
>> w = interp2(service,years,wage,15,1975)
w =
190.6288
? 例
>> [x,y]=meshgrid(-3:.6:3,-2:.4:2);
>> z=(x.^2-2*x).*exp…
(-x.^2-y.^2-x.*y);
>> surf(x,y,z),
axis([-3,3,-2,2,-0.7,1.5])
? 选较密的插值点，用默认的线性插值算
法进行插值
>> [x1,y1]=meshgrid(-3:.2:3,-2:.2:2);
>> z1=interp2(x,y,z,x1,y1);
>> surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])
? 立方和样条插值,
>> z1=interp2(x,y,z,x1,y1,'cubic');
>> z2=interp2(x,y,z,x1,y1,'spline');
>> surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])
>> figure;surf(x1,y1,z2),axis([-3,3,-2,2,-0.7,1.5])
? 算法误差的比较
> z=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
>> surf(x1,y1,abs(z-z1)),axis([-3,3,-2,2,0,0.08])
>> figure;surf(x1,y1,abs(z-z2)),axis([-3,3,-2,2,0,0.025])
二维一般分布数据的插值
? 功能,可对非网格数据进行插值
? 格式,z=griddata(x0,y0,z0,x,y,’method’)
’ v4 ’:MATLAB4.0提供的插值算法，公认效果较好；
’ linear’：双线性插值算法（缺省算法）；
’ nearest’：最临近插值；
’ spline’：三次样条插值；
’ cubic’：双三次插值。
? 例,
在 x为［ 3,3］,y为［－ 2,2］矩形区域随机
选择一组坐标，用 ’ v4 ’与’ cubic’插值法进
行处理，并对误差进行比较。
222(,) ( 2 ) x y x yz f x y x x e ? ? ?? ? ?
>> x=-3+6*rand(200,1);y=-2+4*rand(200,1);
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> [x1,y1]=meshgrid(-3:.2:3,-2:.2:2);
>> z1=griddata(x,y,z,x1,y1,'cubic');
>> surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])
>> z2=griddata(x,y,z,x1,y1,'v4');
>> figure;surf(x1,y1,z2),axis([-3,3,-2,2,-0.7,1.5])
误差分析
>> z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
>> surf(x1,y1,abs(z0-z1)),axis([-3,3,-2,2,0,0.15])
>> figure;surf(x1,y1,abs(z0-z2)),axis([-3,3,-2,2,0,0.15])
? 例,
在 x为［ 3,3］,y为［－ 2,2］矩形区域随机选择一组坐标中，
对分布不均匀数据，进行插值分析。
>> x=-3+6*rand(200,1); y=-2+4*rand(200,1);
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); % 生成已知数据
>> plot(x,y,'x') % 样本点的二维分布
>> figure,plot3(x,y,z,'x'),axis([-3,3,-2,2,-0.7,1.5]),grid
222(,) ( 2 ) x y x yz f x y x x e ? ? ?? ? ?
? 去除在 (-1,-1/2)点为圆心，以 0.5为半径的圆内的点。
>> x=-3+6*rand(200,1); y=-2+4*rand(200,1); % 重新生成样本点
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> ii=find((x+1).^2+(y+0.5).^2>0.5^2); % 找出不满足条件的点坐标
>> x=x(ii); y=y(ii); z=z(ii); plot(x,y,'x')
>> t=[0:.1:2*pi,2*pi]; x0=-1+0.5*cos(t); y0=-0.5+0.5*sin(t);
>> line(x0,y0) % 在图形上叠印该圆，可见，圆内样本点均已剔除
? 用新样本点拟合出曲面
>> [x1,y1]=meshgrid(-3:.2:3,-2:.2:2);
>> z1=griddata(x,y,z,x1,y1,'v4');
>> surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])
? 误差分析
>> z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
>> surf(x1,y1,abs(z0-z1)),axis([-3,3,-2,2,0,0.1])
>> contour(x1,y1,abs(z0-z1),30); hold on,plot(x,y,‘x’);
line(x0,y0) ％误差的二维等高线图
– 命令 3 interp3
三维网格生成用 meshgrid( )函数,调用格式,
[x,y,z]=meshgrid(x1,y1,z1)
其中 x1,y1,z1为这三维所需要的分割形式，应以向
量形式给出，返回 x,y,z为网格的数据生成，均为
三维数组。
griddata3( ) 三维非网格形式的插值拟合
– 命令 4 interpn
n维网格生成用 ndgrid( )函数,调用格式,
[x1,x2,…,xn]=ndgrid[v1,v2,…,vn]
griddatan( ) n维非网格形式的插值拟合
interp3 ( ),interpn( )调用格式同 interp2()函数一致；
griddata3( ),griddatan( )调用格式同 griddata( )函数
一致。
? 例,
通过函数生成一些网格型样本点，试根据样
本点进行拟合，并给出拟合误差。
>> [x,y,z]=meshgrid(-1:0.2:1); [x0,y0,z0]=meshgrid(-1:0.05:1);
>> V=exp(x.^2.*z+y.^2.*x+z.^2.*y…
).*cos(x.^2.*y.*z+z.^2.*y.*x);
>> V0=exp(x0.^2.*z0+y0.^2.*x0…
+z0.^2.*y0).*cos(x0.^2.*y0.*z0+z0.^2.*y0.*x0);
>> V1=interp3(x,y,z,V,x0,y0,z0,'spline'); err=V1-V0; max(err(:))
ans =
0.0419
2 2 2 22(,,) c o s ( )x z y x z yV x y z e x y z z y x????
2.6样条插值的 MATLAB表示
定义一个三次样条函数类,
S=csapi(x,y)
其中 x=[x1,x2,….,xn],y=[y1,y2,…,yn] 为样本
点。 S返回样条函数对象的插值结果，包括子
区间点、各区间点三次多项式系数等。
? 可用 fnplt()绘制出插值结果，其调用格式,
fnplt(S)
? 对给定的向量 xp，可用 fnval()函数计算，其调
用格式,
yp=fnval(S,xp)
其中得出的 yp是 xp上各点的插值结果。
? 例,
>> x0=[0,0.4,1,2,pi]; y0=sin(x0);
>> sp=csapi(x0,y0),fnplt(sp,':'); hold on,
sp =
form,'pp'
breaks,[0 0.4000 1 2 3.1416]
coefs,[4x4 double]
pieces,4
order,4
dim,1
>> ezplot('sin(t)',[0,pi]); plot(x0,y0,'o')
>> sp.coefs
ans =
-0.1627 0.0076 0.9965 0
-0.1627 -0.1876 0.9245 0.3894
0.0244 -0.4804 0.5238 0.8415
0.0244 -0.4071 -0.3637 0.9093
在 (0.4000,1)区间内，插值多项式可以表示为,
32
2 ( ) 0, 1 6 2 7 ( 0, 4 ) 0, 1 8 7 6 ( 0, 4 )
0, 9 2 4 5 ( 0, 4 ) 0, 3 8 9 4
S x x x
x
? ? ? ? ?
? ? ?
例
点，用三次样条插值的方法对这些数据进行拟合
>> x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x);
>> sp=csapi(x,y); fnplt(sp)
>> c=[sp.breaks(1:4)' sp.breaks(2:5)' sp.coefs(1:4,:),..,
sp.breaks(5:8)' sp.breaks(6:9)' sp.coefs(5:8,:) ]
c =
Columns 1 through 7
0 0.1200 24.7396 -19.3588 4.5151 0 0.4800
0.1200 0.2400 24.7396 -10.4526 0.9377 0.3058 0.6000
0.2400 0.3600 4.5071 -1.5463 -0.5022 0.3105 0.7200
0.3600 0.4800 1.9139 0.0762 -0.6786 0.2358 0.8400
Columns 8 through 12
0.6000 -0.2404 0.7652 -0.5776 0.1588
0.7200 -0.4774 0.6787 -0.4043 0.1001
0.8400 -0.4559 0.5068 -0.2621 0.0605
0.9600 -0.4559 0.3427 -0.1601 0.0356
? 格式
S=csapi({x1,x2,…,x n},z)
处理多个自变量的网格数据三次样条插值类,
>> x0=-3:.6:3; y0=-2:.4:2; [x,y]=ndgrid(x0,y0); % 注意这里
只能用 ndgrid,
否则生成的 z 矩阵
顺序有问题
>> z=(x.^2-2*x).*…
exp(-x.^2-y.^2-x.*y);
>> sp=csapi({x0,y0},z);
>> fnplt(sp);
例
函数 spline
– 功能 三次样条数据插值
– 格式 yy = spline(x,y,xx)
– 例,对离散分布在 y=exp(x)sin(x)函数曲线上的数
据点进行样条插值计算,
>>x = [0 2 4 5 8 12 12.8 17.2 19.9 20];
>> y = exp(x).*sin(x);
>>xx = 0:.25:20;
>>yy = spline(x,y,xx);
>>plot(x,y,'o',xx,yy)
