MATLAB的数值运算 (2).ppt

第二章 MATLAB的数值运算,MATLAB具有强大的数值运算能力，它是基于矩阵的运算工具。2.1 基本语法结构MATLAB语言的所有运算都是基于矩阵运算来完成的所有变量都定义为矩阵所有的运算都是关于矩阵的运算。对于传统的标量运算，在MATLAB中单独定义了点运算。2.1.1变量与赋值语句（1）变量名称：字母+数字，以字母开头，长度最大为63个字母，区分大小写（2）变量性质：全为矩阵，size()函数。,2.1 基本语法结构,A=1 2;3 4;5 6;size(A)ans=3 2 矩阵用“”作为标识符，1*1矩阵可省略“”。矩阵的行元素之间用空格或“，”分隔，各行之间用“；”分隔。如：A=1 2;3 4;5 6%或A=1,2;3,4;5,6 A=1 2 3 4 5 6,2.1 基本语法结构,（3）变量赋值 常数赋值：如上例中矩阵A的赋值，x=9;字符串赋值：f1=This is a string;表达式赋值：f2=exp(-2*x)*sin(x/5);(4)变量查询(A)变量值的查询：直接键入变量名（B）变量信息：who,whos命令 whos A Name Size Bytes Class A 3x2 48 double arrayGrand total is 6 elements using 48 bytes,2.1 基本语法结构,2.1.2 函数语句 MATLAB中除赋值语句外的其他语句a=1 2 3;b=4；5；6;%赋值语句c=a*b%乘运算x=0.9;y=sin(x);2.1.3 结构变量1.结构变量的创建 1）直接创建patient.name=hello；patient.bill=50；patient.home=jinan；,patientpatient=name:hello bill:50 home:jinanpatient(2).name=web；patient(2).bill=89；patient(2).home=uk；2）利用结构函数创建S=STRUCT(field1,VALUES1,field2,VALUES2,.)s=struct(type,big,little,color,red,x,3 4),数据的显示格式,2.1.4 变量精度MATLAB中一律使用双精度数可用format命令设置数据的显示格式format只是影响结果的显示，不影响计算与存储。format(short)：短格式（5位定点数）format long：长格式（15位定点数）format short e：短格式e方式format long e：长格式e方式 format bank：2位十进制 99.12（银行货币形式）format hex：十六进制格式,2.1.5 永久变量,matlab中预定义的一些特殊的量。i,j 虚数单位 Realmin 最小的正浮点数，pi 圆周率 Realmax 最大的浮点数，eps 浮点运算的相对精度 Inf 无穷大 NaN not a number，不定值,1/0Warning:Divide by zero.ans=Inf0/0Warning:Divide by zero.ans=NaN,2.2 矩阵运算,2.2.1 矩阵变量赋值方法1.直接赋值 a=1 1+2i;2+i exp(1)a=1.0000 1.0000+2.0000i 2.0000+1.0000i 2.7183 2.增量赋值格式：x=初值：增值：终值 x=1:0.1:1.2x=1.0000 1.1000 1.2000增量缺省时默认为1,2.2.1 矩阵变量赋值方法,3.初等矩阵赋值zeros(m,n)m*n全0矩阵ones(m,n)m*n全1矩阵eyes(m,n)m*n单位矩阵rand(m,n)m*n随机矩阵，01之间均匀分布randn(m,n)m*n随机矩阵，正态分布，期望值0rand(3,4)ans=0.9218 0.4057 0.4103 0.3529 0.7382 0.9355 0.8936 0.8132 0.1763 0.9169 0.0579 0.0099,2.2.1 矩阵变量赋值方法,例2.9 已知控制系统的3个特征根，构造系统的伴随矩阵。sysroot=-1 2 3;%3个特征根a=poly(sysroot);%得到特征方程的系数向量b=a(4),a(3),a(2);comp=zeros(2,1),eye(2);-b;comp=%控制系统的伴随矩阵为 0 1 0%|0 1|0 0 1%|:.|-6-11-6%|0.1|%|-a(n)a(n-1)-a(2)|,2.2.2 矩阵常规运算,矩阵的常规运算应符合矩阵维数的要求，其常规运算符有：+；-；*；.*；;.;/;./a；inv(a)矩阵翻转fliplr,flipud,rot90关于除法左除“”:相当于Ax=B的解，x=A-1B。右除“/”：相当于xA=B的解，x=BA-1 a=1 2;3 4 b=2 3;c=a/bc=0.6154 1.3846,2.2.3 矩阵特征运算,特征值函数 eig()奇异值函数 svd（）范数函数 norm（）秩函数 rank迹函数：矩阵所有对角线上元素的和称为矩阵的迹。trace（）条件数函数：判断矩阵的“病态”程度。cond（）矩阵的行列式运算函数 det（）例：计算矩阵的特征值与奇异值,a=1 2;3 4 eig(a),ans=-0.3723 5.3723,svd(a)ans=5.4650 0.3660,2.2.4 矩阵分解运算,1.奇异值分解:对任意矩阵A，存在酉阵U、V，使得U*S*V=A,其中S=diag(s1,s2,sp),si非负，且 称si为矩阵A的第i个奇异值 U,S,V=svd(X)其中XUSV,a=1;1;U,S,V=svd(a)U=-0.7071-0.7071-0.7071 0.7071S=1.4142 0V=-1,2.2.4 矩阵分解运算,2.LU分解 L,U=lu(A)又称三角分解，目的是分解成一个下三角阵L和一个上三角阵U的乘积，即ALU a=1 2 3;2 4 1;4 6 7;l,u=lu(a),l=0.2500 0.5000 1.0000 0.5000 1.0000 0 1.0000 0 0,u=4.0000 6.0000 7.0000 0 1.0000-2.5000 0 0 2.5000,2.2.4 矩阵分解运算,3.QR分解 Q,R=qr(A)做矩阵的正交三角形分解,将矩阵A做正交化分解，使得Q*R=A，其中Q为正交矩阵（其范数为1，指令norm(Q)=1)，R为对角化的上三角矩阵。,a=1 1 1;2-1-1;2-4 5;q,r=qr(a)q=-0.3333-0.6667-0.6667-0.6667-0.3333 0.6667-0.6667 0.6667-0.3333r=-3 3-3 0-3 3 0 0-3,2.3 基本数学函数,MATLAB的基本数学函数十分丰富，包括：三角函数：sin,sinh(双曲正弦），asin(反正弦），asinh,cos,cosh,acos,acosh，atan2(四象限反正切）指数函数：exp,log,log10,sqrt(平方根）复数函数：abs,angle,congj(共轭复数），image,real数值运算：fix（向0取整），floor(向负无穷取整），ceil(向正无穷取整），round(向最近整数圆整）rem(求余），sign(根据符号取值）矩阵函数：expm,logm,例：sin(2)ans=0.9093,log(8)ans=2.0794,2.4 点运算,a.b a，b两数组必须有相同的行和列两数组相应元素相乘。a./b 都是a的元素被b的对应元素除数组乘方(.)元素对元素的幂，例：a=1 2 3;b=4 5 6;z=a.2z=1.00 4.00 9.00z=a.bz=1.00 32.00 729.00,2.5 逻辑关系运算符,=等于eq=不等于ne 大于gt=大于等于ge&逻辑与and|逻辑或 or 逻辑非not,a=1:3;4:6;7:9;x=5;y=ones(3)*5;xa=x=axa=0 0 0 0 1 1 1 1 1,2.5 逻辑关系运算符,b=0 1 0;1 0 1;0 0 1;a=1:3;4:6;7:9a=1 2 3 4 5 6 7 8 9 ab=a&bab=0 1 0 1 0 1 0 0 1,nb=bnb=1 0 1 0 1 0 1 1 0,2.5 逻辑关系运算符,逻辑关系函数any向量的任意元素不为0则返回真all向量的所有元素不为0则返回真isempty判断空矩阵isequal判断相等矩阵isnumeric判断数值矩阵islogical 判断逻辑数组logical转换数值为逻辑型isnan判断不定数isinf判断无穷大元素isfinite判断有限大元素find寻找非零元素坐标,2.5 逻辑关系运算符,a=magic(5);a2=all(a3)a2=1 1 0 0 0a11=any(a(:,1)10)a11=1a22=any(a10)a22=1 1 1 1 1,a=17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9,matlab语言把多项式表达成一个行向量，该向量中的元素是按多项式降幂排列的。f(x)=anxn+an-1xn-1+a0 可用行向量 p=an an-1 a1+a0表示poly 产生特征多项式系数向量特征多项式一定是n+1维的特征多项式第一个元素一定是1,2.6 多相式运算,例:a=1 2 3;4 5 6;7 8 0;p=poly(a)p=1.00-6.00-72.00-27.00 p是多项式p(x)=x3-6x2-72x-27的matlab描述方法，我们可用：p1=poly2str(p,x)函数文件，显示数学多项式的形式p1=x3-6 x2-72 x-27,2.6 多相式运算,2.6 多相式运算,2.roots 求多项式的根a=1 2 3;4 5 6;7 8 0;p=poly(a)p=1.00-6.00-72.00-27.00r=roots(p)r=12.12-5.73 显然 r是矩阵a的特征值-0.39,用poly令其返回多项式形式p2=poly(r)p2=1.00-6.00-72.00-27.00matlab规定多项式系数向量用行向量表示，一组根用列向量表示。,2.6 多相式运算,3.conv,convs多项式乘运算4.deconv多项式除运算5.polyder函数多项式的微分,2.7、代数方程组求解,matlab中有两种除运算：左除和右除。对于方程ax+b，a 为anm矩阵，有三种情况：(n为方程个数，m为未知数个数）当n=m时，此方程成为“恰定”方程 当nm时，此方程成为“超定”方程 当nm时，此方程成为“欠定”方程 matlab定义的除运算可以很方便地解上述三种方程,1.恰定方程组的解,方程ax=b(a为非奇异)x=a-1 b 矩阵逆两种解:x=inv(a)*b 采用求逆运算解方程 x=ab 采用左除运算解方程,方程ax=ba=1 2;2 3;b=8;13;x=inv(a)*b x=ab x=x=2.00 2.00 3.00 3.00,=,a x=b,例:x1+2x2=8 2x1+3x2=13,2.超定方程组的解,方程 ax=b,mn时此时不存在唯一解。方程解(a a)x=a b x=(a a)-1 a b 求逆法 x=ab matlab用最小二乘法找一 个准确的基本解。,例:x1+2x2=1 2x1+3x2=2 3x1+4x2=3a=1 2;2 3;3 4;b=1;2;3;解1 x=ab x=1.00 0,=,a x=b,3.欠定方程组的解,当方程数少于未知量个数时,即不定情况,有无穷多个解存在。matlab可求出两个解：用除法求的解x是具有最多零元素的解，另一种是具有最小长度或范数的解，这个解是基于伪逆pinv求得的。,x1+2x2+3x3=1 2x1+3x2+4x3=2a=1 2 3;2 3 4;b=1;2;x=ab x=pinv(a)b x=x=1.00 0.83 0 0.33 0-0.17,a x=b,2.8 统计分析,2.8.1 基本统计分析基本统计分析函数位于.toolboxmatlabdatafunmax()找出最大值 min()找出最小值Mean()计算平均值 median()计算中间值Std()计算标准差 sort()元素排序Sun()求和函数 prod()计算元素的积Cumsum()计算元素累计和 cumprod()计算元素累计积例：已知系统的传递函数为，试确定其单位阶跃响应的最大峰值和对应时间？,n=5;d=1 1.2 5;y,x,t=step(n,d),ym=max(y)ym=1.4166 tm=spline(y,t,ym)tm=1.4435,2.8.2 相关分析,主要有两个函数：cov()和corrcoef()1.协方差矩阵函数cov()Cov(X)计算Cov(X)=E(X-EX)T(X-EX)Cov(X,Y)计算两个向量X Y的协方差2.相关矩阵系数P=corrcoef(X)计算相关矩阵P=corrcoef(X,Y)计算两个向量X Y的相关系数,插值的定义是对某些集合给定的数据点之间函数的估值方法。当不能很快地求出所需中间点的函数时，插值是一个非常有价值的工具。Matlab提供了一维、二维、三次样条等许多插值选择,2.9 函数插值运算,table1 table2 interp1 interp2 spline,插值函数,Table1和table2已基本被interp1和interp2代替,2.9 函数插值运算,1.一维插值函数interp1调用格式：YI=interp1(x,y,XI,method)X、y是已知的基准数据向量对，且x必须单调排列，XI为插值点坐标值向量，method为插值算法，默认为线性插值。,x=0:2*pi/10:2*pi;y=sin(x);xi=0:2*pi/180:2*pi;yi=interp1(x,y,xi);yii=sin(xi);stem(x,y);hold on plot(xi,yi)plot(xi,yii,r),2.9 函数插值运算,2.二维插值函数interp2调用格式：ZI=interp1(x,y,z,XI,YI,method)X、y、z是已知的基准数据，XI、YI自变量对组，method为插值算法，默认为线性插值。3.三次样条函数插值Yy=spline(x,y,xx)根据样点数据（x,y），求xx对应的三次样条插值yy。例：x=0:2*pi/360:2*pi;y=sin(x).*cos(x);xi=0:2*pi/45:2*pi;yi=spline(x,y,xi);plot(x,y)plot(xi,yi,r),2.9 函数插值运算,因求f（x）的极大值等价于求-f（x）的极小值，故MATLAB中只有处理极小值的指令X=fminbnd(fun,x1,x2)1.一元函数的极小值点2.多元函数的极小值点X=fminsearch(fun,x0)单纯形法X=fminunc(fun,x0)拟牛顿法X0是一个向量，表示极值点位置的初始猜测例1：f(x)=x2+3x+2在-5 5区间的最小值 f=fminbnd(x2+3*x+2,-5,5)f=-1.5000,2.10 函数优化（函数极值点）,function f=myfun01(x)f=x2+3*x+2;s=fminsearch(myfun01,-5)s=-1.5000f=inline(x2+3*x+2);a=fminsearch(f,-5)a=-1.5000,例2：f(x)=100(y-x2)2+(1-x)2在x=-1.2,y=1的最小值function f=xun(x)f=100*(x(2)-x(1).2).2+(1-x(1).2;X0=-1.2,1;x=fminsearch(xun,x0)x=1.0000 1.0000,2.10 函数优化（函数极值点）,2.11 数值运算,2.11.1 数值积分常用一元函数求闭区间积分运算的函数如下：Quad 采用自适应simpson法计算积分，精度较高，较常用Quadl 采用自适应lobatto法计算积分，精度高，最常用Trapz 采用梯形法求定积分，速度快，精度差Cumsum 等宽矩形法求一个区间上的积分曲线，精度差，一般不用调用格式：Quad（l）（fun，a，b，）%fun可为字符串或m函数例：求,fun=exp(-x.*x);quad(fun,0,1),ans=0.7468 quadl(fun,0,1)ans=0.7468,quad(exp(-x.*x),0,1)ans=0.7468 F=inline(exp(-x.*x);quad(F,0,1)ans=0.7468 quad(myfun01,0,1)ans=0.7468,x=1 2 3 4 5;diff(x)ans=1 1 1 1x=1 2 3 4 5;2 5 4 7 6;0 2 3 4 8;diff(x)ans=1 3 1 3 1-2-3-1-3 2,2.11 数值运算,2.11.2 微分方程数值求解（初值常微分方程）Ode23 采用低阶算法 ode45 采用中阶算法计算步骤：（1）列出微分方程和相应的初始条件（2）运用变量替换，把高阶方程化为一阶微分方程（3）定义描述微分方程的m文件（4）利用ode函数解方程。例：求解微分方程解：（1）降阶,2.11.2 微分方程数值求解,(2)定义m函数function xs=nonel(t,x)xs=zeros(2,1);xs(1)=x(2);xs(2)=-2*x(1)-x(1)2-0.5*x(2);（3）求解,t0=0;tf=20;x0=0 1;tspn=t0 tf;t,x=ode23(nonel,tspn,x0);plot(t,x),2.11 数值运算,2.11.3 非线性方程组数值求解FSOLVE(FUN,X0)使用最小二乘法求解非线性方程组，其中FUN为定义的方程组，X0为初值向量或初值矩阵。例：求解非线性方程组初值为x(0)=Y(0)=z(0)=1先定义方程组的m函数equ.m：function xs=equ(xi)x=xi(1);y=xi(2);z=xi(3);xs=zeros(3,1);xs(1)=sin(x)+y2+log(z)-7;xs(2)=3*x+2y-z3+1;xs(3)=x+y+z-5;,x0=1 1 1;xyz=fsolve(equ,x0)xyz=0.5991 2.3959 2.0050,