第二讲MATLAB的数值计算

1、MATLAB程序设计基础MATLAB的数值计算Matlab的数据类型o变量 变量不需要事先声明，也不需要指定变量类型，它会自动根据所赋予变量的值或对变量的操作来确定变量的类型；赋值过程中，如果变量已存在，则用新值代替旧值，以新的类型代替旧的类型。 变量的命名规则：l变量名区分大小写；l变量名长度不超过31位，第31位之后的字符被忽略；l变量名以字母开头，变量名中可以包含字母、数字、下划线，但不能使用标点。 变量一般为局部变量，即仅在其调用的M文件内部有效；若要定义全局变量，须在变量前加关键字global。o常量 matlab中预定义的一些特殊的量。 i,j 虚数单位Realmin 最小的正浮点

2、数， pi 圆周率Realmax 最大的浮点数， eps 浮点运算的相对精度Inf 无穷大 NaN not a number ，不定值 例如：pians = 3.14161102221023252101/0Warning: Divide by zero.ans = Inf0/0Warning: Divide by zero.ans = NaNo定义变量时应避免与常量名相同，如果改变了某个常量的值，可以用clear命令来恢复。pi=1pi = 1clear pipians = 3.1416o数字变量l数字变量的运算258*369ans = 95202x=258*369x = 952021233a

3、ns = 1860867sqrt(ans)ans = 1.3641e+003l数字的输入输出格式 缺省为实数保留小数点后4位浮点数表示。 其输入格式与C语言一致： 如：9 -73 0.1999 1.475e6 输出格式由format命令控制，只是影响屏幕显示效果，不影响内部存储和计算。format long;pians = 3.14159265358979format long e;pians = 3.141592653589793e+000format long g;pians = 3.14159265358979o字符串l1、字符串的约定字符串用单引号输入或赋值；字符串的每个字符都是都是字

4、符数组的一个元素；字符串和字符数组基本上等价。s=symbolics =symbolic size(s)ans = 1 8 s(3)ans =ml字符串的转换 double 字符串转换为数值代码 num2str 数字转换为字符串 int2str 整数转换为字符串 mat2str 矩阵转换为字符串 str2num 转换字符串为数字 double(s)ans = 49 50 51 50 51 52l字符串操作 strcat strcmp strvcat strncmp findstr upper lower blanks deblank l执行字符串 t=1/(a*b-1);a=2;b=3;c=e

5、val(t)c = 0.2000o结构型变量 由函数struct定义，以指针操作符“.”连接结构型变量名与属性名。 结构型变量名结构型变量名struct（元素名（元素名1，元素值，元素值1，元素名，元素名2，元素值，元素值2，） c=struct(c1,1,c2,1 2 3 4,c3,abcd)c = c1: 1 c2: 1 2 3 4 c3: abcd c.c2ans = 1 2 3 4 c.c3ans =abcdabcdo单元型变量 单元型变量为任意类型的多维数组，其定义需用大括号，元素间用逗号隔开。 a=1,2;3,4a = 1 2 3 4 b=1:4,a,abcdb = 1x4 dou

6、ble 2x2 double abcd?cellplot(b) 单元型变量元素的引用采用大括号为下标标识，用小括号只显示该元素的压缩形式。 b2ans = 1 2 3 4 b(2)ans = 2x2 doubleo向量 向量元素用“”括起来，元素间用空格、逗号或分号分隔； 注意：空格和逗号分隔成行向量，分号分割成列向量。注意：空格和逗号分隔成行向量，分号分割成列向量。l冒号表达式生成向量 基本格式：xx1:step:x2 xx1:x2 a=1:2:12a = 1 3 5 7 9 11 a=12:-2:1a = 12 10 8 6 4 2 a=1:6a = 1 2 3 4 5 6l线性等分向量生

7、成 y=linspace(x1,x2) 生成100维行向量 y=linspace(x1,x2,n) 生成n维行向量 a=linspace(1,100,6)a = 1.0000 20.8000 40.6000 60.4000 80.2000 100.0000l对数等分向量生成 y=logspace(x1,x2)生成50维对数等分向量，y(1)=10 x1 y(50)=10 x2 y=logspace(x1,x2,n)生成n维对数等分向量y(1)=10 x1 y(n)=10 x2 a=logspace(0,5,6)a = 1 10 100 1000 10000 100000o向量的基本运算l与数运

8、算a = 1.0000 20.8000 40.6000 60.4000 80.2000 100.0000 a-1ans = 0 19.8000 39.6000 59.4000 79.2000 99.0000 a*2ans = 2.0000 41.6000 81.2000 120.8000 160.4000 200.0000l点积计算 指两个向量在其中一个向量方向上的投影的乘积。 dot(a,b) a,b必须同维。 a=1 2 3; b=3,4,5; dot(a,b)ans = 26 sum(a.*b)ans = 26l叉积 表示过两相交向量的交点的垂直于两向量所在平面的向量。 cross（a，

9、b） a，b必须为三维向量。l混合积 c=cross(a,b)c = -2 4 -2 dot(a,cross(b,c)ans = 24o矩阵大型矩阵通借助M文件来输入。矩阵987654321A的输入： A = 1,2,3; 4,5,6; 7,8,9A = 1 2 3 4 5 6 7 8 9 a=1 2 34 5 67 8 9a = 1 2 3 4 5 6 7 8 9x=rand(1,5)%产生的均布随机数组产生的均布随机数组 x =0.9501 0.2311 0.6068 0.4860 0.8913 x(3)%寻访数组寻访数组x的第三个元素。的第三个元素。 ans =0.6068 x(1 2

10、5)%寻访数组寻访数组x的第一、二、五个元素组成的子数组。的第一、二、五个元素组成的子数组。 ans =0.9501 0.2311 0.8913 x(1:3)%寻访前三个元素组成的子数组寻访前三个元素组成的子数组 ans = 0.9501 0.2311 0.6068 x(3:end)%寻访除前寻访除前2个元素外的全部其他元素。个元素外的全部其他元素。end是最后一是最后一个元素的下标。个元素的下标。 ans = 0.6068 0.4860 0.8913 常用的特殊矩阵单位矩阵：eye(m,n); eye(m)零 矩 阵：zeros(m,n); zeros(m)一 矩 阵：ones(m,n);

11、ones(m)对角矩阵：对角元素向量 V=a1,a2,an A=diag(V)随机矩阵：rand(m,n)产生一个mn的均匀分别的随机矩阵 eye(2,3)ans= 1 0 0 0 1 0 zeros(2,3)ans= 0 0 0 0 0 0 ones(2,3)ans= 1 1 1 1 1 1 V=5 7 2; A=diag(V)A= 5 0 0 0 7 0 0 0 2 eye(2)ans= 1 0 0 1 zeros(2)ans= 0 0 0 0 ones(2)ans= 1 1 1 1 如果已知A为方阵，则V=diag(A)可以提取A的对角元素构成向量V。其他特殊矩阵 compan 友矩阵函

12、数 magic 魔方矩阵 hankel Hankel矩阵 rosser 对称特征值测试矩阵 hilb Hilbert矩阵 pascal Pascal矩阵 invhilb 反Hilbert矩阵 vander 范德蒙矩阵 o矩阵的基本运算l加减运算 要求两矩阵必须同阶。 a=1 2 3;2 3 4; 3 4 5; b=1 1 1;2 2 2;3 3 3; c=a+bc = 2 3 4 4 5 6 6 7 8l乘法 要求a为ij阶，b为jk阶时，ab才能相乘。 e=b,5 5 5e = 1 1 1 5 2 2 2 5 3 3 3 5 f=a*ef = 14 14 14 30 20 20 20 45

13、26 26 26 60l除法左除“”: 相当于Ax=B的解，x=A-1B。右除“/”：相当于xA=B的解，x=BA-1 A-1B=(BA-1)。通常，右除稍快一些，而左除可以避免奇异性。对于AxB，其中A为（nm）阶矩阵：n=m且非奇异时，方程为恰定方程；nm方程为超定方程；n A=1 2 3;4 5 6;7 8 0;1 3 5; B=1 3 5;2 4 6; A/Bans = 0 0.5000 -3.0000 3.5000 -12.0000 10.2500 1.0000 0.0000(BA)ans = 0 0.5000 -3.0000 3.5000 -12.0000 10.2500 1.00

14、00 0.0000 rand(seed,12); a=rand(100)+1.e8; x=ones(100,1); b=a*x; cond(a)ans = 5.0482e+011 tic;x1=b/a;t1=toct1 = 0.0773 er1=norm(x-x1)er1 = 139.8328 tic;x1=ab;t1=toct1 = 0.0014 er2=norm(x-x1)er2 = 1.1311e-004l矩阵与常数的运算 常数与此矩阵的各元素之间进行运算。 注意：进行数除时，常数通常只能做除数。l矩阵的逆运算 函数 inv A=2 1 -3 -1;3 1 0 7;-1 2 4 -2;1

15、 0 -1 5; inv(A)ans = -0.0471 0.5882 -0.2706 -0.9412 0.3882 -0.3529 0.4824 0.7647 -0.2235 0.2941 -0.0353 -0.4706 -0.0353 -0.0588 0.0471 0.2941l矩阵的行列式运算 函数 det A=2 1 -3 -1;3 1 0 7;-1 2 4 -2;1 0 -1 5; a1=det(A)a1 = -85 a2=det(inv(A)a2 = -0.0118 a1*a2ans = 1l矩阵的幂运算 与数字的幂运算形式相同，用“”算符。l矩阵的指数运算 常用函数 expm e

16、xpm1 expm2 expm3l矩阵的对数运算 函数 logml矩阵的开方运算 函数 sqrtm b=magic(3)b = 8 1 6 3 5 7 4 9 2 sqrtm(b)ans = 2.7065 + 0.0601i 0.0185 + 0.5347i 1.1480 - 0.5948i 0.4703 + 0.0829i 2.0288 + 0.7378i 1.3739 - 0.8207i 0.6962 - 0.1430i 1.8257 - 1.2725i 1.3511 + 1.4155i b0.5ans = 2.7065 + 0.0601i 0.0185 + 0.5347i 1.1480

17、- 0.5948i 0.4703 + 0.0829i 2.0288 + 0.7378i 1.3739 - 0.8207i 0.6962 - 0.1430i 1.8257 - 1.2725i 1.3511 + 1.4155io矩阵的基本函数运算l特征值函数 函数 x,y=eig(A) 可以给出特征值和特征向量的值 x为特征向量矩阵，y为特征值矩阵。 A=7 3 -2;3 4 -1;-2 -1 3;x,y=eig(A)x = 0.5774 0.0988 -0.8105 -0.5774 -0.6525 -0.4908 0.5774 -0.7513 0.3197y = 2.0000 0 0 0 2.3

18、944 0 0 0 9.6056l奇异值函数 函数 svd svds 矩阵非零特征值的平方根称作矩阵的奇异值。 矩阵的奇异值大小通常决定矩阵的性态，如果矩阵的奇异值变化特别大，则矩阵中某个参数的微小变化将严重影响到原矩阵的参数，这样的矩阵又称为病态矩阵，也称奇异矩阵。例：A=1,1;5*eps,0;0,5*eps; sqrt(eig(A*A)l矩阵翻转 函数 fliplr flipud rot90a = 7 3 -2 3 4 -1 -2 -1 3 fliplr(a)ans = -2 3 7 -1 4 3 3 -1 -2 flipud(a)ans = -2 -1 3 3 4 -1 7 3 -2

19、rot90(a)ans = -2 -1 3 3 4 -1 7 3 -2l条件数函数 判断矩阵的“病态”程度。 函数cond 计算矩阵的条件数的值 condest计算矩阵的1范数条件数的估计值 rcond 计算矩阵的条件数的倒数值 矩阵的条件数condA=最大奇异值/最小奇异值，条件数越大，则对参数变化越敏感。l范数函数 函数 norm（X，P） P1 1范数 P2 2范数 Pinf 无穷范数 Pfro F范数 norm（X）norm（X，2）njijniTniijnjaAAAsAaA11max2111|max|maxl秩函数 函数 ranke = 1 1 1 5 2 2 2 2 3 3 3 5

20、 rank(e)ans = 2l迹函数 矩阵所有对角线上元素的和称为矩阵的迹。 函数 tracel正交空间函数 函数 orth 用来求矩阵的一组正交基。l伪逆函数 函数 pinv 求解“病态”问题时，避免产生伪解。 a=magic(4)a = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 b=a*1 1 1 1; inv(a)*bWarning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.567374e-017.ans = 0 8 0 0 p

21、inv(a)*bans = 1.0000 1.0000 1.0000 1.0000o稀疏矩阵 用元素的行列号和元素值来存储一个元素。 a=speye(100)a = (1,1) 1 (2,2) 1 (3,3) 1 (4,4) 1 (99,99) 1 (100,100) 1 b=eye(100); whos a b Name Size Bytes Class a 100 x100 1604 sparse array b 100 x100 80000 double arrayGrand total is 10100 elements using 81604 bytesl常规矩阵转换为稀疏矩阵 sp

22、arse(A) a=1 2 3;4 5 6;7 8 0; b=sparse(a)b = (1,1) 1 (2,1) 4 (3,1) 7 (1,2) 2 (2,2) 5 (3,2) 8 (1,3) 3 (2,3) 6 whos a b Name Size Bytes Class a 3x3 72 double array b 3x3 112 sparse arrayGrand total is 17 elements using 184 bytes sparse(i,j,s,n,m,nzmax) sparse(i,j,s,n,m) sparse(i,j,s,n) n,m为生成稀疏矩阵的行列数，i

23、,j,s为子矩阵，nzmax最多的非零元素数l稀疏矩阵转换为常规矩阵 full（）l通用函数形式 通用函数调用格式 funm（A，funname） funname包括sin sinh asin asinh cos cosh acos acosh tan exp log log2 pow2 sqrt abs funm(a,sqrt)ans = 3.7584 - 0.2071i -0.2271 + 0.4886i 0.3887 + 0.7700i 1.9110 - 1.0514i 0.2745 - 0.0130i 2.3243 + 0.0306i 2.0076 + 0.0483i 1.2246 -

24、 0.0659i 1.3918 - 0.2331i 1.5060 + 0.5498i 1.4884 + 0.8666i 1.4447 - 1.1833i 0.4063 + 0.4533i 2.2277 - 1.0691i 1.9463 - 1.6848i 1.2506 + 2.3006i sqrtm(a)ans = 3.7584 - 0.2071i -0.2271 + 0.4886i 0.3887 + 0.7700i 1.9110 - 1.0514i 0.2745 - 0.0130i 2.3243 + 0.0306i 2.0076 + 0.0483i 1.2246 - 0.0659i 1.39

25、18 - 0.2331i 1.5060 + 0.5498i 1.4884 + 0.8666i 1.4447 - 1.1833i 0.4063 + 0.4533i 2.2277 - 1.0691i 1.9463 - 1.6848i 1.2506 + 2.3006i sqrt(a)ans = 4.0000 1.4142 1.7321 3.6056 2.2361 3.3166 3.1623 2.8284 3.0000 2.6458 2.4495 3.4641 2.0000 3.7417 3.8730 1.0000p矩阵分解函数l特征值分解V,D=eig(X) 矩阵的特征值分解：XVVDV,D=eig

26、(X,nobalance) 关闭平衡算法的求解方法（平衡算法对于某些问题可以得到更高的精度）。V,D=eig(A,B) 广义特征值分解：AVBVD a=-149 -50 -154;537 180 546;-27 -9 -25;v,d=eig(a)v = 0.3162 0.4041 0.1391 -0.9487 -0.9091 -0.9740 0.0000 -0.1010 0.1789d = 1.0000 0 0 0 2.0000 0 0 0 3.0000 b=2 10 2;10 5 -8;2 -8 11;v,d=eig(a,b)v = 0.8211 -0.3138 -0.0191 -0.345

27、2 0.9495 -0.9441 -0.4546 -0.0044 0.3290d = 12.9030 0 0 0 -0.0045 0 0 0 0.0706l奇异值分解 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 lLU分解 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.000

28、0 0 0u = 4.0000 6.0000 7.0000 0 1.0000 -2.5000 0 0 2.5000注意：L实际上是一个“心理上”的下三角矩阵，它事实上是一个置换矩阵P的逆矩阵与一个真正下三角矩阵L1（其对角线元素为1）的乘积。l,u,p=lu(a)l = 1.0000 0 0 0.5000 1.0000 0 0.2500 0.5000 1.0000u = 4.0000 6.0000 7.0000 0 1.0000 -2.5000 0 0 2.5000p = 0 0 1 0 1 0 1 0 0 inv(p)*l*uans = 1 2 3 2 4 1 4 6 7p为置换矩阵，此时满

29、足AP-1LUlChol分解 如果A为n阶对称正定矩阵，则存在一个非奇异下三角实矩阵L，使得ALLT,当限定L的对角元素为正时，这种分解是唯一的。 a=4 -1 1;-1 4.25 2.75;1 2.75 3.5; chol(a)ans = 2.0000 -0.5000 0.5000 0 2.0000 1.5000 0 0 1.0000l正交分解 AQR 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

30、 3 0 0 -3将矩阵A做正交化分 解 ， 使 得Q*R=A，其中Q为正交矩阵（其范 数 为 1 ， 指 令norm(Q)=1)，R为对角化的上三角矩阵。矩阵的特殊操作p变维lreshape(X,M,N) X变为MN维lreshape(X,M,N,P,)X变为MNP 或reshape(X,M N P )l“：” 操作符 a=1:12; b=reshape(a,2,6)b = 1 3 5 7 9 11 2 4 6 8 10 12 c=zeros(4,3); c(:)=a(:)c = 1 5 9 2 6 10 3 7 11 4 8 12p矩阵抽取l对角元素抽取 diag（X,k）抽取X的第k条对

31、角线元素，k=0为主对角线，上对角线为正值，下对角线为负值。 diag（X） 抽取主对角线元素 a=pascal(4)a = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20 v=diag(a)v = 1 2 6 20 v=diag(a,2)v = 1 4l三角阵的抽取 tril(X)提取X的主下三角部分 tril(X,k) 提取X的第k条对角线下面的元素 triu(X)提取X的主上三角部分 triu(X,k) 提取X的第k条对角线上面的元素 al=tril(a,-1)al = 0 0 0 0 1 0 0 0 1 3 0 0 1 4 10 0 au=triu(a,-1)au

32、 = 1 1 1 1 1 2 3 4 0 3 6 10 0 0 10 20数组及其运算 数组的建立、存储完全同于矩阵，由于计算的不同，把相同型矩阵之间的运算称为数组运算。p基本数组运算l数组的四则运算 普通运算同矩阵的运算，另有点运算“.*”、”./”、”.”，即两数组对应元素之间的运算。 a1=1 2 3;2 3 4;3 4 5; b1=1 1 1;2 2 2;3 3 3; a1./b1ans = 1.0000 2.0000 3.0000 1.0000 1.5000 2.0000 1.0000 1.3333 1.6667 a1.b1ans = 1.0000 0.5000 0.3333 1.0000 0.6667 0.5000 1.0000 0.7500 0.6000l