ATLAB矩阵运算及多项式处理
时间:2015-03-28 15:17 来源:自动控制网
一、矩阵的输入 (一)在命令窗口中输入 》a=1;b=2;c=3; 》x=[5bc;a*ba+cc/b] x= 5.0002.0003.000 2.0004.0001.500 》y=[2,4,5 368] y= 245 368 矩阵生成不但可以使用纯数字(含复数),也可以使用变量(或者说采用一个表达式)。矩阵的元素直接排列在方括号内,行与行之间用分号隔开,每行内的元素使用空格或逗号隔开。大的矩阵可以用分行输入,回车键代表分号。 (二)语句生成 1.用线性等间距生成向量矩阵(start:step:end) 》a=[1:2:10] a= 13579 其中start为起始值,step为步长,end为终止值。当步长为1时可省略step参数;另外step也可以取负数。 2.a=linspace(n1,n2,n) 在线性空间上,行矢量的值从n1到n2,数据个数为n,缺省n为100。 》a=linspace(1,10,10) a= 12345678910 3.a=logspace(n1,n2,n) 在对数空间上,行矢量的值从10n1到10n2,数据个数为n,缺省n为50。这个指令为建立对数频域轴坐标提供了方便。 》a=logspace(1,3,3) a= 101001000 4.一些常用的特殊矩阵 单位矩阵:eye(m,n);eye(m) 零矩阵:zeros(m,n);zeros(m) 一矩阵:ones(m,n);ones(m) 对角矩阵:对角元素向量V=[a1,a2,…,an]A=diag(V) 随机矩阵:rand(m,n)产生一个m×n的均匀分别的随机矩阵 》eye(2,3) ans= 100 010 》zeros(2,3) ans= 000 000 》ones(2,3) ans= 111 111 》V=[572];A=diag(V) A= 500 070 002 》eye(2) ans= 10 01 》zeros(2) ans= 00 00 》ones(2) ans= 11 11 如果已知A为方阵,则V=diag(A)可以提取A的对角元素构成向量V。
二、矩阵的运算 (一)转置:对于实矩阵用(’)符号或(.’)求转置结果是一样的;然而对于含复数的矩阵,则(’)将同时对复数进行共轭处理,而(.’)则只是将其排列形式进行转置。 》a=[123;456]' a= 14 25 36
》a=[123;456].' a= 14 25 36
》b=[1+2i2-7i]' b= 1.0000-2.0000i 2.0000+7.0000i 》b=[1+2i2-7i].' b= 1.0000+2.0000i 2.0000-7.0000i
(二)四则运算与幂运算 +;-;*;\和/;^;.*;.\;./;.^ 如:a=[12;34];b=[35;59] 》c=a+bd=a-b 》c=d= 47-2-3 813-2-5 》a*b=[1323;2951] 》a/b=[-0.500.50;3.50–1.50] 》a\b=[-1-1;23] 》a^3=[3754;81118] 》a.*b=[310;1536] 》a./b=[0.330.40;0.600.44] 》a.\b=[3.002.50;1.672.25] 》a.^3=[18;2764] 只有维数相同的矩阵才能进行加减运算。 注意只有当两个矩阵中前一个矩阵的列数和后一个矩阵的行数相同时,才可以进行乘法运算。a\b运算等效于求a*x=b的解;而a/b等效于求x*b=a的解。只有方阵才可以求幂。 点运算是两个维数相同矩阵对应元素之间的运算,在有的教材中也定义为数组运算。
(三)逆矩阵与行列式计算 求逆:inv(A); 求行列式:det(A) 要求矩阵必须为方阵 》a=[123;456;235]; 》b=inv(a) b= -2.33330.33331.0000 2.66670.3333-2.0000 -0.6667-0.33331.0000 》det(a) ans= -3 (四)了解矩阵超越函数 在MATLAB中expsqrt等命令也可以作用到矩阵上,但这种运算是定义在矩阵的单个元素上的,即分别对矩阵的每一个元素进行计算。
超越数学函数可以在函数后加上m而成为矩阵的超越函数,例如:expm,sqrtm。矩阵的超越函数要求运算矩阵为方阵。
三、矩阵的操作 (一)矩阵下标 MATLAB通过确认矩阵下标,可以对矩阵进行插入子块,提取子块和重排子块的操作。 A(m,n):提取第m行,第n列元素 A(:,n):提取第n列元素 A(m,:):提取第m行元素 A(m1:m2,n1:n2):提取第m1行到第m2行和第n1列到第n2列的所有元素(提取子块)。 A(:):得到一个长列矢量,该矢量的元素按矩阵的列进行排列。 矩阵扩展:如果在原矩阵中一个不存在的地址位置上设定一个数(赋值),则该矩阵会自动扩展行列数,并在该位置上添加这个数,而且在其他没有指定的位置补零。 消除子块:如果将矩阵的子块赋值为空矩阵[],则相当于消除了相应的矩阵子块。 (二)矩阵的大小 [m,n]=size(A,x):返回矩阵的行列数m与n,当x=1,则只返回行数m,当x=2,则只返回列数n。 length(A)=max(size(A)):返回行数或列数的最大值。 rank(A):求矩阵的秩 》a=[123;345]; 》[m,n]=size(a) m= 2 n= 3
》length(a) ans= 3 》max(size(a)) ans= 3
》rank(a) ans= 2 (三)了解矩阵操作函数:flipud;fliplr;rot90
四、多项式处理 1.多项式的建立与表示方法 在MATLAB中,多项式使用降幂系数的行向量表示,如:多项式
表示为:p=[1-12025116],使用函数roots可以求出多项式等于0的根,根用列向量表示。若已知多项式等于0的根,函数poly可以求出相应多项式。 p=poly(r) p= 1-12-025116
r=roots(p) r= 11.7473 2.7028 -1.2251+1.4672i -1.2251-1.4672i 2.多项式的运算 相乘conv a=[123];b=[12]c=conv(a,b)=1476 conv指令可以嵌套使用,如conv(conv(a,b),c) 相除deconv [q,r]=deconv(c,b) q=123%商多项式 r=000%余多项式 求多项式的微分多项式polyder polyder(a)=22 求多项式函数值polyval(p,n):将值n代入多项式求解。polyval(a,2)=11 3.多项式的拟合 多项式拟合又称为曲线拟合,其目的就是在众多的样本点中进行拟合,找出满足样本点分布的多项式。这在分析实验数据,将实验数据做解析描述时非常有用。 命令格式:p=polyfit(x,y,n),其中x和y为样本点向量,n为所求多项式的阶数,p为求出的多项式。
4.多项式插值 多项式插值是指根据给定的有限个样本点,产生另外的估计点以达到数据更为平滑的效果。该技巧在信号处理与图像处理上应用广泛。 所用指令有一维的interp(一)二维的interp(二)三维的interp3。这些指令分别有不同的方法(method),设计者可以根据需要选择适当的方法,以满足系统属性的要求。Helppolyfun可以得到更详细的内容。 y=interp1(xs,ys,x,’method’) 在有限样本点向量xs与ys中,插值产生向量x和y,所用方法定义在method中,有4种选择: nearest:执行速度最快,输出结果为直角转折 linear:默认值,在样本点上斜率变化很大 spline:最花时间,但输出结果也最平滑 cubic:最占内存,输出结果与spline差不多
五、MATLAB数据处理 (一)矩阵分解 1.奇异值分解 求矩阵A的奇异值及分解矩阵,满足U*S*V’=A,其中UV矩阵为正交矩阵(U*U’=I),S矩阵为对角矩阵,它的对角元素即A矩阵的奇异值。 [U,S,V]=svd(A) 例:a= 98 68 可以验证: u*u’=I v*v’=I u*s*v’=a
[u,s,v]=svd(a) u= 0.7705-0.6375 0.63750.7705 s= 15.57650 01.5408 v= 0.6907-0.7231 0.72310.6907 2.特征值分解 [V,D]=eig(A) 求矩阵A的特征向量V及特征值D,满足A*V=V*D。其中D的对角线元素为特征值,V的列为对应的特征向量。如果D=eig(A)则只返回特征值。 例:a= 98 68 [v,d]=eig(a) v= 0.7787-0.7320 0.62740.6813 d= 15.44620 01.5538
可以验证:A*V=V*D
3.正交分解 [Q,R]=qr(A) 将矩阵A做正交化分解,使得Q*R=A,其中Q为正交矩阵(其范数为1,指令norm(Q)=1),R为对角化的上三角矩阵。 例:a= 98 68 [q,r]=qr(a) q= -0.8321-0.5547 -0.55470.8321 r= -10.8167-11.0940 02.2188 norm(q) ans= 1
q*r ans= 9.00008.0000 6.00008.0000 4.三角分解 [L,U]=lu(A) 将A做对角线分解,使得A=L*U,其中L为下三角矩阵,U为上三角矩阵。 注意:L实际上是一个“心理上”的下三角矩阵,它事实上是一个置换矩阵P的逆矩阵与一个真正下三角矩阵L1(其对角线元素为1)的乘积。 [L1,U1,P]=lu(A) 例:a=[123;456;789] 比较:[l1,u1,p]=lu(a)[l,u]=lu(a)
l1= 1.0000 0.141.000 0.570.501.00 u1= 7.008.009.00 00.861.71 000.00 p= 001 100 010 l= 0.141.000 0.570.501.00 1.0000 u= 7.008.009.00 00.861.71 000.00 可以验证: u1=u,inv(p)*l1=l a=l*u p*a=l1*u1
(二)数据分析 1.绘制函数图形:fplot() 2.求极值:fmin,fmins 3.求零点:寻找一维函数的过零点fzero() 4.频谱分析(fft):y=FFT(x);unwrap();abs;angle画出幅频和相频曲线 5.了解数据分析函数:max,min,mean,sum,prod等 6.了解积分运算:trap2,quad,quad8 (三)常微分方程数值解 [t,x]=ode23(‘xfun’,t0,tf,x0,tol) [t,x]=ode45(‘xfun’,t0,tf,x0,tol) |