matlab函数、插值和曲线拟合分析
下载
第10章 函数、插值和曲线拟合分析
MATLAB可以处理有估值和没有估值的多项式,还有一些强大的数值分析命令,如求零
值和最小值。 MATLAB中还有数据集合的插值、曲线拟合的命令和函数,还提到了经典的贝
赛尔(Bessel)函数。
10.1 MATLAB中的多项式
MATLAB将阶为n的多项式p(x)存储在长度为 n+1的行向量p中。元素为多项式的系数,并
按x的幂降序排列,表示为:
代表的多项式为:
令A是一个稀疏矩阵,向量 p和q的长度分别为 n+1和m+1,它们分别表示次数为 n和m的多
项式。MATLAB中处理多项式的命令如下:
多项式
polyval(p,x) 计算多项式 p。如果 x是一个标量,则计算出多项式在 x点 的值;如果 x是一个向量或者一个矩阵,则计算出多项式
在x中所有元素上的值。 计算向量x的多项式p的值。同上,计算结果在 y中,同时 [y,err]=
polyval(p,x,E) 还根据polyfit命令给出的矩阵 E返回一个误差估计向量
err。见help polyval和help polyfit,参见10.4节。 直接对矩阵A进行多项式计算。不是象上个命令一样对每个 polyvalm(p,A) 元素进行多项式计算,而是计算p(A)=p1An+p2An,1+…+pn+1I。 计算矩阵A的特征多项式向量。 poly(A)
poly(x) 给出一个长度为n+1的向量,其中的元素是次数为 n的多项
式的系数。这个多项式的根是长度为 n的向量x中元素。 计算带有系数 p的多项式的友矩阵 A,这个矩阵的特征多 compan(p) 项式为p。 计算特征多项式p的根,是一个长度为n的向量,也就是方程 roots(p) p(x)=0的解。表达式poly(roots(p))=p为真。结果可以是复数。 计算多项式p和q的乘积,也可以认为是 p和q的卷积。 conv(p,q)
[k,r]=deconv(p,q) 计算多项式 p除q。k是商多项式, r是残数多项式。这个
计算等价于p和q的逆卷积。
MATLAB 5 手册 136 下载 计算p(x)/q(x)的部分展开式: [u v k]=
residue(p,q) p(x) u(1) u(2) u(j) , , , k(x). ,, , , , q(x) x n(1) x n(2) x n(j)
向量p和q分别是多项式 p(x)和q(x)的系数。得到的余数在
向量u中,极点在列向量v中,商多项式在向量 k中。 [p q]=residue(u,v,x)从同上的部分展开式 u、v和x计算得到多项式p和q。 mpoles 给出极点多样性的相关信息,见help mpole。 s polyder(p) 计算得到长度为的微分多项式向量,多项式的系数在向量中。 n p
polyder(p,q) 返回一个向量,它表示由conv(p,q)定义的多项式微分。 [u,v]=polyder(p,q) 返回两个向量,它们表示由 deconv(p,q)定义的多项式
微分,表达形式为u/v。
? 例10.1
将一些多项式命令应用在下列的多项式上:
MATLAB用下列向量来表示这两个多项式:
(a) 计算多项式在x=1处的值,输入:
结果为:
(b) 很容易对向量或者矩阵计算多项式的值,输入:
结果为:
第10章 函数、插值和曲线拟合分析 137 下载
(c) 两个多项式相乘,得到一个新的多项式:
p5=conv(p2, p3)
给出:
(d) 用roots命令求多项式的根:
roots2=roots(p2),roots3=roots(p3) 给出:
图10-1中显示出两个多项式的图形。
图10-1 多项式p2(x)=3x2+2x,4和p3(x)=2x 3,2
(e) 多项式p(x)的牛顿-拉普森迭代,多项式的系数在向量 p中:
(f) 命令roots(poly(A))求得矩阵A的特征值。假设矩阵 A为:
运行命令:
usedRoots=roots(poly(A)) 结果为:
MATLAB 5 手册 138 下载
然而这样得到的特征值没有用 MATLAB命令eig(A)得到的特征值的精度高,而且有效 性也差些:
usedEig=eig(A)
给出:
结果是一样的,但是顺序正好相反。
(g) 对于所有的矩阵A都有:polyvalm(poly(A), A)。0 =
这是Cayley-Hamilton法则。这个法则对于秩为 5的方阵来说:
?
10.2 函数的零值
MATLAB的M文件可以表示数学函数;参见 2.9节。函数:
如果输入下面的 M文件g.m,这个函数就可以在 MATLAB中调用:
使用元素运算符.*、./、.^、+和,定义MATLAB函数g。结果是如果这个函数被一个向量调用, 那么得到的结果也是一个向量。本章中提到的所有MATLAB函数需要以这种方式来定义数学函数。
用plot命令可以画出函数的图形:
x=linspace(0, 2); % 生成向量x
plot(x,g(x)); % 画g(x)图形
grid; % 画格栅
% title('The g(x) functi'o)n 给出图标题
或者使用fplot命令:
fplot('g', [0 2]); % 画g(x)图形
grid; % 画格栅
title('The g(x) functi'o)n % 给出图标题
第10章 函数、插值和曲线拟合分析 139 下载
结果如图 10-2所示。命令plot和fplot都定义在 13.1 节中。
求函数 f(x)的零值就等于求方程 f(x) = 0的解。单变量函数的零值可以用 M AT L A B命令 fzero来求。对于多项式可以用 roots命令来求,参见 10.1节。fzero用迭代法来求解,使 得初始的估计值接近理想的零值。
g(x)函数
图10-2 用fplot 画的g(x)图形
函数的零值
fzero(fcn,x0) 求函数的一个零值, fcn为函数的名字。要求给出一 个初始值x0,近似值的相对误差为 eps。 求函数的一个零值, fcn为函数的名字。要求给出一 fzero(fcn,x0,tol) 个初始值 x0,由用户定义近似值的相对误差 tol。 求函数的一个零值,同上。如果 pic不为零,则给出迭 fzero(fcn,x0,tol,pic) 代过程。 求多变量函数的零值,fcn=fcn(x0, p1, p2,...)。如果tol和pic fzero(fcn,x0,tol,
pic p1 p2,...) 没有给出,则令它们为空矩阵,如 fzero(fcn,x0,
[],[],p1)。
zerodemo命令给出了一个演示实列。
? 例10.2
(a) 求本节开头定义的函数 g(x)的零值:
MATLAB 5 手册 140 下载
结果为:
(b) 求函数 sinx和2x,2的交集,也就是求方程 sin x=2x-2的解。先定义函数 sinm(x),将它
存放在M文件sinm.m中,如下:
画出曲线是找到初始值的一个好
方法
快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载
,所以:
结果如图 10-3所示。可以看出 2是一个可接受的估计值,输入:
xzero=fzero('sinm', 2)
结果为:
xzero =
1.4987
这就是方程 sinx=2x,2的解。
sin(x)-2*x+2函数
图10-3 sinm(x)函数
?
第10章 函数、插值和曲线拟合分析 141 下载
10.3 函数的最小值和最大值
最优化是求最优解,也就是在某个区间内有条件约束或者无条件约束地找到函数的最大 值或者最小值。 MATLAB使用数字方法求函数的最小值。使用迭代算法,也就是有些步骤要
min。 重复许多次。现在,假设要求函数 f在某个区间内的最小值 x
迭代方法需要一个初始估计值 x0。从x0开始找到一个更接近 x min的新值 x1,这个值的好坏取 决于使用的数学方法。直到找到有足够精度的近似值 xi才停止迭代,也就是绝对值 |x min,xi|足 够小。
如果函数有多个局部最小值, fmin可以找到它们中的一个。也可用 MATLAB的最优化工 具箱来求得,参见附录 C。
这里提到了
标准
excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载
MATLAB系统的两个最优化命令, fmin命令可以求单变量函数的最小 值;fmins命令可以求多变量函数的最小值,同时它还要求有一个初始向量。
没有求函数f的最大值的命令,相反函数 h=,f 的最小值可以求得。
函数的最小值
fmin(fcn,x1,x2) 求函数在区间(x1,x2)内的最小值,fcn是目标函数名。如 果没有局部最小值,则返回区间内的最小 x值。相对误差
小于10,4。 求函数在区间 (x1,x2)内的最小值,fcn是目标函数名。 fmin(fcn,x1,x2,
options) 如果没有局部最小值,则返回区间内的最小 x值。向量
options为控制参数,如 options(1)=1,显示中间结果;
options(2)表示得到的结果 x的精度,缺省为 10,4。输入 help foption可得更多信息。 s 求函数fcn的最小值。由用户自己给出一个初始估计向量 fmins(fcn,x0)
x0,相对误差为 10,4。
fmins(fcn,x0, 带优化参数求函数fcn的最小值,同上。输入help fmins 和
options) help foptions可得更多信息。例如,优化参数可以控
制迭代次数和计算结果的精度。
? 例10.3
(a) 在区间[0,2 ]内求函数 cos的最小值:
cosmin=fmin('cos', 0, pi) % 求cos的最小值 * 2
cosmin=
3.1416
这就是期望得到的结果。
(b) 同样可以简单地求高级函数的最小值。对定义在 10.2节中的函数g(x),求在区间 [0,2] 内的最小值。
MATLAB 5 手册 142 下载
注意,这是一个局部最小值,不一定是函数在这个区间内的最小值。从图 10-2上可以看
出在一个更小区间内可以得到第二个最小值,这个值比第一个值还小:
(c) 还可以用fmin命令来求函数的最大值,但是要先编写一个返回, g(x)的函数,这个函
数保存在 M文件minusg.m中。
求这个函数的最小值就等于求函数 g的最大值。
gmax=fmin('minusg', 0, 2)
结果为:
在这个区间内有若干个最大值, MATLAB求出的最大值不一定是函数的全局最大值。
(d) 用fmins命令来求多个变量函数的最小值,假设函数为:
编写 M文件fx1x2.m:
函数fmins要求有一个初始估计向量,假设给 (1,0):
结果为:
用下面的程序画出函数的图形来:
x=linspace(,1, 1 50); % 新建向量x,假设y=x
for i=1: 50 % 计算fx1x2在每一点的值
for j=1: 50
Z(i, j)=fx1x2([x(i) x(j)]); end end
meshc(x ,x, Z); % 带有基本等值线的网格图
view(80, 10); % 指定观察点
第10章 函数、插值和曲线拟合分析 143 下载
命令meshc画出函数的表面图形,同时在xy平面画出图形的等值线。命令meshc和view定 义在13.5节中,在4.2节中提到了命令linspace。结果如图10-4所示,从图上可以看出最小值。
图10-4 函数x12 +x22 - 0 . 5 x 1x2-sinx1在区间[,1,1]×[,1,1]上的图形
?
10.4 插值、曲线拟合和曲面拟合
如果在有限个数据点内给出函数,那么利用插值的方法就可以找到中间点的近似值。最 简单的插值就是对两个相邻数据点进行线性插值。 interp1和interp2命令用特殊的算法来 进行等距离数据点的快速插值。使用时,必须在方法的名字前加上一个星号, ’ * ’,如 interp1(x, Y, xx ,’*cubic’) 。
MATLAB中有几个函数可以用不同的方法来进行数据插值。
插值 interp1(x,y,xx) 返回一个长度和向量 xx相同的向量f(xx)。函数f 由向量 x和y定义,形式为 y=f(x),用线性插值 的方法来计算值。为了得到正确的结果,向量 x 必须按升序或降序排列。 返回一个相应向量的矩阵 F(xx),同上。矩阵 Y interp1(x,Y,xx) 的每列是一个关于 x的函数,对于每个这样的函
数xx的值通过插值得到。矩阵 F(xx)的行数和向 量xx的长度相同,列数和矩阵 Y的相同。 进行一维插值,字符串metstr规定不同的插值方 interp1(x,y,xx,
metstr) 法,可用的方法有:
MATLAB 5 手册 144 下载
‘linear’
线性插值。
‘nearest’
最邻近插值。
‘spline’
三次样条插值。也叫外推法。
‘cubic’ 三次插值,要求 x的值等距离。
所有插值方法均要求 x是单调的。 interp1q(x,y,xx) 和interp1相同,但是对于非均匀空间的数据插
值更快。 进行矩阵Xx和Yy的二维插值,并由X、Y和Z interp2(X,Y,Z,Xx,
Yy) 所描述的函数Z=f(X,Y)内的插值所决定。如果X、
Y和Z中任何一个是一个向量,则它的元素被认
为应用于相应的行和列。 进行二维插值,字符串metstr规定了不同的插 interp2(X,Y,Z,Xx,
Yy,metstr) 值方法,可用的方法有:
‘linear’ 线性插值。
‘nearest’ 最邻近插值。
‘spline’ 三次样条插值。
‘cubic’ 三次插值。
VV=interp3(X,Y,Z, 进行由X,Y和Z所描述的函数V的插值,XX、 V,XX,YY,ZZ, YY和ZZ是插值点。字符串metstr规定了不同 metstr) 的插值方法,可用的方法有:
‘nearest’ 使用最邻近点的值。
‘linear’ 使用8个最邻近点进行线性插值。
‘三次样条插值。 spline’
‘cubic’ 使用64个最邻近点进行三次插值。
b i c ’ VV=interpn(X1,X2, 和interp3相同,但是 V和VV可以是多维数组。 X3,...,V,Y1,,Y2, 如果X 1,X2,X3,...是等距离的,使用星号*如 Y3,...,method) ’ * c u 可以加快计算速度。
Interpft(y,n) 快速傅利叶变换插值,返回一个长度为n,从y计
算得到的向量。要求 y的值是等距离的,结果在 与y相同区间内计算。 返回相同大小的矩阵 Xx和Yy,它们表示一个网 griddata(x,y,z,
Xx,Yy,’method’) 格,由函数 z=f(x, y)的插值得到。向量 x、y和z
包含的是三维空间的 x 、 y 和 z 坐标。字符串
method规定了不同的插值方法,可用的方法如
下:
‘linear’ 基于三角的线性插值。
‘nearest’ 最邻近插值。
第10章 函数、插值和曲线拟合分析 145 下载
‘cubic’
基于三角的三次插值。
‘v4’使用MATLAB 4的插值方法。
[X1,X2,X3,...]= 变换由向量x1,x2,x3...给出的域。对于矩阵
ndgrid(x1,x2,x3, X1,X2,X3,...来说,可以用做多变量函数的
...) 估计值和多维插值。矩阵 Xn的第n维和向量 xn
的元素相同。 等同于[X1,X2,...]= ndgrid(x,x,x,...)。 [X1X2,...]
=ndgrid(x)
? 例10.4
做一个函数 sinx2在区间[0,2 ]上40个函数值的表。
(a) 用interp1来计算中间点的 sinx 2函数值,而不用sin求,命令为:
结果为:
和用sin计算得到的正确结果比较:
在表中用更多的数据点可以使精度更好。
(b) 用样条插值可得更高精度的结果。假设向量 x和y定义如上,那么:
结果为:
这是一个更好的近似值。
?
命令griddata可以在三维空间内建立任意数据点外的函数。
? 例10.5
先生成三个有 10个元素的向量,它们的值随机分布在 0,1之间:
建立一个网格来计算内部曲面之后,可以用命令 griddata来对这些点之间的曲面进行 插值。下面的图 10-5给出了不同的插值方法之间的区别:
MATLAB 5 手册 146 下载
图10-5 用griddata 对10个随机点的插值,左上图用的是 ‘linear ’,右上图用
的是‘cubic ’,右下图用的是 ‘nearest ’,右下图用的是 ‘v4’
stps=0:0.03:1; % 向量A的值在[0,1]之间 [X,Y]=meshgrid(stps); % 生成一个[0,1]×[0,1]坐标网格 Z1=griddata(x, y, z, % X, Y); 线性插值
Z2=griddata(x, y, z, % X,’Yc,ubic’); 三次插值
Z3=griddata(x, y, z, % X,’Yn,earest’); 最邻近插值 Z4=griddata(x, y, z, % X,’Yv,4’); MATLAB 4中的插值方法
% subplot(2, 2, 1); 画第1个子图 mesh(X, Y, Z1); % 画曲面网格
% hold on 保持当前图形
% plot3(x, y, z’o’); 画出数据点 ,
% hold off 释放图形
subplot(2 2 2); % 画第2个子图 mesh(X, Y ,Z2); % 画曲面网格 hold on % 保持当前图形
% plot3(x, y, z’,o’); 画出数据点
% hold off 释放图形
subplot(2, 2, 3); % 画第3个子图 mesh(X, Y, Z3); % 画曲面网格 hold on; plot3(x, y, z’o’); , % 画出数据点 hold off subplot(2, 2, 4); % 画第4个子图 mesh(X, Y, Z4); % 画曲面网格 hold on plot3(x, y, z’,o’); % 画出数据点 hold off
命令hold和subplot在13.3节进行了说明,命令 meshgrid在13.4节中,命令mesh和
plot3可以在13.5节中找到。结果如图 10-5所示。
?
第10章 函数、插值和曲线拟合分析 147 下载
用spline命令可以求得三次样条的近似值,还可以求得三次样条插值 pp形式的向量,向 量的元素是三次样条函数的系数。用 ppval命令可以求得三次样条函数的估计值。
三次样条数据插值 spline(x,y,xx) 等同于interp1(x,y,xx, ’spline’),但是参数必须是向量。
spline(x,y) 返回三次样条插值向量的 pp形式,它是函数 y=f(x)的近似值。
pp是’piecewise polynomial’的缩写,得到的向量元素包含
计算的三次样条系数。这个命令可以被 ppval函数使用。 在细胞数组X的函数Y内进行n维数据的插值,得到一个新集 YI=
splncore(X,Y,XI) 合XI。函数可以被interp2、interp3和interpn使用。
ppval(pp,xx) 计算三次样条函数。如果三次样条定义为pp=spline(x, y) ,
那么ppval(pp,xx)得到的结果和spline(x,y,xx)一样。 通过用给定点建立 pp函数来求分段多项式,其中 coeff(i, :)包 p=mkpp(points, 含第i个多项式的系数。 多项式的个数由 l=length(points),1确定,第 i个多项 coeff,d)
式的阶数为n=length(coeff(:))/l。
[points,coeff,l, 给出分段多项式相关信息,见上。
n,d]=unmkpp(p)
多项式可以使用最小二乘法来进行数据拟合 (也可见7.7节),用的命令是 polyfit。
多项式曲线拟合
polyfit(x,y,n) 找到次数为n的多项式系数,对于数据集合 {(xi, yi)}, 满足差的平方和最小。 返回同上的多项式 P和矩阵E。多项式系数在向量 p [p,E]=polyfit(x,y,n)
中,矩阵E用在polyval函数中来计算误差。
? 例10.6
下面给出了对向量x和y拟合不同阶的多项式图形,阶分别为 3,4,5。
p3=polyfit(x, y, 3); % 用向量x和y中元素拟合不同
p4=polyfit(x, y, 4); % 次数的多项式
p5=polyfit(x, y, 5); xcurve= -3.5:0.1:7.2; % 生成x值
p3curve=polyval(p3, xcurve); % 计算在这些x点的多项式值
p4curve=polyval(p4, xcurve);
p5curve=polyval(p5, xcurve);
MATLAB 5 手册 148 下载
结果如图 10-6所示。
图10-6 不同次数的多项式拟合曲线图
可以看出,次数越高的多项式精度越好。 5次的多项式经过所有的 6个点。对于这个特殊 的数据集合来说,这个多项式可称为内部插值多项式。
?
在MATLAB中,用legendre命令来计算标量或者向量的勒让德函数,它是在选定的区 间内形成完全直交集合的直交多项式系统。勒让德多项式是阶为零的勒让德函数,可以在给 定的数据集合内进行曲线拟合。贝塞尔函数是一个经典的特殊函数,它用在数学物理学中。
勒让德函数和贝塞尔函数
legendre(n,x) 返回一个矩阵,它的元素是在 x 内计算得到的 n 阶,
m=0,1, ...,n的相关勒让德函数值。 x的元素要求在区间
[,1,1]内。矩阵的第 1行对应m=0,并包含根据 x计算得
到的n阶勒让德多项式。
besselj(order,z) 计算第 1类贝塞尔函数,变量 order指定函数的阶, z用来
进行函数运算。 计算第 2类贝塞尔函数,变量 order指定函数的阶, z用来 bessely(n,x)
进行函数运算。
besselh(order,k,z) 计算向量z中元素的 Hankel函数值 (第3类贝塞尔函数 ),变
量k定义使用哪一类 Hankel函数。
besseli(order,z) 计算改良型第 1类贝塞尔函数,变量 order指定函数的阶,
z用来进行函数运算。
第10章 函数、插值和曲线拟合分析 149 下载
besselk(order,z) 计算改良型第 2类贝塞尔函数,变量 order指定函数的阶, z用来进行函数运算。 k=0或者不给出 k,计算 Airy函数 w=Ai(z); k=1,计算导 w=airy(k,z) 数Ai?(z);k=2,计算第 2类Airy函数 Bi(z ); k=3,计算导
数Bi?(z)。 返回一个错误标志数组 err。 [w,err]=airy(...)
10.5 信号分析
这里简短地介绍一些MATLAB中用作信号分析的命令。通过 help和demo命令可以得到更多
信息,也可参见‘信号过程工具箱’和它的手册。复数命令(2.4节)和卷积命令(10.1节)也要涉及。
信号分析
fft(x) 进行向量 x的离散傅立叶变换。如果 x的长度是 2的幂,则用快速 傅立叶变换, FFT。注意,变换没有规格化。 得到一个长度为n的向量。它的元素是x中前n个元素离散傅立叶变 fft(x,n) 换值。如果x有m
2),则返回一个同等大小的数组,
将A的每一维内的左右半部互换。
MATLAB 5 手册 150 下载
ifftshift(A) fftshift(A)命令的逆变换。
filter(b,a,x) 由向量 a和b所描述形成的滤波器对 x向量进行数字滤波,产生滤
波后的数据。输入help filte可得更多信息。 r
Y=filter2(h,X) 用在矩阵 h中的FIR滤波器处理 x中的数据。结果 Y由二维卷积计
算得到,包含卷记的中心部分,且与 X的大小相同。 Y=filter2(h, Y由二维卷积计算得出,其维数由参数 form规定,form是下列字 X,form) 符串之一:
‘full’ 返回二维卷积的全部,这样 Y就比X大。
‘same’ 等同于Y=filter2(h,X)。
‘valid’ 仅仅返回卷积的一部分,这部分是不带零插值边缘
计算的。这样Y就比X小。
? 例10.7
新建一个名为 ‘hat funtion’ 的函数,之后对其进行傅立叶变换。这个函数在 0,1处值为 0, 在0.5处值为1。
用linspace建立这个函数:
用下列命令画出函数的图形,如图 10-7所示:
进行傅立叶变换:
这个傅立叶变换得到复数值,但是只显示出实数部分。最后,对它进行逆傅立叶变换得 到原来函数:
图10-7 函数的傅立叶变换图
所有的画图命令都定义在第 13章中。
?