第6讲：MATLAB数据分析与应用

第6讲：MATLAB数据分析与应用6.1多项式运算6.2线性方程组的求解6.3非线性方程组的数值求解6.4函数的零点和极值6.5基本统计处理6.6数据插值和曲线拟合16.1多项式运算matlab语言把多项式表达成一个行向量，该向量中的元素是按多项式降幂排列的。f(x)=anxn+an-1xn-1+……+a0

可用行向量p=[anan-1……a1a0]表示。1.poly(x)——x是矩阵或向量，产生首元素是1的n+1维的特征多项式的系数向量。2例1:

a=[12.1229-5.7345-0.3884];p=poly(a)

p=1.00-6.00-72.00-27.00例2:

a=[123;456;780];p=poly(a)

p=1.00-6.00-72.00-27.00友矩阵：poly(compan(p))3p是多项式p(x)=x3-6x2-72x-27的matlab描述方法，可用：p1=poly2str(p,'x')—函数，显示多项式的数学形式：

p1=x^3-6x^2-72x-27matlab指令为

p=[1,-6,-72,-27];

p1=poly2str(p,'x')4例3:

a=[123;456;780];p=poly(a)

p=1.00-6.00-72.00-27.00r=roots(p)

r=12.12-5.73——显然r是矩阵a的特征值-0.392.roots(A)——A是多项式的系数，求多项式的根

5当然我们可用poly令其返回多项式形式p2=poly(r)

p2=1.00-6.00-72.00-27.00

matlab规定多项式系数向量用行向量表示，一组根用列向量表示。6Y=polyval(P,x)：代数多项式求值，若x为一数值，则求在该点的值；若x为向量或矩阵，则对向量或矩阵中的每个元素求其值。polyvalm(P,A)：矩阵多项式求值，A为方阵，它以方阵为自变量求多项式的值。设A为方阵，P代表多项式x3-5x2+8，则后者含义是：A^3-5*A^2+8*eye(size(A))，而前者是：A.^3-5*A.^2+8*ones(size(A))。实验：p=[1,-5,0,8];A=compan(p);polyvalm(p,A)3.求多项式的值7例4:

a(x)=x2+2x+3;b(x)=4x2+5x+6;c=(x2+2x+3)(4x2+5x+6)a=[123];b=[456];c=conv(a,b)=conv([123],[456])

c=4.0013.0028.0027.0018.00p=poly2str(c,'x')

p=4x^4+13x^3+28x^2+27x+184.conv(),convs()——多项式乘运算8a=[1,2,3];c=[4.00,13.00,28.00,27.00,18.00]d=deconv(c,a)

d=4.005.006.00[d,r]=deconv(c,a)余数c除a后的整数5.deconv()——多项式除运算9命令格式：

polyder(p):求p的微分polyder(a,b):求多项式a,b乘积的微分[p,q]=polyder(a,b):求多项式a,b商的微分例5：a=[12345];poly2str(a,'x')

ans=x^4+2x^3+3x^2+4x+5b=polyder(a)

b=4664poly2str(b,'x')

ans=4x^3+6x^2+6x+46.polyder函数——多项式微分10

1.利用左除运算符的直接解法对于方程ax=b，a为an×m矩阵，有三种情况：当n=m时，此方程成为“恰定”方程当n>m时，此方程成为“超定”方程当n<m时，此方程成为“欠定”方程matlab定义的左除运算可以很方便地解上述三种方程。6.2线性方程组的求解11方程ax=b(a为非奇异)，x=a-1b矩阵逆两种解:x=inv(a)

b—求逆运算解方程x=a\b—左除运算解方程（等价）例6:x1+2x2=82x1+3x2=13a=[12;23];b=[8;13];

x=inv(a)*b

x=a\bx=x=2.002.003.003.00

1）“恰定”方程组12方程ax=b,m<n时，不存在唯一解。方程解(a'a)x=a'b,x=(a'a)-1a'b——求逆法x=a\b——最小二乘法

(等价)例7:

x1+2x2=12x1+3x2=23x1+4x2=3，a=[12;23;34];b=[1;2;3];

x=a\b

x=inv(a'

a)

a'

b

x=x=1.001.0000.00

2）“超定”方程组13方程ax=b,m>n时，即不定情况，有无穷多个解存在。matlab可求出两个解：1）具有最多零元素的解——除法求解2）具有最小长度或范数的解——伪逆pinv例8:x1+2x2+3x3=12x1+3x2+4x3=2，a=[123;234];b=[1;2];

x=a\b

x=pinv(a)

b

x=1.00x=0.8300.330-0.173）“欠定”方程组14

矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。

(1)LU分解：将非奇异的矩阵分解为一个下三角矩阵L和一个上三角矩阵U的乘积，即A=L*U。Ax=b可变换为x=U\(L\b)或x=U\(L\Pb)。

分解格式(矩阵A必须是方阵)：[L,U]=lu(A)：产生一个上三角阵U和一个下三角阵L。[L,U,P]=lu(A)：产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P2.利用矩阵分解的直接解法15命令如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';[L,U]=lu(A);x=U\(L\b)或采用LU分解的第2种格式，命令如下：[L,U,P]=lu(A);x=U\(L\P*b)例9

用LU分解求解例7-1中的线性方程组。16

(2)QR分解：把矩阵A分解为一个正交矩阵Q和一个上三角矩阵R的乘积，即A=Q*R。Ax=b可变换为x=R\(Q\b)或x=E(R\(Q\b))。

分解格式：[Q,R]=qr(A)：产生一个正交矩阵Q和一个上三角矩阵R。[Q,R,E]=qr(A)：产生一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E。QR分解17(3)Cholesky分解：把对称正定的矩阵A分解成一个上三角矩阵R和其转置矩阵R'的乘积，即A=R'R。Ax=b可变成x=R\(R'\b)。分解格式：R=chol(A)：产生一个上三角阵R。[R,p]=chol(A)：A对称正定时，p=0；否则p为一个正整数。如果A满秩，则R为一个满足R‘R=A(1:p-1,1:p-1)的阶为p-1的上三角阵。Cholesky分解18@函数Matlab6.0以后的版本中可以用@函数进行函数调用，其调用格式：handle=@function这类似于C++语言中的引用。实例：利用句柄传递数据。fhandle=@humps;%为hump函数创建一个函数句柄，并将它指定为fhandle变量x=fminbnd(@humps,0,1)%传递humps句柄给fminbnd函数，传递所有变量匿名函数19对于非线性方程组F(X)=0，用fsolve函数求其数值解。fsolve函数的调用格式为：X=fsolve('fun',X0)

X=fsolve('fun',X0,options)其中X为返回的解，fun是用于定义需求解的非线性方程组的函数文件名，X0是初值，options为选项设定，默认缺省。如果想改变其中某个选项，则可以调用optimset()函数来完成。6.3非线性方程组的求解20例10：求解方程在[-1,10]内的解。

首先作图如下：

fplot(‘[5*x^2*sin(x)-exp(-x),0]’,[0,10])发现在[0,10]区间中有4个解，分别在0,3,6,9附近，所以用命令：

fun=inline(‘5*x.^2.*sin(x)-exp(-x)’);fsolve(fun,[0369],1e-6)

ans=0.50183.14076.28329.4248注意：用inline建立的函数fun调用时不加单引号；而用M函数文件建立的函数则要加单引号或用@函数。21

首先编制函数文件fc.m如下：functiony=fc(x)y=[x(1)-0.7*sin(x(1))-0.2*cos(x(2));x(2)-0.7*cos(x(1))+0.2*sin(x(2))];

在MATLAB命令窗口输入：x0=[0.50.5];fsolve(‘fc’,x0)或fsolve(@fc,x0)ans=0.52650.5079例11：求非线性方程组在(0.5,0.5)附近的解。22对于函数f(x)的零点，用fzero函数求解。基本格式为：x=fzero('fun',X0)[x,fval]=fzero(fun,X0,options)求函数零点x和该点的函数值fval，fun可以是字符串、内联函数或M函数文件名的函数句柄。X0是初值，options为选项设定，默认缺省。6.4函数的零点23例12：求的零点。首先作图确定零点范围：x=-5:0.1:5;f=x.^3-2*x-5;plot(x,f,x,0);发现零点大致在2附近，编写M文件，函数名zero1.m

functionf=zeros(x)f=x.^3-2*x-5;计算2附近的零点和该点的函数值[xz,yz]=fzero(@zero1,2)

xz=yz=2.0946-8.8818e-01624对于函数f(x)的极值，用fmidbnd函数求解极小值。基本格式为：

x=fmidbnd('fun',x1,x2)[x,fval]=fmidbnd(fun,x1,x2,options)求函数fun在(x1,x2)间的极小点x及该点函数值fval，用法同fzero。6.5函数的极值25例13：求的极小值。首先作图确定极小值范围：x=-5:0.1:5;f=(x-3).^2-1;plot(x,f);发现在(2,4)间有极小值，编M文件，函数名为m1.m

functionf=m1(x)

f=(x-3).^2-1;计算(2,4)间的极小点和该点的函数值[xm,ym]=fminbnd(‘m1’,2,4)xm=ym=3-1261、查取最大值——max函数[Y,I]=max(X)：返回矩阵X的各列中的最大元素值及其该元素的位置并赋予行向量Y与I；当X为向量时，则Y与I为单变量。[Y,I]=max(X,[],DIM)：按数组X的第DIM维的方向查取其最大的元素值及其该元素的位置赋予向量Y与I。6.6基本统计处理27U=max(A,B)：A,B是两个同型的向量或矩阵，结果U是与A,B同型的向量或矩阵，U的每个元素等于A,B对应元素的较大者。U=max(A,n)：n是一个标量，结果U是与A同型的向量或矩阵，U的每个元素等于A对应元素和n中的较大者。28例14

查找下面数列x的最大值。

x=[359618]%产生数列xx=359618

y=max(x)%查出数列x中的最大值赋予yy=9

[y,l]=max(x)%查出数列x中的最大值及其该元素的位置赋予y,ly=9l=329例15

分别查找下面3×4的二维数组x中各列和各行元素中的最大值。

x=[1842;9625;3671]%产生二维数组xx=184296253671

y=max(x)%查出二维数组x中各列元素的最大值产生赋予行向量yy=987530

[y,l]=max(x)%查出二维数组x中各列元素的最大值，元素的行下标赋予y,ly=9875l=2132[y,l]=max(x,[],1)%本命令的结果同上

y=9875l=2132

[y,l]=max(x,[],2)%DIM=2，在各行中查找y=8l=29173

31例16

试取下面两个2×3的二维数组x、y所有同一位置上的元素值大者构成一个新矩阵p。

x=[456;148];%产生二维数组x

y=[175;457];%产生二维数组y

p=max(x,y)%在x,y同一位置上的两个元素中查找出最大值，赋予与x,y同样大小的二维数组pp=476458322、查取最小值MIN函数用来查取数据序列的最小值。3、求中值Y=median(X)：将median(X)返回矩阵X各列元素的中值；若X为向量，则Y为单变量。Y=median(X,DIM)：按数组X的第DIM维方向的元素求其中值赋予向量Y。若DIM=1，为按列操作；若DIM=2，为按行操作。若X为二维数组，则Y为一个向量；若X为一维数组，则Y为单变量。33例17

对二维数组x试从不同维方向求出其中值。x=[1842;9625;3671];y0=median(x)%按列操作y0=3642y1=median(x,1)%按列操作,结果为行向量y1=3642y2=median(x,2)%按行操作,结果为列向量y2=3.00005.50004.5000344、求和Y=sum(X)：将sum(X)返回矩阵X各列元素之和赋予行向量Y；若X为向量，则Y为单变量。Y=sum(X,DIM)：按数组X的第DIM维的方向的元素求其和赋予Y。若DIM=1，为按列操作；若DIM=2，为按行操作。若X为二维数组，Y为一个向量；若X为一维数组，则Y为单变量。355、求平均值Y=mean(X)：将mean(X)返回矩阵X各列元素之和的平均值赋予行向量Y。若X为向量，则Y为单变量。Y=mean(X,DIM)：按数组X的第DIM维的方向的元素求其平均值赋予向量Y。若DIM=1，为按列操作；若DIM=2，为按行操作。若X为二维数组，Y为一个向量；若X为一维数组，则Y为单变量。36例18x=[456;148];y=sum(x,1)y=5914y=sum(x,2)y=1513y1=mean(x,1)y1=2.504.507.00y2=mean(x,2)y2=5.00004.3333376、求积Y=prod(X)：将prod(X)返回矩阵X各列元素之积赋予行向量Y。若X为向量，则Y为单变量。Y=prod(X,DIM)：按数组X的第DIM维的方向的元素求其积赋予向量Y。若DIM=1，为按列操作；若DIM=2，为按行操作。若X为二维数组，Y为一个向量；若X为一维数组，则Y为单变量。38例19x=[456;148];y1=prod(x,1)y1=42048y2=prod(x,2)

y2=12032397、累加和、累乘积cumsum(X)：求向量或矩阵X的累加和。cumprod(X)：求向量或矩阵X的累加积。cumsum(A,dim)：当dim为1时，该函数等同于cumsum(A)；当dim为2时，其第i行是A的第i行的累加和向量。cumprod(A,dim)：同上。408、标准方差std(X)：求向量或矩阵X的标准方差。Y=std(A,flag,dim)：当dim=1时，求各列元素的标准方差；当dim=2时，求各行元素的标准方差。flag取0或1，当flag=0时，按σ=1计算；当flag=1时，按σ=2计算。缺省flag=0，dim=1。419、相关系数

corrcoef(X)：返回从矩阵X形成的一个相关系数矩阵。此相关系数矩阵的大小与矩阵X一样。它把矩阵X的每列作为一个变量，然后求它们的相关系数。

corrcoef(X,Y)：在这里，X,Y是向量，它们与corrcoef([X,Y])的作用一样。

42例20

生成满足正态分布的10000×5随机矩阵，然后求各列元素的均值和标准方差，再求这5列随机数据的相关系数矩阵。命令如下：X=randn(10000,5);M=mean(X)D=std(X)R=corrcoef(X)4310、升序排序sort(X)：返回对向量X中的元素按升序排列的新向量。[Y,I]=sort(A,dim)：对矩阵A的各列或各行重新排序，其中dim指明对A的列还是行进行排序。若dim=1，则按列排；若dim=2，则按行排。Y是排序后的矩阵，而I记录Y中的元素在A中位置。441、一维数据插值Y1=interp1(X,Y,X1,‘method’)：根据X,Y的值计算函数在X1处的值。X,Y是两个等长的已知向量，分别描述采样点和样本值，X1是一个向量或标量，描述欲插值的点，Y1是一个与X1等长的插值结果。method是插值方法，允许的取值有linear(线性插值)、nearest(最近邻点)、cubic(三次函数)、spline(三次样条)。6.4数据插值和曲线拟合(cftool)45例21

某观测站测得某日6:00时至18:00时之间每隔2小时的室内外温度(℃)，用3次样条插值分别求得该日室内外6:30至17:30时之间每隔2小时各点的近似温度(℃)。h=6:2:18;%时间变量t=[18,20,22,25,30,28,24;15,19,24,28,34,32,30]’;%已知数据XI=6.5:2:17.5%欲插值点YI=interp1(h,t,XI,‘spline’)%3次样条插值46Z1=interp2(X,Y,Z,X1,Y1,‘method’)：向量X,Y分别是描述两个参数的采样点，Z是与采样点对应的函数值，X1,Y1是两个向量或标量，描述欲插值的点。Z1是根据相应的插值方法得到的插值结果。method的取值与一维插值函数相同。X,Y,Z也可以是矩阵形式。同样，X1,Y1的取值范围不能超出X,Y的给定范围，否则，会给出“NaN”错误。2、二维数据插值47例22

某实验对一根长10米的钢轨进行热源的温度传播测试。x表示测量点0:2.5:10(米)，h表示测量时间0:30:60(秒)，T表示测试所得各点的温度(℃)。试用线性插值求出在一分钟内每隔20秒、钢轨每隔1米处的温度TI。x=0:2.5:10;h=[0:30:60]';T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];xi=[0:10];hi=[0:20:60]';TI=interp2(x,h,T,xi,hi)483、曲线拟合[P,S]=polyfit(X,Y,m)：根据采样点X和采样点函数值Y，产生一个m次多项式P及其在采样点的误差向量S。其中X,Y是等长的向量，P是一个长