China-bub.coM 下载 第12章 MATLAB程序设计 MATLAB有一些命令可以来控制MATLAB语句的执行，如条件语句、循环语句和支持用 户交互的命令，本章将介绍这些命令。MATLAB是一种高级的程序设计语言，能帮助用户解 决矩阵问题或其他问题。那些熟悉其他编程语言的用户，如熟悉Pascal、C++、FORTRAN等， 对理解本章内容有一定的优势。但是确信这部分内容能够让所有的读者理解和掌握。 12.1条件控制语句 MATLAB中由if语句做出判断。If语句的基本格式如下： if logical expression statements end 注意，在if和logical expression(逻辑表达式)之间要有一个空格。statement(程序语句)可 以是一个命令，也可以是由逗号、分号隔开的若千命令或者是‘returns'。只有当逻辑表达式 为tre(真)时，才能执行这些命令。逻辑表达式可以是一个标量、一个向量或者一个矩阵。如 果逻辑表达式的所有元素为非零值，它才为true(。 1f语句也可以写成一行。 if logical expression,statements,end 当然，通常前一种形式使得MATLAB程序更加结构化和易读。 ■例12.1 假设定义m×n的矩阵A。下面的语句是判断矩阵A的第1列元素是否为0，若全为0，则从 矩阵A中删除第1列： 1fA(:,1)=0 A=A(1:m,2:n) end 或者写成一行： 1fA(:,1)=0,A=A(1:m,2:n),end ■ if语句可以与elseif或else组合起来用于更复杂的上下文语句中。可能有如下的结构存在： if logical expression statements 1 else statements 2 end 如果逻辑表达式为true,则执行statements 1中的命令语句：如果为false则执行statements2 中的语句。 考虑下面的1f语句：
下载 第1 2章 M AT L A B程序设计 M AT L A B有一些命令可以来控制 M AT L A B语句的执行，如条件语句、循环语句和支持用 户交互的命令，本章将介绍这些命令。 M AT L A B是一种高级的程序设计语言，能帮助用户解 决矩阵问题或其他问题。那些熟悉其他编程语言的用户，如熟悉 P a s c a l、C + +、F O RT R A N等， 对理解本章内容有一定的优势。但是确信这部分内容能够让所有的读者理解和掌握。 12.1 条件控制语句 M AT L A B中由i f语句做出判断。I f语句的基本格式如下： 注意，在i f和logical expre s s i o n(逻辑表达式)之间要有一个空格。 s t a t e m e n t(程序语句)可 以是一个命令，也可以是由逗号、分号隔开的若干命令或者是‘ r e t u r n s’。只有当逻辑表达式 为t r u e (真)时，才能执行这些命令。逻辑表达式可以是一个标量、一个向量或者一个矩阵。如 果逻辑表达式的所有元素为非零值，它才为 t r u e (。 i f语句也可以写成一行。 当然，通常前一种形式使得 M AT L A B程序更加结构化和易读。 ■ 例1 2 . 1 假设定义m×n的矩阵A。下面的语句是判断矩阵 A的第1列元素是否为0，若全为0，则从 矩阵A中删除第1列： 或者写成一行： i f语句可以与e l s e i f或e l s e组合起来用于更复杂的上下文语句中。可能有如下的结构存在： 如果逻辑表达式为t r u e , 则执行s t a t e m e n t s 1中的命令语句；如果为f a l s e则执行s t a t e m e n t s 2 中的语句。 考虑下面的i f语句： ■
166 MATLAB5手册 China-bub.coM 下载 if logical expression I statements 1 elseif logical expression 2 statements 2 end 当logical expression 1为true时，执行statements1中的命令：如果logical expression1为 false并且logical expression2为truel时，执行statements2。 注意，elseif必须写成一个单词，如果分开写成e1seif,将会被解释成不同的意思。 命令elseif不像else if一样需要一个额外的end。 另外1f语句可以被嵌套成下面的形式： if logical expression I statements I elseif logical expression 2 statements 2 else statements 3 end 更复杂的情况如下： if logical expression I statements 1 if logical expression 2 statements 2 else statements 3 end else statements 4 end ■例12.2 ()如果A为非奇异矩阵，就能解方程Ax-b:否则要取决于扩展矩阵(Ab)的梯形形式行的 个数。提示：如果一个矩阵是方阵或为满秩的，则它为非奇异矩阵。 号给出矩阵A和方程右边b。 s size(A) if(s(1)==s(2)&(rank(A)=s(1)) x =A\b else rref ([A b]) end (b)如果矩阵A的行列式为0，则计算特征值为0的个数： if det(A)==0 length(find(eig(A)==0)) end ■
当logical expre s s i o n 1为t r u e时，执行s t a t e m e n t s 1中的命令；如果logical expre s s i o n 1为 f a l s e并且 logical expre s s i o n 2为t r u e时，执行s t a t e m e n t s 2。 注意，e l s e i f必须写成一个单词，如果分开写成 else if，将会被解释成不同的意思。 命令e l s e i f不像else if一样需要一个额外的e n d。 另外i f语句可以被嵌套成下面的形式： 更复杂的情况如下： ■ 例1 2 . 2 (a) 如果A为非奇异矩阵，就能解方程 A x = b；否则要取决于扩展矩阵(A b)的梯形形式行的 个数。提示：如果一个矩阵是方阵或为满秩的，则它为非奇异矩阵。 % 给出矩阵A和方程右边b。 (b) 如果矩阵A的行列式为0，则计算特征值为0的个数： 1 6 6 M ATLAB 5 手册 下载 ■
China-pub.CoM 第I2章MATLAB程序设计 167 下载 另一种条件语句是switch-case语句，如下： switch logical expression (scalar or string) case value 1 expression I case value 2 expression 2 ，，+ otherwise expression end logical expression经过计算给出一个标量或字符串作为结果。将这个结果与alue1,aue2,.进行比较， 如果它们匹配，则执行相应ase下的语句xpression。如果没有匹配的，则执行herwise下的语句。 如果expression的结果是一个标量，将通过检查：expression==value来决定执行的表达式。 如果表达式的结果是一个字符串，那么用strcmp(expression,value)来检查。测试结果为真， 则执行相应的表达式，而其他case语句中的表达式将不会被执行。 通过将不同的值放入细胞矩阵，就能用case语句与不同的值进行比较：见例12.3。 ■例12.3 检测掷一次骰子所得的点数是单数还是双数： function dicetest(result) switch result case{1,3,5] disp('odd number of eyes') case[2,4,6] disp('even number of eyes') otherwise disp('What kind of dice do you have?') 运行这个函数可以得到如下结果： dicetest(1) odd number of eyes dicetest(4) even number of eyes dicetest(7) What kind of dice do you have? ■ 如果表达式出错，可以使用try/catch组合，其形式如下： try expression I catch expression 2 end
另一种条件语句是s w i t c h - c a s e语句，如下： logical expre s s i o n经过计算给出一个标量或字符串作为结果。将这个结果与v a l u e1, value 2, ...进行比较， 如果它们匹配，则执行相应的c a s e下的语句e x p re s s i o n。如果没有匹配的，则执行o t h e r w i s e下的语句。 如果e x p re s s i o n的结果是一个标量，将通过检查：e x p re s s i o n= =v a l u e来决定执行的表达式。 如果表达式的结果是一个字符串，那么用 s t r c m p(e x p ression, value)来检查。测试结果为真， 则执行相应的表达式，而其他 c a s e语句中的表达式将不会被执行。 通过将不同的值放入细胞矩阵，就能用 c a s e语句与不同的值进行比较；见例 1 2 . 3。 ■ 例1 2 . 3 检测掷一次骰子所得的点数是单数还是双数： 运行这个函数可以得到如下结果： 如果表达式出错，可以使用 t r y / c a t c h组合，其形式如下： 第1 2章 M AT L A B程序设计 1 6 7 下载 ■
168 MATLAB5手册 China-pub.com 下载 MATLAB开始执行expression1,但如果有错误，错误信息将被存储在lasterr中，并且 执行expression2。 12.2 循环语句 MATLAB有两个命令for和while能反复执行语句。在逻辑控制下，这些命令能灵活地一次 或多次执行语句。 命令for与大多数的程序设计语言中的do或for命令一样。这个命令就是反复执行一条 语句或一组语句，而执行的次数已预先定义好。以end结束这组语句。 for循环通常的语法为： for variable expression statements end 象if语句一样，for语句也能写在一行上： for variable expression,statements,end 在for和variable之间需要有一个空格。这里的variable是循环变量名。在表达式中给出循 环的初始值、步长和终值。这个步长可为负数或单位值。如果为单位值，循环变量每次迭代 将增加1。通常我们用冒号来定义expression,例如i:j:k或i:j,参见4.3节。 表达式中的列值被一个一个地存放在循环变量中。因此，可以用一个矩阵来代替表达式。 例如下面的语句： for v =A,...,end 就等价于： for j=1:n,v=A(:,j);...,end 当表达式用冒号来表示时，那么列值都是标量，例如MATLAB中的语句：forv=i:j:k。 循环是可以嵌套的： for variable I =expression A statements 1 for variable II =expression B statements 2 end statements 3 end ■例12.4 (a)下列矩阵有三个非零对角值（这是一个三对角阵）： 100 0 5 10 0 0 15 0 0 0151 0 0015 这个矩阵可循环使用命令￡ox来创建。这种方法在任何标准的程序设计语言中都是一样的
M AT L A B开始执行e x p re s s i o n 1，但如果有错误，错误信息将被存储在 l a s t e r r中，并且 执行e x p re s s i o n 2。 12.2 循环语句 M AT L A B有两个命令f o r和w h i l e能反复执行语句。在逻辑控制下，这些命令能灵活地一次 或多次执行语句。 命令f o r与大多数的程序设计语言中的 d o或f o r命令一样。这个命令就是反复执行一条 语句或一组语句，而执行的次数已预先定义好。以 e n d结束这组语句。 f o r循环通常的语法为： 象i f语句一样，f o r语句也能写在一行上： 在f o r和v a r i a b l e之间需要有一个空格。这里的 v a r i a b l e是循环变量名。在表达式中给出循 环的初始值、步长和终值。这个步长可为负数或单位值。如果为单位值，循环变量每次迭代 将增加1。通常我们用冒号来定义e x p re s s i o n，例如i : j : k或i : j，参见4 . 3节 。 表达式中的列值被一个一个地存放在循环变量中。因此，可以用一个矩阵来代替表达式。 例如下面的语句： 就等价于： 当表达式用冒号来表示时，那么列值都是标量，例如M AT L A B中的语句：for v=i:j: k。 循环是可以嵌套的： ■ 例1 2 . 4 (a) 下列矩阵有三个非零对角值 (这是一个三对角阵) : 这个矩阵可循环使用命令f o r来创建。这种方法在任何标准的程序设计语言中都是一样的。 1 6 8 M ATLAB 5 手册 下载
China-pub.CoM 第I2章MATLAB程序设计 169 下载 A=[门； for k 1:5 for j=1:5 if k =j A(k,k)=5; elseif abs(k-j)==1 A(k,j)=1; else A(k,j=; end end end 这里的分号’；是非常重要的。如果这些赋值语句没有分号，矩阵A将在屏幕上输出25次， 每一次A中的元素将被赋值一次。 同样也会遇见由于不注意使用￡。r循环而导致无效操作的例子。通过定时器时钟就能清 楚地计算出花费的时间。例如，for循环，见例12.21。下面的命令能完成上面同样的事情， 并且更加有效： A zeros(5); for k 1:4 A(k,k)=5; A(k,k+1)=1; A(k+1,k)=1; end A(5,5)=5; 这个矩阵能通过更加快速有效的方法得到，但是使用命令diag更加清楚。 A=门； A diag(5*ones(5,1))+diag(ones(4,1),1)+... diag(ones(4,1),-1); 这种结构的大矩阵应该创建成稀疏矩阵：参见第9章。 (b)在区间[一2，一0.75]内，步长为0.25，对函数=fx)=1+1/x求值，并列表。将所得x值和 y值分别存入向量r和s中，并列表显示： r=门；8=[门； forx=-2.0:0.25:-0.75 y=1+1/x; r=[r x]; 8=[sy]; end [r;s]
这里的分号’ ; ’是非常重要的。如果这些赋值语句没有分号，矩阵 A将在屏幕上输出 2 5次， 每一次A中的元素将被赋值一次。 同样也会遇见由于不注意使用 f o r循环而导致无效操作的例子。通过定时器时钟就能清 楚地计算出花费的时间。例如， f o r循环，见例 1 2 . 2 1。下面的命令能完成上面同样的事情， 并且更加有效： 这个矩阵能通过更加快速有效的方法得到，但是使用命令 d i a g更加清楚。 这种结构的大矩阵应该创建成稀疏矩阵；参见第 9章。 (b) 在区间[－2 ,－0 . 7 5 ]内，步长为0 . 2 5，对函数y=f(x) = 1 + 1 /x求值，并列表。将所得 x值和 y值分别存入向量r和s中，并列表显示： 第1 2章 M AT L A B程序设计 1 6 9 下载
170 MATLAB5手册 China-bub.coM 下载 此表也能够不用for循环语句创建，其结果为： ans -2.0000 0.5000 -1.7500 0.4286 -1.5000 0.3333 -1.2500 0.2000 -1.0000 0 -0.7500 -0.3333 (c)MATLAB命令sum(A)给出一行向量，向量的元素是矩阵A每一列元素的和。用下面 的程序能得到相近的结果： A=[123;456] sum_v=[] for v=A sum_v [sum_v sum(v)] end disp('Compare w/sum(A):') disp(sum(A)); 其结果是： A= 1 2 3 4 5 6 sum_v 5 sun_v 5 7 sumv 7 9 Compare w/sum(A): 5 7 9 (d)将下列MATLAB命令保存在文件qrmethod.m中： 号矩阵A及整数m和n应在调用该文件之前定义好。 号在变换成上海森伯形式之后使用Q方法。 告n步以后结束。 号每隔m步输出结果。 A hess(A); for i =:n [Q,R]gr(A); A =R*Q; nd =norm(diag(A,-1));
此表也能够不用f o r循环语句创建，其结果为： (c) MATLAB 命令s u m ( A )给出一行向量，向量的元素是矩阵 A每一列元素的和。用下面 的程序能得到相近的结果： 其结果是： (d) 将下列M AT L A B命令保存在文件q r m e t h o d . m中： % 矩阵A及整数m和n应在调用该文件之前定义好。 % 在变换成上海森伯形式之后使用 Q R方法。 % n步以后结束。 % 每隔m步输出结果。 1 7 0 M ATLAB 5 手册 下载
China-pub.CoM 第I2章MATLAB程序设计 171 下载 if rem(i,m)==0 A,i,nd end end 下面的程序是用来完成非移位的QR方法（见8.2节）30次迭代，每隔15次输出结果： A0=[-9-3-16;13716;3310]; m=15;n=30; format long; A=A0; qrmethod; A= 9.9899746707437722.62301237506363-15.53274662438004 0.00708686385759-5.985685125529255,77401643542405 00.00741470005235 3.99571045478546 i= 15 nd 0.01025677416162 A= 10.0000047162466022.6274399374496715.51339551121122 -0.00000333488655-6.00001449452640-5.77348898879412 0 0.000016936546124.00000977827978 i= 30 nd 1.726175143943722e-05 (e)在7.5节中定义了命令planerot。这个算法是使用该命令来返回一个矩阵，该矩阵的 0元素都在任何作为参数输入的m×n矩阵的主对角线的下方。 function B Givens(A) 号如果将A矩阵乘函数B返回的矩阵，就能将大小为m×的A矩阵变成上三角阵。也就是说，根据Q=和· R=B+A,其中满足Q+R=A,该函数可以对A进行QR分解。 [m,n]size(A); =eye(m); for j=1:n for i j:m for k (i+1):m G =eye(m); Plan planerot([A(j,j)A(k,j)]'); %找出2×2的矩阵 G(j,j)=P1an(1,1);G(k,j)=P1an(2,1); G(j,k)=P1an(1,2);G(k,k)=P1an(2,2);
下面的程序是用来完成非移位的 Q R方法(见8 . 2节) 3 0次迭代，每隔1 5次输出结果： (e) 在7 . 5节中定义了命令p l a n e r o t。这个算法是使用该命令来返回一个矩阵，该矩阵的 0元素都在任何作为参数输入的 m×n矩阵的主对角线的下方。 function B = Givens(A) % % 如果将 A矩阵乘函数 B返回的矩阵，就能将大小为 m×n的A矩阵变成上三角阵。也就是说，根据 Q = B ´ 和 R = B*A，其中满足Q*R = A，该函数可以对A进行Q R分解。 第1 2章 M AT L A B程序设计 1 7 1 下载 ％找出2×2的矩阵
172 MATLAB5手册 China-bub.coM 下载 %在mXm矩阵中正确定位 B=G*B; A G*A %<-To see the step-by-step reduction of AA.remove this semicolon. end end end 在这个算法中，每一步的循环将A中的两个元素取出构成2×2的矩阵，然后将矩阵主对角 线下的一个元素赋值为零。这个结果矩阵可用来创建QR因式分解。 定义一个检测矩阵Atest: 2312 4 1 2 Atest 6 32 2100 下面的命令为： Atest=[123123;441221; 763211;121002]; Giv Givens(Atest); Q Giv',R=Giv*Atest QR =Q*R %仅仅是检测 给出MATLAB输出： Q= 0.1222 0.6630 0.6674 0.3162 0.4887 0.1842 -0.5721 0.6325 0.8552 -0.2947 0.2860 -0.3162 0.1222 0.6630 -0.3814 -0.6325 R= 8.1854 7.5745 3.5429 2.8099 2.0769 1.9547 0.0000 1.6208 1.9523 0.4420 1.3997 3.2047 -0.0000 -0.0000 1.9069 0.0953 0.4767 0.9535 0.0000 0.0000 0.0000 0.9487 1.5811 0.0000 QR 1.0000 2.0000 3.0000 1.0000 2.0000 3.0000 4.0000 4.0000 1.0000 2.0000 2.0000 1.0000 7.0000 6.0000 3.0000 2.0000 1.0000 1.0000 1.0000 2.0000 1.0000 0 0 2.0000 注意，通过乘Q和R我们又能得到初始矩阵Atest,因此QR=Atest。. ()以下程序通过使用两个for循环和平面组合来画出雪花图形。这个算法生成Helge von Koch曲线，这是一个不规则碎片例子。程序中使用的图形命令定义在第13章中，这里的注释 主要说明其功能。该算法将当前几何图形每一面分成了相同的三个部分。第一部分和最后部
在这个算法中，每一步的循环将 A中的两个元素取出构成2×2的矩阵，然后将矩阵主对角 线下的一个元素赋值为零。这个结果矩阵可用来创建 Q R因式分解。 定义一个检测矩阵A t e s t： 下面的命令为： 给出M AT L A B输出： 注意，通过乘Q和R我们又能得到初始矩阵A t e s t，因此Q R=A t e s t。 (f) 以下程序通过使用两个 f o r循环和平面组合来画出雪花图形。这个算法生成 Helge von K o c h曲线，这是一个不规则碎片例子。程序中使用的图形命令定义在第 1 3章中，这里的注释 主要说明其功能。该算法将当前几何图形每一面分成了相同的三个部分。第一部分和最后部 1 7 2 M ATLAB 5 手册 下载 % 在m×m矩阵中正确定位 ％ 仅仅是检测
China-pub.Com 第I2章MATLAB程序设计 173 下载 分是新几何学的两个方面。中间的部分用等边的三角形的两个边替代，如图12-1所示。 如果将迭代进行下去，几乎平面的每一部分都将被覆盖到。事实上，不规则碎片的尺寸 为1.2619，比1大一点而比2小。 %文件：Koch.m %该程序画出Helge von Koch雪花，一个不规则碎片图形 clear; %时除旧变量 %向量在平面中新定义一个三角形。这是开始的几何状态。 new=[0.5+(sqrt(3)/2)*i,-0.5+(sqrt(3)/2)*1,.. 0,0.5+(sqrt(3)/2)*i]; plot(new); %画出三角形并等待0.5秒 pause(0.5); %迭代5次： 向量old是前一次迭代 for k=1:5; old new; [m,n]size(old); n=n-1; %o1d定义了图中的n-1条边。 %对每条边：定义4个新点（其中一个是'o1d')。 for j=0:n-1; diff=(o1d(j+2)-o1d(j+1)/3; new(4*j+1)=o1d(j+1); new(4*j+2)=old(j+1)+diff; new(4*j+3)=new(4*j+2)+diff*((1-sqrt(3)*i)/2); new(4*j+4)=old(j+1)+2*diff; end; %向量new的最后一个元素与向量old的最后一个元素相同。 new(4*n+1)÷o1d(n+1); plot(new); %画出新图并等待0.5秒。 pause(0.5); end; %移开坐标轴，并使其等长度，图形会更匀称。 axis off;axis square; 执行程序，就可得到一个逐渐复杂的图形。 图12-2给出了最后的图形
分是新几何学的两个方面。中间的部分用等边的三角形的两个边替代，如图 1 2 - 1所示。 如果将迭代进行下去，几乎平面的每一部分都将被覆盖到。事实上，不规则碎片的尺寸 为1 . 2 6 1 9，比1大一点而比2小。 执行程序，就可得到一个逐渐复杂的图形。 图1 2 - 2给出了最后的图形。 第1 2章 M AT L A B程序设计 1 7 3 ％ 文件 ％ 该程序画出Helge von Koch 雪花，一个不规则碎片图形 ％删除旧变量 ％向量在平面中新定义一个三角形。这是开始的几何状态。 ％画出三角形并等待0.5秒 ％ 迭代5次： 向量old是前一次迭代 ％ old定义了图中的n-1条边。 ％ 对每条边：定义4个新点(其中一个是'old')。 ％ 向量new的最后一个元素与向量old的最后一个元素相同。 ％画出新图并等待0.5秒。 ％ 移开坐标轴，并使其等长度，图形会更匀称。 下载
174 MATLAB5手份 China-bub.coM 下载 图12-25次迭代后的Helge von Koch不规则碎片 图I2-I von Koch算法对直线的两次迭代的结果图 图形，其原始几何图形为三角形 ■ 只要逻辑表达式为真，whi1e命令将反复执行程序语句。像for语句一样程序体由一个 end来结束。使用while循环来表示整个while语句，具体形式如下： while,statements,end 通常，while循环有如下形式： while logical expression statements end 将其写在一行的形式为： while logical expression,statemenesd while循环能够像for循环一样嵌套： while logical expression A statements 1 while logical expression B statements 2 end statements 3 end ■例12.5 (a)构造一个特征值在1和一1之间的2×2的随机矩阵，可以用下面的迭代来实现： A=rand(2); 构造一个特征值在1和-1之间的矩阵。 while max(abs(eig(A)))>=1 A rand(2); end e eig(A); TheText ['lambda_1=',num2str(e(1)),... ’,1 ambda_2=’,num28tr(e(2)]; A 输出显示矩阵及其特征值
只要逻辑表达式为真， w h i l e命令将反复执行程序语句。像 f o r语句一样程序体由一个 e n d来结束。使用w h i l e循环来表示整个w h i l e语句，具体形式如下： while, s t a t e m e n t s, end 通常，w h i l e循环有如下形式： 将其写在一行的形式为： while logical expression, statements , end w h i l e循环能够像f o r循环一样嵌套： ■ 例1 2 . 5 (a) 构造一个特征值在1和－1之间的2×2的随机矩阵，可以用下面的迭代来实现： A = r a n d ( 2 ) ; % 构造一个特征值在1和- 1之间的矩阵。 A % 输出显示矩阵及其特征值。 1 7 4 M ATLAB 5 手册 下载 图12-1 von Koch算法对直线的两次迭代的结果图 图12-2 5次迭代后的Helge von Koch不规则碎片 图形，其原始几何图形为三角形 ■