《matlab举例.doc》由会员分享,可在线阅读,更多相关《matlab举例.doc(17页珍藏版)》请在三一文库上搜索。
1、Matlab特色举例(2009-03-11 16:05:31)标签:杂谈 分类:科学研究考虑两个矩阵 A 和 B 的乘积问题,在 C 语言中要实现两个矩阵的乘积并不仅仅是一组双重循环的问题。双重循环当然是矩阵乘积所必需的,除此之外要考虑的问题很多。例如:A 和 B 有一个是复数矩阵怎么考虑;其中一个是复数矩阵时怎么考虑;全部是实系数矩阵时又怎么管理;这样就要在一个程序中有 4 个分支,分别考虑这 4 种情况。然后还得判断这两个矩阵是否可乘。而考虑两个矩阵是否可乘也并不仅仅是判断 A 的列数是否等于 B 的行数这么简单。其中一个若为标量,则它们可以无条件地相乘。其中有标量时又得考虑实数与复数的问
2、题等。所以说,没有几十分钟的时间,用 C 语言并不可能编写出考虑各种情况的子程序。有了 MATLAB 这样的工具,A 和 B 矩阵的乘积用 A*B 这样简单的算式就能表示了。例 1-1矩阵生成与运算。考虑金庸作品中经常提及的一个“数学问题”, 该问题用半数学语言描述就是:如何生成一个 3x3 矩阵, 并将自然数 1, 2, ., 9 分别置成这 9 个矩阵元素,才能使得每一行、每一列、且主、反对角线上元素相加都等于一个相同的数。 这样的矩阵称为“魔方矩阵”。用 MATLAB 的 magic() 函数,我们可以由下面的命令立即生成这样的矩阵: A=magic(3)A = 8 1 6 3 5 7
3、4 9 2还可以由 B=magic(10) 一次生成 10x10 的魔方矩阵。如果想求出矩阵的行列式和特征值,可以分别由 det(B) 与 eig(B) 立即得出结果,而同样的工作在 C 下并不是很简单就可以得出的,算法选择不好,还可能得出错误的结果。例 1-2考虑一个二元函数如何用三维图形的方式表现出这个曲面? 用 C 这类语言,绘制图形是一个难点,且从一个机器移植程序到另一个机器,大部分调试程序时间都花在这上。但使用 MATLAB 这类高级语言,完成这样的工作就是几个直观语句的事。且得出的图形美观准确、可以将语句毫不变化地移植到另外的机器上,得出完全一致的结果,如下所示。 x,y = me
4、shgrid(-3:1/8:3); z = 3*(1-x).2.*exp(-(x.2) - (y+1).2)- 10*(x/5 - x.3 - y.5). .*exp(-x.2-y.2)- 1/3*exp(-(x+1).2 - y.2); surf(x,y,z), shading interp; colorbar例 1-3微分方程的数值解法是在科学与工程计算中经常遇到的问题。假设著名的 Lorenz 模型的状态方程表示为:若令 且初值为,e 为一个小常数,假设 则我们可以由下面的几个语句就可以描述微分方程:function xdot = lorenzeq(t,x)xdot=-8/3*x(1)+
5、x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3);这样下面几个语句就能求解该微分方程,绘制出时间曲线与相空间曲线,如下所示。 t_final=100; x0=0;0;1e-10; t,x=ode45(lorenzeq,0,t_final,x0); plot(t,x), figure; plot3(x(:,1),x(:,2),x(:,3); axis(10 40 -20 20 -20 20);例 1-5 (注,这里的编号采用作者书中的序号) 设有解析函数,利用 MATLAB 的符号运算工具箱可以对该函数进行解析推导,得出诸如高阶导数、积分、
6、Taylor 幂级数展开等。 syms x; f=x2*(sin(x)2;diff(f); f1=simple(ans)f1 = x-x*cos(2*x)+x2*sin(2*x) diff(f,x,2); f2=simple(ans)f2 = 1-cos(2*x)+4*x*sin(2*x)+2*x2*cos(2*x) diff(f,x,3); f3=simple(ans)f3 = 6*sin(2*x)+12*x*cos(2*x)-4*x2*sin(2*x) diff(f,x,4); f4=simple(ans)f4 = 24*cos(2*x)-32*x*sin(2*x)-8*x2*cos(2*
7、x) int(f4,x)ans = 6*sin(2*x)+12*x*cos(2*x)-4*x2*sin(2*x) taylor(x2*(sin(x)2,15,x)ans = x4-1/3*x6+2/45*x8-1/315*x10+2/14175*x12-2/467775*x14Matlab支持的数据结构 MATLAB 语言的赋值语句有两种: 变量名 = 运算表达式 返回变量列表 = 函数名(输入变量列表) MATLAB 支持变量和常量,其中 pi 为圆周率 p, 更重要的,MATLAB 支持 IEEE 标准的运算符号,如 Inf 表示无穷大,NaN (Not a Number) 为 0/0,
8、0*Inf 或 Inf/Inf 等运算结果。MATLAB 变量名应该由字母引导,后面可以跟数字、字母或下划线等符号。MATLAB 是区分变量名字母大小写的。(1) 矩阵 MATLAB 最基本的数据结构是复数矩阵。输入一个复数矩阵是很简单的事。例如可以给出下面的语句: B=1+9i,2+8i,3+7j; 4+6j 5+5i,6+4i; 7+3i,8+2j 1i其中 为 MATLAB 的提示符。矩阵各行元素由分号分隔,而同行不同元素由逗号或空格分隔。给出了上面的命令,则可以给出下面的结果。B = 1.0000 + 9.0000i 2.0000 + 8.0000i 3.0000 + 7.0000i
9、4.0000 + 6.0000i 5.0000 + 5.0000i 6.0000 + 4.0000i 7.0000 + 3.0000i 8.0000 + 2.0000i 0 + 1.0000i其中,元素 1+9i 表示复数项。有这样的表述方法,实矩阵、向量或标量均可以更容易地输入了。如果赋值表达式末尾有分号,则其结构将不显示,否则将显示出全部结果。 MATLAB 和其他语言不同,它无需事先声明矩阵的维数。下面的语句可以建立一个更大的矩阵 B(2,5)=1B = 1.0000 + 9.0000i 2.0000 + 8.0000i 3.0000 + 7.0000i 0 0 4.0000 + 6.0
10、000i 5.0000 + 5.0000i 6.0000 + 4.0000i 0 1.0000 7.0000 + 3.0000i 8.0000 + 2.0000i 0 + 1.0000i 0 0 冒号表达式是 MATLAB 里最具特色的表示方法。其调用格式为 a=s1:s2:s3; 这一语句可以生成一个行向量,其中 s1 为向量的起始值,s2 为步距,而 s3 为向量的终止值。例如 S=0:.1:2*pi; 将产生一个起始于 0, 步距为 0.1, 而终止于 6.2 的向量 (pi 为 MATLAB 保留常量p), 而不是终止于2p。如果写成 S=0:-0.1:2*pi; 则不出现错误,而返回
11、一个空向量。 冒号表达式可以用来提取矩阵元素,例如 B(:,1) 将提取 B 矩阵的第 1 列而 B(1:2,1:2:3) 将提取 B 的前 2 行与 1,3,5 列组成的子矩阵。在矩阵提取时还可以采用end 这样的算符。如 B(2:end,:) 将提取 B 矩阵的后 2 列构成的子矩阵。(2) 多维数组 多维数组是 MATLAB 在其 5.0 版本开始提供的。假设有 2 个 3x3 矩阵 A1, A23,则可以由下面的命令建立起一个 3x3x2 的数组:A=cat(3,A1,A2)。试验 A1=cat(2,A1,A2) 和 A2=cat(1,A1A2) 将得到什么结果。 对矩阵或多维数组 A
12、 可以使用 size(A) 来测其大小,也可以使用 reshape() 函数重新按列排列。对向量来说,还可以用 length(A) 来测其长度。 不论原数组 A 是多少维的,A(:) 将返回列向量。(3) 字符串与字符串矩阵 MATLAB 的字符串是由单引号括起来的。如可以使用下面的命令赋值 strA=This is a string. 多个字符串可以用 str2mat() 函数构造出字符串矩阵。如 B=str2mat(strA, ksa saj,aa);字符串变量可以由下表中的命令进行操作:命令意义命令意义strcmp(A,B)比较A和B字符串是否相同。findstr(A,B)测试A是否为B
13、的子字符串,或反过来strrep(A,s1,s2)在A中用s2替换s1length(A)字符串A的长度deblank(A)删除A字符串尾部的空格double(A)字符串转换双精度数据(4) 单元数据结构 用类似矩阵的记号将给复杂的数据结构纳入一个变量之下。和矩阵中的圆括号表示下标类似,单元数组由大括号表示下标。 B=1,Alan Shearer,180,100, 80, 75; 77, 60, 92; 67, 28, 90; 100, 89, 78B = 1 Alan Shearer 180 4x3 double访问单元数组应该由大括号进行,如第 4 单元中的元素可以由下面的语句得出 B4an
14、s = 100 80 75 77 60 92 67 28 90 100 89 78(5) 结构体 MATLAB 的结构体有点象 C 语言的结构体数据结构。每个成员变量用点号表示,如 A.p 表示 A 变量的 p 成员变量。获得该成员比 C 更直观,仍用 A.p 访问,而不用 A-p。用下面的语句可以建立一个小型的数据库。 student_rec.number=1;student_rec.name=Alan Shearer;student_rec.height=180;student_rec.test=100, 80, 75; 77, 60, 92; 67, 28, 90; 100, 89, 7
15、8; student_recstudent_rec = number: 1 name: Alan Shearer height: 180 test: 4x3 double其中 test 成员为单元型数据。删除成员变量可以由 rmfield() 函数进行,添加成员变量可以直接由赋值语句即可。另外数据读取还可以由 setfield 和 getfield 函数完成。(6) 类与对象 类与对象是 MATLAB 5.* 开始引入的数据结构。在 MATLAB 手册中定义了一各很好的类 - 多项式类。该例子值得细读,去体会类和对象的定义,重载函数编写等信息。事实上,在实际工具箱设计中,用到了很多的类,例如在
16、控制系统工具箱中定义了 LTI (线性时不变系统) 类,并在此基础上定义了其子类:传递函数类 TF, 状态方程类 SS, 零极点类 ZPK 和频率响应类 FR。举例:我们将通过一个例子来介绍类的构造。 在 MATLAB 语言使用手册中给出了一个很有代表性的例子:多项式类的建立问题。假设我们想为多项式建立一个单独的类,重新定义加、减、乘及乘方等运算,并定义其显示方式。那么建立一个类至少应该执行下面的步骤:(这个例子更详细的情况请参考 MATLAB 手册) 首先应该选定一个恰当的名字,例如这里的多项式类可选择为 polynom。 以这个名字建立一个子目录,目录的名字前加 。对本例来说,即应该在当前
17、的工作目录下建立 polynom 子目录,而这个目录无需在 MATLAB 路径下再指定。 编写一个引导函数,函数名应该和类同名。定义类的使用方法:function p = polynom(a)if nargin = 0p.c = ; p = class(p,polynom);elseif isa(a,polynom), p = a;else,p.c = a(:).; p = class(p,polynom);end可以看出,本函数分三种情况加以考虑: 如果不给输入变量,则建立一个空的多项式; 如果输入变量 a 已经为多项式类,则将它直接传送给输出变量 p; 如果 a 为向量,则将此向量变换成行
18、向量,再构造成一个多项式对象。 如果想正确地显示新定义的类,则必需首先定义 display() 函数,并对新定义的类重新定义其基本运算。对多项式来说,我们可以如下定义有关的函数: 要改变显示函数的定义,则需在此目录下重新建立一个新函数 display()。这种重新定义函数的方法又称为函数的重载。显示函数可以如下地重载定义。function display(p)disp( ); disp(inputname(1), = )disp( ); disp( char(p); disp( );注意,这里应该定义的是 display() 而不是 disp()。 从上面的定义可见,显示函数要求重载定义 ch
19、ar() 函数,用于把多项式转换成可显示的字符串。 该函数的定义为function s=char(p)if all(p.c=0), s =0;elsed=length(p.c)-1; s=;for a=p.c;if a=0;ifisempty(s)if a0, s=s, + ;else, s=s, - ; a = -a; endendif a=1 | d=0, s=s, num2str(a);if d0, s=s, *; endendif d=2, s=s, x, int2str(d);elseif d=1, s=s x; endendd=d-1;end, end 仔细研究此函数,可以发现,该
20、函数能自动地按照多项式显示的格式构造字符串。比如,多项式各项用加减号连接,系数与算子之间用乘号连接,而算子的指数由 表示。再配以显示函数,则可以将此多项式以字符串的形式显示出来。 双精度处理:双精度转换函数的重载定义是很简单的。function c = double(p)c = p.c; 加运算:两个多项式相加,只需将其对应项系数相加即可。这样,加法运算的重载定义可由下面的函数实现。注意,这里要对 plus() 函数进行重载定义。function p=plus(a,b)a=polynom(a); b=polynom(b);k=length(b.c)-length(a.c);p=polynom(
21、zeros(1,k) a.c+zeros(1,-k) b.c);同理,还可以重载定义多项式的减法运算:function p=minus(a,b)a=polynom(a); b=polynom(b);k=length(b.c)-length(a.c);p=polynom(zeros(1,k) a.c-zeros(1,-k) b.c); 乘法运算:多项式的乘法实际上可以表示为系数向量的卷积,可以由 conv() 函数直接获得。故可以如下重载定义多项式的乘法运算。function p=mtimes(a,b)a=polynom(a); b=polynom(b); p=polynom(conv(a.c,
22、b.c); 乘方运算: 多项式的乘方运算只限于正整数乘方的运算,其 n 次方相当于将该多项式自乘 n 次。若 n=0,则结果为 1。这样我们就可以重载定义多项式的乘方运算为:function p=mpower(a,n)if n=0, n=floor(n); a=polynom(a); p=1;if n=1,for i=1:n, p=p*a; endendelse, error(Power should be a non-negative integer.)end 多项式求值问题:可以对多项式求值函数 polyval() 进行重载定义。function y=polyval(a,x)a=polyn
23、om(a); y=polyval(a.c,x);定义了此类之后,我们就可以方便地进行多项式处理了。例如我们可以建立两个多项式对象 P(s)=x3+4x2-7 和 Q(s)=5x4+3x3-1.5x2+7x+8 其相应的MATLAB 语句为 P=polynom(1,4,0,-7), Q=polynom(5,3,-1.5,7,8)P =x3 + 4*x2 - 7Q =5*x4 + 3*x3 - 1.5*x2 + 7*x + 8然后调用下面函数就可以得出相应的计算结果 P+Qans =5*x4 + 4*x3 + 2.5*x2 + 7*x + 1 P-Qans =-5*x4 - 2*x3 + 5.5*
24、x2 - 7*x - 15 P*Qans =5*x7 + 23*x6 + 10.5*x5 - 34*x4 + 15*x3 + 42.5*x2 - 49*x - 56 X=P3X =x9 + 12*x8 + 48*x7 + 43*x6 - 168*x5-336*x4+147*x3+588*x2-343 y=polyval(X,1 2 3 4 5 6)y =-8 4913 175616 1771561 10360232 43986977由于前面的重载定义,下面的表达式也能得出期望的结果 P+1 2 3ans =x3 + 5*x2 + 2*x - 4 使用 methods() 函数可以列出一个新的类
25、已经定义的方法函数名。 methods(polynom)Methods for class polynom:char double mpower plus polyvaldisplay minus mtimes polynom变量的运算(1) MATLAB 变量的代数运算 如果给定两个矩阵 A 和 B, 则我们可以用 A+B, A-B, A*B 可以立即得出其加、减和乘运算的结果。若这两个矩阵数学上不可以这样运算,则将得出错误信息,并终止正在运行的程序。 在 MATLAB 下,如果 A 和 B 中有一个是标量,则可以无条件地进行这样的运算。MATLAB 不介意这些变量是纯实数还是含有虚部的复数
26、。 矩阵的除法实际上就是线性方程的求解,如 Ax=B 这一线性方程的解即为 x=inv(A)*B, 或更简单地 x=AB。这又称为矩阵的左除,而 x=B/A 称为矩阵的右除。 方阵的乘方可以由 算符直接得出,如 An。用 MATLKAB 这样的语言,你可以轻易地算出 A0.1, 亦即 A 矩阵开 10 次方得出的主根。 矩阵的点运算也是相当重要的。所谓点运算即两个矩阵相应元素的元素,如 A.*B 得出的是 A 和 B 对应元素的积,故一般情况下 A*B 不等于 A.*B。矩阵的点乘又称为其 Hadamard 积。点运算的概念又可以容易地用到点乘方上,例如 A.2, A.A 等都是可以接受的运算
27、式子。 Kronecker 乘积是 MATLAB 在矩阵运算中的另一个有意义的问题,用 kron(A,B) 立即可以得出两个矩阵的 Kronecker 乘积。(2) 逻辑运算 MATLAB 并没有单独定义逻辑变量。在 MATLAB 中,数值只有 0 和“非 0” 的区分。非 0 往往被认为是逻辑真,或逻辑 1。除了单独两个数值的逻辑运算外,还支持矩阵的逻辑运算,如 A & B, A | B, 和 A 分别表示逻辑与、或、非的运算。例如,下面的 A 和 B 矩阵与运算将得出如下结果 A=0 2 3 4;1 3 5 0; B=1 0 5 3;1 5 0 5; A&B ans = 0 0 1 1 1
28、 1 0 0(3) 关系表达式与表达式函数 MATLAB 的大于、小于和等于等关系分别由 、 A=0 2 3 4;1 3 5 0; B=1 0 5 3;1 5 0 5; A=B ans = 0 0 0 0 1 0 0 0确实使得 A 和 B 对应元素相等的位将返回 1,否则返回 0。MATLAB 还可以用 = 和 mysum=0; for i=1:1:100,mysum=mysum+i; end; mysummysum =5050在上面的式子中,可以看到 for 循环语句中 s3 的值为 1。在 MATLAB 实际编程中,如果 s3 的值为 1,则可以在该语句中省略,故该语句可以简化成 for
29、 i=1:100。在实际编程中,在 MATLAB 下采用循环语句会降低其执行速度,所以前面的程序可以由下面的命令来代替: i=1:100; mysum=sum(i)。在这一语句中,首先生成了一个向量 i, 然后用内部函数 sum() 求出 i 向量的各个元素之和,或更简单地,该语句还可以写成 sum(1:100)。如果前面的 100 改成 10000, 再运行这一程序,则可以明显地看出,后一种方法编写的程序比前一种方法快得多。MATLAB 并不要求循环点等间距,假设 V 为任意一个向量,则可以用 for i=V 来表示循环。同样的问题在 while 循环结构下可以表示为mysum = 0; i
30、=1; while (i 下键入该 M 文件的文件名,这样 MATLAB 就会自动执行该 M 文件中的各条语句,并将结果直接返回到 MATLAB 的工作空间。M 函数格式是 MATLAB 程序设计的主流,一般情况下,不建议您使用 M 脚本文件格式编程。 MATLAB 的 M 函数是由 function 语句引导的,其基本格式如下:function 返回变量列表 = 函数名 (输入变量列表)注释说明语句段, 由 % 引导输入、返回变量格式的检测函数体语句 这里输入和返回变量的实际个数分别由 nargin 和 nargout 两个 MATLAB 保留变量来给出,只要进入该函数,MATLAB 就将自
31、动生成这两个变量,不论您是否直接使用这两个变量。返回变量如果多于 1 个,则应该用方括号将它们括起来,否则可以省去方括号。输入变量和返回变量之间用逗号来分割。注释语句段的每行语句都应该由百分号 % 引导,百分号后面的内容不执行,只起注释作用。用户采用 help 命令则可以显示出来注释语句段的内容。此外,正规的变量个数检测也是必要的。如果输入或返回变量格式不正确,则应该给出相应的提示。我们将通过下面的例子来演示函数编程的格式与方法。例 3-假设我们想生成一个 nxm 阶的 Hilbert 矩阵, 它的第 i 行第 j 列的元素值为 1/(i+j-1)。我们想在编写的函数中实现下面几点: 如果只给
32、出一个输入参数,则会自动生成一个方阵,即令 m=n 在函数中给出合适的帮助信息,包括基本功能、调用方式和参数说明 检测输入和返回变量的个数,如果有错误则给出错误信息 如果调用时不要求返回变量,则将显示结果矩阵。其实在编写程序时养成一个好的习惯,无论对程序设计者还是对程序的维护者、使用者都是大有裨益的。 采用 MATLAB 函数编写格式和上述要求,我们可以编写出一个函数function A=myhilb(n, m)%MYHILB a demonstrative M-function.% A=MYHILB(N, M) generates an N by M Hilbert matrix A.% A
33、=MYHILB(N) generates an N by N square Hilbert matrix.% MYHILB(N,M) displays ONLY the Hilbert matrix, but do not return any% matrix back to the calling function.%See also: HILB.% Designed by Professor Dingyu XUE, Northeastern University, PRC% 5 April, 1995, Last modified by DYX at 21 March, 2000if na
34、rgout1, error(Too many output arguments.); endif nargin=1, m=n;elseif nargin=0 | nargin2error(Wrong number of iutput arguments.);endA1=zeros(n,m);for i=1: nfor j=1:mA1(i,j)=1/(i+j-1);end, endif nargout=1, A=A1; elseif nargout=0, disp(A1); end这样规范编写的函数用 help 命令可以显示出其帮助信息: help myhilbMYHILB a demonstr
35、ative M-function.A=MYHILB(N, M) generates an N by M Hilbert matrix A.A=MYHILB(N) generates an N by N square Hilbert matrix.MYHILB(N,M) displays ONLY the Hilbert matrix, but do not return anymatrix back to the calling function.See also: HILB. 有了函数之后,可以采用下面的各种方法来调用它,并产生出所需的结果。 A=myhilb(3,4)A =1.0000 0
36、.5000 0.3333 0.25000.5000 0.3333 0.2500 0.20000.3333 0.2500 0.2000 0.1667 A=myhilb(4)A =1.0000 0.5000 0.3333 0.25000.5000 0.3333 0.2500 0.20000.3333 0.2500 0.2000 0.16670.2500 0.2000 0.1667 0.1429 myhilb(4)1.0000 0.5000 0.3333 0.25000.5000 0.3333 0.2500 0.20000.3333 0.2500 0.2000 0.16670.2500 0.2000
37、 0.1667 0.1429MATLAB 工具箱编写技巧 放入一个目录中的为某种目的专门编写的一组 MATLAB 函数就可以组成一个工具箱。从某种意义上说,任何一个 MATLAB 语言的使用者都可以是工具箱的作者。在一个工具箱中,应该有一个名为 Contents.m 的文件,用来描述工具箱中所有 MATLAB 函数的名称和意义。在该文件中第 1 行应该给出该工具箱的名称,在第 2行中给出该工具箱的版本与修改时间等信息。然后分类地给出该工具箱中各类函数的最基本功能。注意,本文件中所有的语句都应该是注释语句,由百分号 % 引导,空行也应该由 % 引导。 另外,因为 MATLAB是一种解释性语言,所以即使在某个或某些函数中存在语法错误,但如果没执行到该语句时可能就不会发现该错误,这在一个成功的程序设计中是不能容许的。要查出某目录中所有的M函数语法错误,首先应该用 cd 命令 进入该目录,然后运行 pcode * 命令进行伪代码转换。因为该命令会将 MATLAB 函数转换成伪代码,而在转换过程中该程序将自动翻译每一条语句,所以一旦发现有语法错误,将会停止翻译,给出错误信息。改正了该语法错误后,再重新执行 pcode 命令,直到没有错误为止。至少这样会保证目录下所有的程序不含有语法错误。
链接地址:https://www.31doc.com/p-2726310.html