MATLAB的使用

 

绪  论

Matlab是“Matrix Laboratory”的缩写，意为“矩阵实验室”，是当今美国很流行的科学计算软件．信息技术、计算机技术发展到今天，科学计算在各个领域得到了广泛的应用．在许多诸如控制论、时间序列分析、系统仿真、图像信号处理等方面产生了大量的矩阵及其相应的计算问题．自己去编写大量的繁复的计算程序，不仅会消耗大量的时间和精力，减缓工作进程，而且往往质量不高．美国Mathwork软件公司推出的Matlab软件就是为了给人们提供一个方便的数值计算平台而设计的．

Matlab是一个交互式的系统，它的基本运算单元是不需指定维数的矩阵，按照IEEE的数值计算标准（能正确处理无穷数Inf(Infinity)、无定义数NaN(not-a-number)及其运算）进行计算．系统提供了大量的矩阵及其它运算函数，可以方便地进行一些很复杂的计算，而且运算效率极高．Matlab命令和数学中的符号、公式非常接近，可读性强，容易掌握，还可利用它所提供的编程语言进行编程完成特定的工作．除基本部分外，Matlab还根据各专门领域中的特殊需要提供了许多可选的工具箱，如应用于自动控制领域的Control System工具箱和神经网络中Neural Network工具箱等．

第一节　Matlab的安装及使用

§1.1 Matlab的安装

版本Matlab 6.x for Windows．因为它使用方便，界面美观，我们选择它作为主要讲解版本．Matlab还有许多附加的部分，最常见的部分称为Simulink，是一个用作系统仿真的软件包，它可以让您定义各种部件，定义各自对某种信号的反应方式及与其它部件的连接方式．最后选择输入信号，系统会仿真运行整个模拟系统，并给出统计数据．Simulink有时是作为Matlab的一部分提供的，称为Matlab with Simulink版本．Matlab还有许多工具箱，它们是根据各个特殊领域的需要，用Matlab自身的语言编写的程序集，使用起来非常方便．您可以视工作性质和需要购买相应的工具箱．常见的工具箱有：

Signal Process	信号处理	System Identification	系统辨识
Optimization	优化	Neural Network	神经网络
Control System	自动控制	Spline	样条
Symbolic Math	符号代数	Image Process	图像处理
Nonlinear Control	非线性控制	Statistics	统计

§1.2 Matlab基本用法

从Windows中双击Matlab图标，会出现Matlab命令窗口（Command Window），在一段提示信息后，出现系统提示符“>>”．您可以在提示符后键入各种命令，通过上下箭头可以调出以前打入的命令，用滚动条可以查看以前的命令及其输出信息．

如果对一条命令的用法有疑问的话，可以用Help菜单中的相应选项查询有关信息，也可以用help命令在命令行上查询，您可以试一下help、help help和help eig（求特征值的函数）命令．

下面我们先从输入简单的矩阵开始掌握Matlab的功能．

§1.2.1输入简单的矩阵

输入一个小矩阵的最简单方法是用直接排列的形式．矩阵用方括号括起，元素之间用空格或逗号分隔，矩阵行与行之间用分号分开．例如输入：

A=[1  2  3 ;  4  5  6 ;  7  8  0]

系统会回答

A =

1      2     3

4      5     6

7      8     0

表示系统已经接收并处理了命令，在当前工作区内建立了矩阵A．

大的矩阵可以分行输入，用回车键代替分号，如：

A=[ 1     2     3

4    5     6

7    8     0]

结果和上式一样，也是

A =

1      2     3

4      5     6

7      8     0

§1.2.2矩阵元素

Matlab的矩阵元素可以是任何数值表达式．如：

x=[ -1.3   sqrt(3)   (1+2+3)*4/5]

结果：

x =

 -1.3000   1.7321    4.8000

在括号中加注下标，可取出单独的矩阵元素．如：

x(5)=abs(x(1))

结果

x =

 -1.3000   1.7321    4.8000    0     1.3000

注：结果中自动产生了向量的第5个元素，中间未定义的元素自动初始为零．

大的矩阵可把小的矩阵作为其元素来完成，如：

A=[A; [10  11  12]]

结果

A =

1      2     3

4      5     6

7      8     0

10    11   12

小矩阵可用“:”从大矩阵中抽取出来，如：

A=A(1:3,:);

即从A中取前三行和所有的列，重新组成原来的A. (详细介绍参见第二节的相关内容）

§1.2.3语句和变量

Matlab的表述语句、变量的类型说明由Matlab系统解释和判断．Matlab语句通常形式为：

变量=表达式

或者使用其简单形式为：

表达式

表达式由操作符或其它特殊字符、函数和变量名组成．表达式的结果为一个矩阵，显示在屏幕上，同时保存在变量中以留用．如果变量名和“=”省略，则具有ans名(意思指回答)的变量将自动建立．例如：

键入1900/81

结果为：

ans =

23.4568

需注意的问题有以下几点：

l 语句结束键入回车键，若语句的最后一个字符是分号，即“;”，则表明不输出当前命令的结果．

l 如果表达式很长，一行放不下，可以键入“ …”（三个点，但前面必须有个空格，目的是避免将形如“数2 …”理解为“数2.”与“..”的连接，从而导致错误），然后回车．

l 变量和函数名由字母加数字组成，但最多不能超过63个字符，否则系统只承认前63个字符．

l Matlab变量字母区分大小写，如A和a不是同一个变量，函数名一般使用小写字母，如inv(A)不能写成INV(A)，否则系统认为未定义函数．

§1.2.4 who和系统预定义变量

输入who命令可检查工作空间中建立的变量，键入：

who

系统输出为：

Your variables are:

A  ans  x

这里表明三个变量已由前面的例子产生了．

但列表中列出的并不是系统全部的变量，系统还有以下内部变量：

eps、pi、Inf、NaN

变量eps在决定诸如矩阵的奇异性时，可作为一个容许差，容许差的初值为1.0到1.0以后计算机所能表示的下一个最大浮点数，IEEE在各种计算机、工作站和个人计算机上使用这个算法．用户可将此值置为任何其它值(包括0值)．

变量pi是p.

Inf表示无穷大．如果您想计算1/0

S=1/0

结果会是

Warning：Divide by zero

S=Inf

具有IEEE规则的机器，被零除后，并不引出出错条件或终止程序的运行，而产生一个警告信息和一个特殊值在计算方程中列出来．

变量NaN表示它是个不定值．由Inf/Inf或0/0运算产生．

要了解当前变量的信息请键入whos，屏幕将显示：

  Name           Size               Bytes             Class

  A                4x3                96                  double array

  S                 1x1                8                  double array

  ans                     1x1                8                  double array

  x                 1x5                40                  double array

Grand total is 19 elementsusing 152 bytes

从size及bytes项目可以看出，每一个矩阵实元素需8个字节的内存．4×3的矩阵使用96个字节，全部变量的使用内存总数为152个字节．自由空间的大小决定了系统变量的多少，如计算机上有虚拟内存的话，其可定义的变量个数会大大增加．

§1.2.5数和算术表达式

Matlab中数的表示方法和一般的编程语言没有区别．如：

3                    -99                       0.0001

9.63972           1.6021E-20           6.02252e23

在计算中使用IEEE浮点算法其舍入误差是eps．浮点数表示范围是10^-308~10³⁰⁸．

数学运算符有：

+            加

-             减

*            乘

/             右除

\             左除

^            幂

这里1/4和4\1有相同的值都等于0.25（注意比较：1\4=4）．只有在矩阵的除法时左除和右除才有区别．

§1.2.6复数与矩阵

在Matlab中输入复数首先应该建立复数单位：

i=sqrt(-1)

及j=sqrt(-1)

之后复数可由下面语句给出：

Z=3+4i    (注意： 在4与i之间不要留有任何空间!)

输入复数矩阵有两个方便的方法，如：

A=[1  2; 3  4] + i*[5 6; 7  8]

和 A=[1+5i  2+6i;  3+7i  4+8i]

两式具有相等的结果．但当复数作为矩阵的元素输入时，不要留有任何空间，如1+5i，如在“+”号左右留有空格，就会被认为是两个分开的数．

不过实际使用复数时并没有这么麻烦，系统有一个名为startup.m的Matlab命令文件，建立复数单位的语句也放在其中．当Matlab启动时，此文件自动执行，i和j将自动建立．

§1.2.7输出格式

任何Matlab语句执行结果都可在屏幕上显示，同时赋给指定的变量，没有指定变量时赋给ans．数字显示格式可由format命令来控制（Windows系统下的Matlab系统的数字显示格式可以由Option菜单中的NumericalFormat菜单改变）．format仅影响矩阵的显示，不影响矩阵的计算与存贮．(Matlab以双精度执行所有的运算)

首先，如果矩阵元素是整数则矩阵显示就没有小数，如x=[-1 0 1]，结果为：

        x=

               -1    0     1

如果矩阵元素不是整数则输出形式有：（用命令：format 格式进行切换）

格式	中文解释	说明
format	短格式 （缺省格式）	Default. Same as SHORT
format short	短格式 （缺省格式）	Scaled fixed point format with 5 digits （只显示五位十进制数）
format long	长格式	Scaled fixed point format with 15 digits
format short e	短格式e方式	Floating point format with 5 digits
format long e	长格式e方式	Floating point format with 15 digits
format short g	短格式g方式	Best of fixed or floating point format with 5 digits
format long g	长格式g方式	Best of fixed or floating point format with 15 digits
format hex	16进制格式	Hexadecimal format
format +	+格式	The symbols +, - and blank are printed                         for positive, negative and zero elements.                        Imaginary parts are ignored
format bank	银行格式	Fixed format for dollars and cents
format rat	有理数格式	Approximation by ratio of small integers
format compact	压缩格式	Suppress extra line-feeds
format loose	自由格式	Puts the extra line-feeds back in

例如：

x=[4/3  1.2345e-6]

在不同的输出格式下的结果为：

短格式                   1.3333                                       0.0000

短格式e方式         1.3333e+000                               1.234e-006

长格式                   1.333333333333333                    0.000001234500000

长格式e方式         1.333333333333333e-000            1.23450000000000e-006

有理数格式            4/3                                             1/810045

16进制格式           3ff5555555555555                      3eb4b6231abfd271

+格式                    +                                               +

对于短格式，如果矩阵的最大元素比数999999999大，或者比数0.0001小，则在打印时，将加入一个普通的长度因数．如y=1.e20*x，意为x被10²⁰乘，结果为：

y=

         1.0e+020*

1.3333             0.0000

“+”格式是显示大矩阵的一种紧凑方法，“+”，“-”和空格显示正数、负数和零元素．

最后format compact命令压缩显示的矩阵，以允许更多的信息显示在屏幕上．

§1.2.8 Help求助命令和联机帮助

Help求助命令很有用，它对Matlab大部分命令提供了联机求助信息．您可以从Help菜单中选择相应的菜单，打开求助信息窗口查询某条命令，也可以直接用help命令．

键入help

得到help列表文件，键入“help 指定项目”，如：

键入help eig

则提供特征值函数的使用信息．

键入help [

显示如何使用方括号等．

键入help help

显示如何利用help本身的功能．

还有，键入lookfor  <关键字>：可以从m文件的help中查找有关的关键字．

§1.2.9 退出和存入工作空间

退出Matlab可键入quit或exit或选择相应的菜单．中止Matlab运行会引起工作空间中变量的丢失，因此在退出前，应键入save命令，保存工作空间中的变量以便以后使用．

键入 save

则将所有变量作为文件存入磁盘Matlab.mat中，下次Matlab启动时，

键入load

将变量从Matlab.mat中重新调出．

save和load后边可以跟文件名或指定的变量名，如仅有save时，则只能存入Matlab.mat中．如save temp命令，则将当前系统中的变量存入temp.mat中去，命令格式为：

savetemp x 仅仅存入x变量．

savetemp X Y Z 则存入X、Y、Z变量．

load  temp可重新从temp.mat文件中提出变量，load也可读ASCII数据文件．详细语法见联机帮助．

第二节　向量与矩阵运算

Matlab能处理数、向量和矩阵．但一个数事实上是一个1×1的矩阵，1个n维向量也不过是一个1×n或n×1的矩阵．从这个角度上来讲，Matlab处理的所有的数据都是矩阵．Matlab的矩阵处理能力是非常灵活、强大的．以下我们将从矩阵的产生、基本运算、矩阵函数等几个方面来说明．

§2.1向量及矩阵的生成

除了我们在上节介绍的直接列出矩阵元素的输入方法，矩阵还可以通过几种不同的方式输入到Matlab中．

§2.1.1 通过语句和函数产生

1. 向量的产生

除了直接列出向量元素（即所谓的“穷举法”）外，最常用的用来产生相同增量的向量的方法是利用“:”算符（即所谓的“描述法”）．在Matlab中，它是一个很重要的字符．如：

z=1:5

z =

     1    2     3    4     5

即产生一个1~5的单位增量是1的行向量，此为默认情况.

用“:”号也可以产生单位增量不等于1的行向量，语法是把增量放在起始量和结尾量的中间．如：

x=0:pi/4:pi

即产生一个由0~pi的行向量，单位增量是pi/4=3.1416/4=0.7854．

x =

 0     0.7854           1.5708           2.3562           3.1416

也可以产生单位增量为负数的行向量．如：

y=6:-1:1

y =

     6    5     4     3    2     1

2. 矩阵的产生

Matlab提供了一批产生矩阵的函数：

zeros	产生一个零矩阵	diag	产生一个对角矩阵
ones	生成全1矩阵	tril	取一个矩阵的下三角
eye	生成单位矩阵	triu	取一个矩阵的上三角
magic	生成魔术方阵	pascal	生成PASCAL矩阵

例如：

ones(3)

ans =

 1    1     1

 1    1     1

 1    1     1

eye(3)

ans =

 1    0     0

 0    1     0

 0    0     1

除了以上产生标准矩阵的函数外，Matlab还提供了产生随机（向量）矩阵的函数rand和randn，及产生均匀级数的函数linspace、产生对数级数的函数logspace和产生网格的函数meshgrid等等．详细使用请查阅随机文档．

“ : ”冒号可以用来产生简易的表格，为了产生纵向表格形式，首先用冒号“ : ”产生行向量，再进行转置，计算函数值的列，然后形成有二列的矩阵．例如命令：

x=(0.0:0.2:3.0)';

y=exp(-x).*sin(x);

[x y]

产生结果为：

ans =

0             0

0.2000      0.1627

0.4000      0.2610

0.6000      0.3099

0.8000      0.3223

1.0000      0.3096

1.2000      0.2807

1.4000      0.2430

1.6000      0.2018

1.8000      0.1610

2.0000      0.1231

2.2000      0.0896

2.4000      0.0613

2.6000      0.0383

2.8000      0.0204

3.0000      0.0070

§2.1.2 通过后缀为.m的命令文件产生

如有文件data.m，其中包括正文：

A=[  1 2 3

4 5 6

7 8 0]

则用data命令执行data.m，可以产生名为A的矩阵．

§2.2 矩阵操作

在Matlab中可以对矩阵进行任意操作，包括改变它的形式，取出子矩阵，扩充矩阵，旋转矩阵等．其中最重要的操作符为“:”， 它的作用是取出选定的行与列．

例如：

A(:,:) 代表A的所有元素；试比较A(:),将A按列的方向拉成长长的1列（向量）；

A(:,J) 代表A的第J列；

A(J:K) 代表 A(J), A(J+1), …, A(K)，如同A(:)的第J到第K个元素；

A(:,J:K) 代表A(:,J), A(:,J+1), …, A(:,K)，如此类推.

对矩阵可以进行各种各样的旋转、变形、扩充：

 

矩阵的转置用符号“ ' ”表示：

如A=[1  2  3;  4  5  6 ;  7  8  0]

那么：计算B=A'

B =

 1    4     7

 2    5     8

 3    6     0

符号“ ' ”为矩阵的转置，如果Z为复矩阵，则Z'为它的复数共轭转置，非共轭转置使用Z.' 或conj(Z')求得．

reshape改变矩阵的形状，这是什么意思呢？可举一个例子来说明．

A=[A;[10 11 12]]

A =

 1    2     3

 4    5     6

 7    8     0

 10 11   12

则 reshape(A,2,6)

ans =

 1    7     2     8    3     0

 4    10   5     11   6     12

可见，reshape 是将矩阵元素以列为单位进行重组，原来4×3的矩阵变为了2×6的矩阵．

第三节　矩阵的基本运算

§3.1 加和减

如矩阵A和B的维数相同，则A+B与A-B表示矩阵A与B的和与差．如果矩阵A和B的维数不匹配，Matlab会给出相应的错误提示信息．如：

A=                                B=

1     2     3                   1     4     7

4     5     6                   2     5     8

7     8     0                   3     6     0

C =A+B返回：

C =

    2     6    10

    6    10    14

   10    14     0

如果运算对象是个标量（即1×1矩阵），可和其它矩阵进行加减运算．例如：

x=    -1                  y=x-1=    -2

0                                 -1

2                                 1

§3.2矩阵乘法

Matlab中的矩阵乘法有通常意义上的矩阵乘法，也有Kronecker乘法，以下分别介绍．

§3.2.1 矩阵的普通乘法

矩阵乘法用“ * ”符号表示，当A矩阵列数与B矩阵的行数相等时，二者可以进行乘法运算，否则是错误的．计算方法和线性代数中所介绍的完全相同．

如：A=[1  2 ; 3  4]; B=[5  6 ; 7  8];  C=A*B，

结果为

C= × = =

即Matlab返回：

C =

    19   22

    43   50

如果A或B是标量，则A*B返回标量A（或B）乘上矩阵B（或A）的每一个元素所得的矩阵．

§3.3 矩阵除法

在Matlab中有两种矩阵除法符号：“＼”即左除和“／”即右除．如果A矩阵是非奇异方阵，则A\B是A的逆矩阵乘B，即inv(A)*B；而B/A是B乘A的逆矩阵，即B*inv(A)．具体计算时可不用逆矩阵而直接计算．

通常：

x=A\B就是A*x=B的解；

x=B/A就是x*A=B的解.

当B与A矩阵行数相等可进行左除．如果A是方阵，用高斯消元法分解因数．解方程：A*x(:, j)=B(:, j)，式中的(:, j)表示B矩阵的第j列，返回的结果x具有与B矩阵相同的阶数，如果A是奇异矩阵将给出警告信息．

如果A矩阵不是方阵，可由以列为基准的Householder正交分解法分解，这种分解法可以解决在最小二乘法中的欠定方程或超定方程，结果是m×n的x矩阵．m是A矩阵的列数，n是B矩阵的列数．每个矩阵的列向量最多有k个非零元素，k 是A的有效秩．

右除B/A可由B/A=(A'\B')'左除来实现．

§3.4矩阵乘方

A^P意思是A的P次方．如果A是一个方阵，P是一个大于1的整数，则A^P表示A的P次幂，即A自乘P次．如果P不是整数，计算涉及到特征值和特征向量的问题，如已经求得：[V,D]=eig(A)，则：

A^P=V*D.^P/V（注：这里的.^表示数组乘方，或点乘方，参见后面的有关介绍）

如果B是方阵， a是标量，a^B就是一个按特征值与特征向量的升幂排列的B次方程阵． 如果a和B都是矩阵，则a^B是错误的．

§3.5 矩阵的超越函数

在Matlab中解释exp(A)和sqrt(A)时曾涉及到级数运算，此运算定义在A的单个元素上． Matlab可以计算矩阵的超越函数，如矩阵指数、矩阵对数等．

一个超越函数可以作为矩阵函数来解释，例如将“m”加在函数名的后边而成expm(A)和sqrtm(A)，当Matlab运行时，有下列三种函数定义：

expm           矩阵指数

logm        矩阵对数

sqrtm           矩阵开方

所列各项可以加在多种m文件中或使用funm．请见应用库中sqrtm.m，1ogm.m，funm.m文件和命令手册．

§3.6数组运算

数组运算由线性代数的矩阵运算符“*”、“/”、“\”、“^”前加一点来表示，即为“.*”、“./”、“.\”、“.^”．注意没有“.+”、“.-”运算．

§3.6.1数组的加和减

对于数组的加和减运算与矩阵运算相同，所以“+”、“-”既可被矩阵接受又可被数组接受．

§3.6.2数组的乘和除

数组的乘用符号.*表示，如果A与B矩阵具有相同阶数，则A.*B表示A和B单个元素之间的对应相乘．例如x=[1  2  3]; y=[ 4  5  6];

计算z=x.*y

结果z=4  10  18

数组的左除（.\）与数组的右除（./），由读者自行举例加以体会．

§3.6.3 数组乘方

数组乘方用符号.^表示．

例如：键入：

x=[ 1  2  3]

y=[ 4  5  6]

则z=x.^y=[1^4  2^5  3^6]=[1  32  729]

(1) 如指数是个标量，例如x.^2，x同上，则：

z=x.^2=[1^2  2^2  3^2]=[ 1  4  9]

(2) 如底是标量，例如2 .^[x y] ，x、y同上，则：

z=2 .^[x y]=[2^1  2^2  2^3  2^4 2^5  2^6]=[2  4  8  16  32  64]

从此例可以看出Matlab算法的微妙特性，虽然看上去与其它乘方没什么不同，但在2和“．”之间的空格很重要，如果不这样做，解释程序会把“．”看成是2的小数点． Matlab看到符号“^”时，就会当做矩阵的幂来运算，这种情况就会出错，因为指数矩阵不是方阵．

§3.7 矩阵函数

Matlab的数学能力大部分是从它的矩阵函数派生出来的，其中一部分装入Matlab本身处理中，它从外部的Matlab建立的M文件库中得到，还有一些由个别的用户为其自己的特殊的用途加进去的．其它功能函数在求助程序或命令手册中都可找到．手册中备有为Matlab提供数学基础的LINPACK和EISPACK软件包，提供了下面四种情况的分解函数或变换函数：

(1)三角分解；(2)正交变换；(3) 特征值变换；(4)奇异值分解．

§3.7.1三角分解

最基本的分解为“LU”分解，矩阵分解为两个基本三角矩阵形成的方阵，三角矩阵有上三角矩阵和下三角矩阵．计算算法用高斯变量消去法．

从lu函数中可以得到分解出的上三角与下三角矩阵，函数inv得到矩阵的逆矩阵，det得到矩阵的行列式．解线性方程组的结果由方阵的“\”和“/”矩阵除法来得到．

例如：

A=[ 1      2     3

4      5     6

7      8     0]

LU分解，用Matlab的多重赋值语句

[L,U]=lu(A)

得出

L =

0.1429	1.0000	0
0.5714	0.5000	1.0000
1.0000	0	0

U =

7.0000	8.0000	0
0	0.8571	3.0000
0	0	4.5000

注：L是下三角矩阵的置换，U是上三角矩阵的正交变换，分解作如下运算，检测计算结果只需计算L*U即可．

求逆由下式给出： x=inv(A)

x =

-1.7778	0.8889	-0.1111
1.5556	-0.7778	0.2222
-0.1111	0.2222	-0.1111

从LU分解得到的行列式的值是精确的，d=det(U)*det(L)的值可由下式给出：

d=det(A)

d =

27

直接由三角分解计算行列式：d=det(L)*det(U)

d =

27.0000

为什么两种d的显示格式不一样呢? 当Matlab做det(A)运算时，所有A的元素都是整数，所以结果为整数．但是用LU分解计算d时，L、U的元素是实数，所以Matlab产生的d也是实数．

例如：线性联立方程取 b=[ 1

                                           3

                                           5]

解Ax=b方程，用Matlab矩阵除得到

x=A\b

结果x=

0.3333

0.3333

0.0000

由于A=L*U，所以x也可以有以下两个式子计算：y=L\b，x=U\y．得到相同的x值，中间值y为：

y =

5.0000

0.2857

0.0000

Matlab中与此相关的函数还有rcond、chol和rref．其基本算法与LU分解密切相关．chol函数对正定矩阵进行Cholesky分解，产生一个上三角矩阵，以使R'*R=X．rref用具有部分主元的高斯－约当消去法产生矩阵A的化简梯形形式．虽然计算量很少，但它是很有趣的理论线性代数．为了教学的要求，也包括在Matlab中．

§3.7.2正交变换

“QR”分解用于矩阵的正交－三角分解．它将矩阵分解为实正交矩阵或复酉矩阵与上三角矩阵的积，对方阵和长方阵都很有用．

例如A=[  1     2     3

4     5     6

7     8     9

10    11    12]

是一个降秩矩阵，中间列是其它二列的平均，我们对它进行QR分解：

[Q,R]=qr(A)

Q =

-0.0776	-0.8331	0.5444	0.0605
-0.3105	-0.4512	-0.7709	0.3251
-0.5433	-0.0694	-0.0913	-0.8317
-0.7762	0.3124	0.3178	0.4461

R =

-12.8841	-14.5916	-16.2992
0	-1.0413	-2.0826
0	0	0.0000
0	0	0

可以验证Q*R就是原来的A矩阵．由R的下三角都给出0，并且R(3,3)=0.0000，说明矩阵R与原来矩阵A都不是满秩的．

下面尝试利用QR分解来求超定和降秩的线性方程组的解．

例如：

b=[ 1

3

5

7]

讨论线性方程组Ax=b，我们可以知道方程组是超定的，采用最小二乘法的最好结果是计算x=A\b．

结果为：

Warning: Rank deficient, rank = 2 tol =   1.4594e-014

x =

    0.5000

         0

    0.1667

我们得到了缺秩的警告．用QR分解法计算此方程组分二个步骤：

y=Q'*b

x=R\y

求出的y值为

y =

-9.1586

-0.3471

0.0000

0.0000

x的结果为

Warning: Rank deficient, rank = 2 tol =   1.4594e-014

x =

    0.5000

         0

    0.1667

用A*x来验证计算结果，我们会发现在允许的误差范围内结果等于b．这告诉我们虽然联立方程Ax=b是超定和降秩的，但两种求解方法的结果是一致的．显然x向量的解有无穷多个，而“QR”分解仅仅找出了其中之一．

§3.7.3奇异值分解

在Matlab中三重赋值语句

[U,S,V]=svd(A)

在奇异值分解中产生三个因数：

A=U*S*V '

U矩阵和V矩阵是正交矩阵，S矩阵是对角矩阵，svd(A)函数恰好返回S的对角元素，而且就是A的奇异值（其定义为：矩阵A'*A的特征值的算术平方根）．注意到A矩阵可以不是方的矩阵．

奇异值分解可被其它几种函数使用，包括广义逆矩阵pinv(A)、秩rank(A)、欧几里德矩阵范数norm(A,2)和条件数cond(A)．

§3.7.4 特征值分解

如果A是n×n矩阵，若l满足Ax=lx，则称l为A的特征值，x为相应的特征向量．

函数eig(A)返回特征值列向量，如果A是实对称的，特征值为实数．特征值也可能为复数，例如：

A=[  0     1

-1     0]

eig(A)

产生结果

ans =

0 + 1.0000i

0 - 1.0000i

如果还要求求出特征向量，则可以用eig(A)函数的第二个返回值得到：

[x,D]=eig(A)

D的对角元素是特征值．x的列是相应的特征向量，以使A*x=x*D．

计算特征值的中间结果有两种形式：

Hessenberg形式为hess(A)，Schur形式为schur(A)．

schur形式用来计算矩阵的超越函数，诸如sqrtm(A)和logm(A)．

如果A和B是方阵，函数eig(A,B)返回一个包含一般特征值的向量来解方程

Ax=lBx

双赋值获得特征向量

[X,D]=eig(A,B)

产生特征值为对角矩阵D．满秩矩阵X的列相应于特征向量，使A*X=B*X*D，中间结果由qz(A,B)提供．

§3.7.5秩

Matlab计算矩阵A的秩的函数为rank(A)，与秩的计算相关的函数还有：rref(A)、orth(A)、null(A)和广义逆矩阵pinv(A)等．

利用rref(A)，A的秩为非0行的个数．rref方法是几个定秩算法中最快的一个，但结果上并不可靠和完善．pinv(A)是基于奇异值的算法．该算法消耗时间多，但比较可靠．其它函数的详细用法可利用Help求助．

第四节　Matlab中的图形

§4.1　二维作图

绘图命令plot绘制x-y坐标图；polor命令绘制极坐标图．

§4.1.1 基本形式

如果y是一个向量，那么plot(y)绘制一个y中元素的线性图．假设我们希望画出

y=[0., 0.48,  0.84,  1., 0.91,  6.14 ]

则用命令：plot(y)

它相当于命令：plot(x, y)，其中x=[1,2,…,n]或x=[1;2;…;n]，即向量y的下标编号, n为向量y的长度

Matlab会产生一个图形窗口，显示如下图形，请注意：坐标x和y 是由计算机自动绘出的．

图4.1.1.1plot([0.,0.48,0.84,1.,0.91,6.14])

 

上面的图形没有加上x轴和y轴的标注，也没有标题．用xlabel，ylabel，title命令可以加上．

如果x，y是同样长度的向量，plot(x,y)命令可画出相应的x元素与y元素的x-y坐标图．例：

x=0:0.05:4*pi;  y=sin(x);  plot(x,y)

grid on, title(' y=sin( x ) 曲线图' )

xlabel(' x = 0 : 0.05 : 4Pi ')

结果见下图．

图4.1.1.2 y=sin(x)的图形

 

title	图形标题
xlabel	x坐标轴标注
ylabel	y坐标轴标注
text	标注数据点
grid	给图形加上网格
hold	保持图形窗口的图形

表4.1.1.1 Matlab图形命令

§4.1.2 多重线

在一个单线图上，绘制多重线有三种办法.

第一种方法是利用plot的多变量方式绘制：

plot(x1,y1,x2,y2,...,xn,yn)

x1,y1,x2,y2,...,xn,yn是成对的向量，每一对x, y在图上产生如上方式的单线．多变量方式绘图是允许不同长度的向量显示在同一图形上．

第二种方法也是利用plot绘制，但加上hold on/off命令的配合：

plot(x1,y1)

hold on

plot(x2,y2)

hold off

第三种方法还是利用plot绘制，但代入矩阵：

如果plot用于两个变量plot(x,y)，并且x，y是矩阵，则有以下情况：

（1）如果y是矩阵，x是向量，plot(x,y)用不同的画线形式绘出y的行或列及相应的x向量，y的行或列的方向与x向量元素的值选择是相同的．

（2）如果x是矩阵，y是向量，则除了x向量的线族及相应的y向量外，以上的规则也适用．

（3）如果x，y是同样大小的矩阵，plot(x,y)绘制x的列及y相应的列．

还有其它一些情况，请参见Matlab的帮助系统．

§4.1.3 线型和颜色的控制

如果不指定划线方式和颜色，Matlab会自动为您选择点的表示方式及颜色．您也可以用不同的符号指定不同的曲线绘制方式．例如：

plot(x,y,'*')                    用'*'作为点绘制的图形．

plot(x1,y1,':',x2,y2,'+')    用':'画第一条线，用'+'画第二条线．

线型、点标记和颜色的取值有以下几种：

线型		点标记		颜色
-	实线	.	点	y	黄
:	虚线	o	小圆圈	m	棕色
-.	点划线	x	叉子符	c	青色
--	间断线	+	加号	r	红色
		*	星号	g	绿色
		s	方格	b	蓝色
		d	菱形	w	白色
		^	朝上三角	k	黑色
		v	朝下三角
		>	朝右三角
		<	朝左三角
		p	五角星
		h	六角星

表4.1.3.1线型和颜色控制符

如果你的计算机系统不支持彩色显示，Matlab将把颜色符号解释为线型符号，用不同的线型表示不同的颜色．颜色与线型也可以一起给出，即同时指定曲线的颜色和线型．

例如：     t=-3.14:0.2:3.14;

x=sin(t);  y=cos(t);

plot(t,x, '+r',t,y, '-b')

图4.1.3.1不同线型、颜色的sin,cos图形

 

§4.1.4对数图、极坐标图及条形图

 

 

例6：极坐标绘图

t=0:0.01:2*pi;

polar(t,sin(6*t))   %% (图4.1.4.3)

    

             图4.1.4.3 极坐标绘图         图4.1.4.4正态分布的统计直方图与其正态分布拟合曲线

 

例7：正态分布图

我们可以用命令normrnd生成符合正态分布的随机数．

normrnd(u,v,m,n)

其中，u表示生成随机数的期望，v代表随机数的方差．

运行：

a=normrnd(10,2,10000,1);

histfit(a)   %% (图4.1.4.4)

我们可以得到正态分布的统计直方图与其正态分布拟合曲线．

 

例8：比较正态分布（图4.1.4.5（1））与平均分布（图4.1.4.5（2））的分布图：

yn=randn(30000,1);    %% 正态分布

x=min(yn): 0.2 : max(yn);

subplot(121)

hist(yn,x)

yu=rand(30000,1);     %% 平均分布

subplot(122)

hist(yu,25)

4.1.4.5(1)                          4.1.4.5(2)

 

图4.1.4.5 正态分布与平均分布的分布图

§4.1.5 子图

在绘图过程中，经常要把几个图形在同一个图形窗口中表现出来，而不是简单地叠加（例如上面的例8）．这就用到函数subplot．其调用格式如下：

subplot(m,n,p)

subplot函数把一个图形窗口分割成m×n个子区域，用户可以通过参数p调用个各子绘图区域进行操作．子绘图区域的编号为按行从左至右编号．

 

例9：绘制子图

x=0:0.1*pi:2*pi;

subplot(2,2,1)

plot(x,sin(x),'-*');

title('sin(x)');

subplot(2,2,2)

plot(x,cos(x),'--o');

title('cos(x)');

subplot(2,2,3)

plot(x,sin(2*x),'-.*');

title('sin(2x)');

subplot(2,2,4);

plot(x,cos(3*x),':d')

title('cos(3x)')

得到图形如下：

图4.1.5.1子图

§4.1.6 填充图

利用二维绘图函数patch，我们可绘制填充图．绘制填充图的另一个函数为fill．

下面的例子绘出了函数humps（一个Matlab演示函数）在指定区域内的函数图形．

例10：用函数patch绘制填充图

fplot('humps',[0,2],'b')

hold on

patch([0.5 0.5:0.02:1 1],[0humps(0.5:0.02:1) 0],'r');

hold off

title('A region under an interesting function.')

grid

图4.1.6.1填充图

 

我们还可以用函数fill来绘制类似的填充图．

例11：用函数fill绘制填充图

x=0:pi/60:2*pi;

y=sin(x);

x1=0:pi/60:1;

y1=sin(x1);

plot(x,y,'r');

hold on

fill([x1 1],[y10],'g')

图4.1.6.2填充图

§4.2 三维作图

§4.2.1 mesh(Z)语句

mesh(Z)语句可以给出矩阵Z元素的三维消隐图，网络表面由Z坐标点定义，与前面叙述的x-y平面的线格相同，图形由邻近的点连接而成．它可用来显示用其它方式难以输出的包含大量数据的大型矩阵，也可用来绘制Z变量函数．

显示两变量的函数Z=f(x,y)，第一步需产生特定的行和列的x-y矩阵．然后计算函数在各网格点上的值．最后用mesh函数输出．

下面我们绘制sin(r)/r函数的图形．建立图形用以下方法：

x=-8:.5:8;

y=x';

x=ones(size(y))*x;

y=y*ones(size(y))';

R=sqrt(x.^2+y.^2)+eps;

z=sin(R)./R;

mesh(z)       %%  试运行mesh(x,y,z)，看看与mesh(z)有什么不同之处？

各语句的意义是：首先建立行向量x，列向量y；然后按向量的长度建立1-矩阵；用向量乘以产生的1-矩阵，生成网格矩阵，它们的值对应于x-y坐标平面；接下来计算各网格点的半径；最后计算函数值矩阵Z．用mesh函数即可以得到图形．

图4.2.1三维消隐图

第一条语句x的赋值为定义域，在其上估计函数；第三条语句建立一个重复行的x矩阵，第四条语句产生y的响应，第五条语句产生矩阵R（其元素为各网格点到原点的距离）．用mesh方法结果如上．

另外，上述命令系列中的前4行可用以下一条命令替代：

[x, y]=meshgrid(-8:0.5:8)

§4.2.2与mesh相关的几个函数

(1) meshc与函数mesh的调用方式相同，只是该函数在mesh的基础上又增加了绘制相应等高线的功能．下面来看一个meshc的例子：

[x,y]=meshgrid([-4:.5:4]);

z=sqrt(x.^2+y.^2);

meshc(z)       %% 试运行 meshc(x,y,z)，看看与meshc(z)有什么不同之处？

我们可以得到图形：

图4.2.2.1  meshc 图

 

地面上的圆圈就是上面图形的等高线．

(2) 函数meshz与mesh的调用方式也相同，不同的是该函数在mesh函数的作用之上增加了屏蔽作用，即增加了边界面屏蔽．例如：

[x,y]=meshgrid([-4:.5:4]);

z=sqrt(x.^2+y.^2);

meshz(z)       %% 试运行 meshz(x,y,z)，看看与meshz(z)有什么不同之处？

我们得到图形：

图4.2.2.2  meshz 图

 

§4.2.3 其它的几个三维绘图函数

(1) 在Matlab中有一个专门绘制圆球体的函数sphere，其调用格式如下：

[x,y,z]=sphere(n) 

此函数生成三个(n+1)×(n+1)阶的矩阵，再利用函数surf(x,y,z)可生成单位球面．

[x,y,z]=sphere   此形式使用了默认值n=20

sphere(n)   只绘制球面图，不返回值．

运行下面程序：

sphere(30);

axis square;

我们得到球体图形：

图4.2.3.1  球面图

 

若只输入sphere画图，则是默认了n=20的情况．

 

(2) surf函数也是Matlab中常用的三维绘图函数．其调用格式如下：

surf(x,y,z,c)

输入参数的设置与mesh相同，不同的是mesh函数绘制的是一网格图，而surf绘制的是着色的三维表面．Matlab语言对表面进行着色的方法是，在得到相应网格后，对每一网格依据该网格所代表的节点的色值（由变量c控制），来定义这一网格的颜色．若不输入c，则默认为c=z．

我们看下面的例子：

%绘制地球表面的气温分布示意图．

[a,b,c]=sphere(40);

t=abs(c);    %求绝对值

surf(a,b,c,t);

axis equal

colormap('hot')

我们可以得到图形如下：

图4.2.3.2  等温线示意图

 

 

§4.2.4图形的控制与修饰

(1) 坐标轴的控制函数axis，调用格式如下：

axis([xmin,xmax,ymin,ymax,zmin,zmax])

用此命令可以控制坐标轴的范围．

与axis相关的几条常用命令还有：

axis auto                自动模式，使得图形的坐标范围满足图中一切图元素

axis equal                     严格控制各坐标的分度使其相等

axis square             使绘图区为正方形

axis on                  恢复对坐标轴的一切设置

axis off                 取消对坐标轴的一切设置

axis manual           以当前的坐标限制图形的绘制

 

（2）grid on  在图形中绘制坐标网格．

grid off  取消坐标网格．

 

（3）xlabel,ylabel, zlabel分别为x轴, y轴, z轴添加标注．title为图形添加标题．

以上函数的调用格式大同小异，我们以xlabel为例进行介绍：

xlabel('标注文本'，'属性1'，'属性值1'，'属性2'，'属性值2'，…)

这里的属性是标注文本的属性，包括字体大小、字体名、字体粗细等．

 

例如：

[x, y]=meshgrid(-4:.2:4);

R=sqrt(x.^2+y.^2);

z=-cos(R);

mesh(x,y,z)

xlabel('x\in[-4,4]','fontweight','bold');

ylabel('y\in[-4,4]','fontweight','bold');

zlabel('z=-cos(sqrt(x^2+y^2))','fontweight','bold');

title('旋转曲面','fontsize',15,'fontweight','bold','fontname','隶书');

 

 

图4.2.4.1 添加标注

 

以上各种绘图方法的详细用法，请看联机信息．

§4.3 统计回归图

对平面上n个点：

在平面直线族{ 为实数}中寻求一条直线 ，使得散点到与散点相对应的在直线上的点之间的纵坐标的误差的平方和最小，用微积分的方法可得：

 

所求得的这条直线： 称为回归直线．

例：已知如下点列，求其回归直线，并计算最小误差平方和．

x	0.1	0.11	.12	.13	.14	.15	.16	.17	.18	.2	.21	.23
y	42	43.5	45	45.5	45	47.5	49	53	50	55	55	60

参考的程序如下：

x=[0.10.11 .12 .13 .14 .15 .16 .17 .18 .2.21 .23];

y=[4243.5 45 45.5 45 47.5 49 53 50 55 55 60];

n=length(x);

xb=mean(x);  

yb=mean(y);

x2b=sum(x.^2)/n;

xyb=x*y'/n;

b=(xb*yb-xyb)/(xb^2-x2b);

a=yb-b*xb;

y1=a+b.*x;

plot(x,y,'*',x,y1);

serror=sum((y-y1).^2)

 

 

图4.3.1回归直线

第五节　Matlab编程

§5.1关系运算

§5.1.1 比较运算

比较两个同阶矩阵有下面六种相关操作符：

相关操作符
<	小于
<=	小于等于
>	大于
>=	大于等于
==	等于
~=	不等于

表5.1.1.1相关操作符

比较两个元素的大小，结果是“1”表明为真，结果是“0”表明为假．

例如2+2~=4

结果是“0”，表明为假．

例如一个6阶魔术方阵，矩阵元素计算满足各种条件：

A=magic(6)

ans =

    35    1     6    26   19    24

     3   32     7    21   23    25

    31    9     2    22   27    20

     8   28    33    17   10    15

    30    5    34    12   14    16

     4   36    29    13   18    11

阶数为n的魔术方阵，即n×n矩阵，是由1～n²的整数组成(n=6)．仔细观察这个矩阵，我们会发现任何行和、任何列和都相等．另外，每个3×3子行列式的对角线元素和，都可被3整除．为了显示这一特性，键入：

p=(rem(A,3)==0)

p =

     0    0     1     0    0     1

     1    0     0     1    0     0

     0    1     0     0    1     0

     0    0     1     0    0     1

     1    0     0     1    0     0

     0    1     0     0    1     0

为了再仔细地观察这个模式，可以用format+格式画出矩阵的压缩格式．此格式用“+”代表正元素，“-”代表负元素，空格代表0．

format+

p =

  +  +

+  + 

 +  +

  +  +

+  + 

 +  +

find 函数在关系运算中很有用，它可以在0-1矩阵中找非零元素的下标．

若y是一个向量，例如：y=[1   3   2   4   3.5   2.9]，则find(y<3.0)，将指出y的分量在哪些位置上小于3.0．

ans=     1     3    6

即：向量y的第1、3、6位置上的元素小于3.0．

当输入x==NaN时结果为NaN，因为按照IEEE算法规定任何具有NaN的操作，结果都是NaN．调试NaN很有用，例如测试x，输入isnan(x)函数，如果x元素是不定值则得1，否则得0．isfinite(x)更有用，如-¥<x<¥时则得1．

§5.1.2逻辑运算

&	与
|	或
~	非

表5.1.2.1逻辑运算符

 

“＆”和“|”操作符可比较两个标量或两个同阶矩阵．对于矩阵来说必须符合规则，如果A和B都是0-1矩阵，则A＆B或A|B也都是0-1矩阵，这个0-1矩阵的元素是A和B对应元素之间逻辑运算的结果，逻辑操作符认定任何非零元素都为真，给出“1”，任何零元素都为假，给出“0”．

非(或逻辑非)是一元操作符，即~A：当A是非零时结果为“0”；当A为“0”时，结果为“1”．因此下列两种表示：

p |(~p)             结果为1．

p ＆ (~p)          结果为0．

any和all函数在连接操作时很有用，设x是0-1向量，如果x中任意有一元素非零时，any(x)返回“1”，否则返回“0”；all(x)函数当x的所有元素非零时，返回“1”，否则也返回“0”．这些函数在if语句中经常被用到．如：

ifall(A<5)

do something

end

§5.2控制流

Matlab与其它计算机语言一样，也有控制流语句．控制流语句可使原本简单地在命令行中运行的一系列命令或函数，组合成为一个整体——程序，从而提高工作效率．

§5.2.1 for 循环

Matlab与其它计算机语言一样有do或for循环，完成一个语句或一组语句在一定时间内反复运行的功能．例如：

for  i = 1:n ,  x( i )=0,  end

x的第一个元素赋0值，如果n＜1，结构上合法，但内部语句不运行，如果x不存在或比n元素小，额外的空间将会自动分配．

多重循环写成锯齿形是为了增加可读性．例如：

m=9;n=9;

for i= 1:m

for j=1:n

               A( i, j ) = 1/( i + j - 1);

end

end

A

程序的说明：

（1）事实上，上述程序给出了Hilbert矩阵的构造过程，可参见函数hilb(n)．

（2）语句内部使用分号，表示计算过程不输出中间结果．

（3）循环后的A命令表示显示矩阵A的结果．

（4）每个for语句必须以end语句结束，否则是错误的．

for循环的通用形式为：

for v=expression

statements

end

其中expression表达式是一个矩阵，因为Matlab中都是矩阵，矩阵的列被一个接一个的赋值到变量v，然后statements语句运行．

通常expression是一些m:n或m:k:n仅有一行的矩阵，并且它的列是个简单的标量．但如注意到expression可以为矩阵，即v可以为向量，对某些问题的处理将大大简化．

§5.2.2  while 循环

Matlab中的while 循环语句为一个语句或一组语句在一个逻辑条件的控制下重复未知的次数．

它的一般形式为：

while  expression

statements

end

当expression的所有运算为非零值时，statements语句组将被执行．如果判断条件是向量或矩阵的话，可能需要all或any函数作为判断条件．

例如计算expm(A)，在A并不是太大时，直接计算expm(A)是可行的．

expm(A)=I+A+A^2/2!+A^3/3!+…   注意：这里的I表示单位矩阵．

程序为：

E = 0*A; F = E + eye(size(E)); N = 1;

while norm(F,1) > 0,

F = A*F/N;

E = E + F;

N = N + 1;

end

§5.2.3 if和break语句

下面介绍if语句的二个例子．

(1) 一个计算如何被分成三个部分，用符号校验：

if n<0

A=negative(n)

elseif mod(n,2)==0

A=even(n)

else

A=odd(n)

end

其中的三个函数negative(n)、even(n)、odd(n)是自编的输出函数．参见下面的函数文件．

 

(2) 这个例子涉及数论中一个很有趣的问题，取任何的正整数，如果是偶数，用2除；如果是奇数，用3乘，并加上1，反复这个过程，直到你的整数成为1．这个极有趣不可解的问题是：有使这个过程不中止的整数吗?

%classic "3n+1"problem from number theory

while 1

n=input('Enter n, negative quits:');

if n<=0 break,end

        while n>1

               if rem(n,2) == 0  %% 是连续的2个等号

                      n=n/2

               else

                      n=3*n+1

               end;

        end

end

这个过程能永远进行吗?

程序的说明：

（1）本程序用到了if语句与while语句，过程比较复杂；

（2）使用input函数，可使程序在执行过程中，从键盘输入一个数（矩阵）；

（3）break语句提供了程序跳出死循环的途径．

§5.3 M文件、命令文件及函数文件

§5.3.1 M文件

Matlab通常使用命令驱动方式，当单行命令输入时，Matlab立即处理并显示结果，同时将运行说明或命令存入文件．

Matlab语句的磁盘文件称作M文件，因为这些文件名的未尾是.M形式，例如一个文件名为bessel.m，提供bessel函数语句．

一个M文件包含一系列的Matlab语句，一个M文件可以循环地调用它自己．

M文件有两种类型：

第一类型的M文件称为命令文件，它是一系列命令、语句的简单组合．

第二类型的M文件称为函数文件，它提供了Matlab的外部函数．用户为解决一个特定问题而编写的大量的外部函数可放在Matlab工具箱中，这样的一组外部函数形成一个专用的软件包．

这两种形式的M文件，无论是命令文件，还是函数文件，都是普通的ASCII文本文件，可选择编辑或字处理文件来建立．

§5.3.2命令文件

当一个命令文件被调用时，Matlab运行文件中出现的命令而不是交互地等待键盘输入，命令文件的语句在工作空间中运算全局数据，对于进行分析解决问题及做设计中所需的一长串繁杂的命令和解释是很有用的．

例如：一个自编的命令文件fibo.m，用于计算Fibonnaci数列

% An M-file to calculate Fibonnaci numbers

f=[1, 1 ]; i = 1;

while f(i)+f(i+1)<1000

f(i+2)=f(i)+f(i+1);

i=i+1;

end

plot(f)

在Matlab命令窗口中键入fibo命令，并回车执行，将计算出所有小于1000的Fibonnaci数，并绘出图形．

要注意的是：文件执行后，f和i变量仍然留在工作空间．

§5.3.3函数文件

如果M文件的第一行包含function，这个文件就是函数文件，它与命令文件不同，所定义变量和运算都在文件内部，而不在工作空间．函数被调用完毕后，所定义变量和运算将全部释放．函数文件对扩展Matlab函数非常有用．

例如：一个自编的函数文件mean.m，用于求向量的（或矩阵按列的）平均值

function  y=mean(x)

% MEAN Average or mean value，For Vectors,

% MEAN (x) returns the mean value

% For matrix MEAN (x) is a row vector

% containing the mean value of each column

[m,n]=size(x);

if m==1

m=n;

end

y=sum(x)/m;

磁盘文件中定义的新函数称为mean函数，它与Matlab函数一样使用，例如z为从1到99的实数向量：

z=1:99;

计算均值：mean(z)

ans=

50

mean.m程序的说明：

（1）第一行的内容：函数名，输入变量，输出变量，没有这行这个文件就是命令文件，而不是函数文件．

（2）％：表明％右边的行是说明性的内容注释．前一小部分行来确定M文件的注释，并在键入helpmean后显示出来．显示内容为连续的若干个％右边的文字．

（3）变量m,n和y是mean的局部变量，在mean运行结束后，它们将不在工作空间z中存在．如果在调用函数之前有同名变量，先前存在的变量及其当前值将不会改变．

 

再例如：一个计算标准差的函数文件stat.m

function [mean,stdev]=stat(x)

[m,n]=size(x);

if m==1

m=n

end

mean=sum(x)/m;

stdev=sqrt(sum(x.^2)/m-mean.^2);

stat表明返回多输出变量是可能的．

 

又如：使用多输入变量计算矩阵秩函数

function r=rank(x,tol)

% rank of a matrix

s=svd(x);

if(nargin==1)

tol=max(size(x))*s(1)*eps;

end

r=sum(s>tol);

这个变量说明利用永久变量nargin确定输入变量的个数，变量nargout虽然这里没有使用，但它包含有输出变量的个数．

一些有用的说明：

当M函数文件第一次在Matlab运行时，它被编译并放入内存，以后使用时不用重新编译即可得到．

what命令：显示磁盘当前目录中的M文件，

dir命令：列出所有文件．

一般而言，输入一个名字到Matlab，例如键入whoopie命令，Matlab用以下步骤解释：

(1) 看whoopie是否为变量．

(2) 检验whoopie是否为在线函数．

(3) 检验whoopie文件的当前目录．

(4) 将whoopie看成Matlab的PATH中的一个文件，在Matlab PATH目录中搜索．

如果whoopie存在，Matlab首先将其作为变量而不是作为函数．

§5.4 字符串、输入及输出

§5.4.1 echo、input、pause、keyboard

一般来说，当一个M文件运行时，文件的命令不在屏幕上显示，而echo命令则使M文件运行时，命令在屏幕上显示，这对于调试、演示相当有用．

input功能：输入Input('How many apples')给用户一个提示串，等待，然后显示用户通过键盘输入的大量表达式．可以用input命令建立驱动M文件的菜单．

与input功能相同，但功能更强的keyboard命令将计算机作为一个命令文件来调用，放入M文件中，此特性对调试或正在运行期间修改变量很有用．

pause命令：使用户暂停运行一个程序，当再按任一键时恢复执行，pause(n)等待n秒钟后再继续执行．

§5.4.2串和宏串

字符串用单个引号输入到Matlab中，例如：

s='Hello'

结果显示为：

s =

Hello

字符存在向量中，每个元素就是一个字符，如：

size(s)

ans =

1       5

表明S为一个1×5的矩阵，有五个元素．字符以ASCII值存入，abs函数或double函数将显示以下值(即Hello的ASCII值)

abs(s)

ans =

   72   101  108   108   111

getstr函数，使向量作为字符显示，而不显示ASCII值．

disp可在变量中显示字符．sprintf, num2str和int2str可以将数字转换成串．

字符变量通过括号连成大串．例如：

s='hello';

s=[s,'world']

s =

hello world

eval是与字符变量—起工作的函数，执行简单字符宏调用．eval( t )执行包含在t内的字符．如果t是任何Matlab表达式或语句的源字符，则字符串被解释执行．例如：

t='eye(2)', eval(t)

结果为：

ans=

1     0

0     1

又例如，给矩阵元素赋值

t='1/(i+j-1)';

for i=1:n

for j=1:n

               a( i,j)=eval(t);

end

end

这儿有一个例子，介绍如何一起使用eval与load命令，装入十个具有顺序文件名的文件中的数据：

fname='mydata';

for i = 1:10

eval([ 'load ', fname, int2str( i )])

end

§5.4.3外部程序

Matlab与外部独立程序的通讯方式可以是多种多样的，下面介绍其中的一个办法：

(1) Matlab中将变量存入磁盘

(2) 运行外部程序（读数据文件，进行处理），将结果写到磁盘上

(3) 将处理后的文件装回到工作空间中

例如：用外部程序gareqn找garfield方程的结果：

function y=garfield(a,b,q,r)

save gardata a,b,q,r

! gareqn

load gardata

使用FORTRAN或其它语言写gareqn程序，使其可以读gardata.mat，进行处理，将结果存入文件中．

这个程序可将计算机的“连接码”提供给Matlab，在许多系统中它将新的目标码连接到程序中比物理联接要方便得多．

§5.4.4输入输出数据

可使用各种方法将其它程序和外部世界的数据送入Matlab，同样可把Matlab数据输送到外部世界，使你的程序以Matlab使用的文件形式直接计算数据．

最好的方法取决于多少数据，数据是否可读，什么形式等：

(1) 清晰的元素表输入：

如果你有少量数据，比如说小于10~15个元素，使用方括号[]输入．

(2) 使用文本编辑建立命令文件，将数据列为清晰的元素表输入．如果数据不是可读形式，又不得不以一种方法键入，可以重复运行M文件，重复修改数据．

(3) 如果数据以ASCII形式存贮，并有固定长度，行尾有回车符，各数间有空格的文件称为flatfile(ASCII的flat file可由普通文本编辑来编辑)，flat file通过load命令直接读进Matlab，结果存入名为文件名的变量中去．

(4) 将数据文件译成Matlab文件形式，使用load命令，translate程序由Matlab中的应用程序库支持，translate程序将ASCII文件、二进制文件、FORTRAN非格式文件和DIF文件转换为Matlab使用的特定的MAT文件，当磁盘文件中存有大量数据时，这个方法输入最好．

Matlab数据输出到外部世界的方法：

(1) 小矩阵时：使用diary命令建立日志文件，在文件中列出变量，用文本编辑处理日志文件，日志的输出包括运行中的Matlab命令．

(2) 使用save命令存入变量，退出Matlab，用translate程序将MAT文件转换成任一种其它文件形式．

第六节　Matlab符号运算

 

Matlab本身并没有符号计算功能，1993年通过购买Maple的使用权后，开始具备符号运算的功能．符号运算的类型很多，几乎涉及数学的所有分支．

§6.1 Matlab符号运算的工作流程

§6.1.1工作过程

 

§6.1.2 核心工具

sym函数  VS    syms 语句

sym函数：构造符号变量和表达式:  a=sym('a')

（Construct symbolic numbers,variables and objects）

syms语句：构造符号对象的简捷方式（Short-cut forconstructing symbolic objects）

§6.1.3 符号变量确定原则

（1）除了i 和j之外，字母位置最接近x的字母；若距离相等，则取ASCII码大的；

（2）若没有除了i 与j以外的字母，则视x为默认的符号变量；

（3）可利用函数findsym(string,N)来询问在众多符号中，哪N个为符号变量．例如：键入findsym(3*a*b+y^2,1)，即可得到答案y．更多的例子见下表：

 

符号表达式	默认符号变量
a*x^2+b*x+c	x
1/(4+cos(t))	t
4*x/y	x
2*a+b	b
2*i	x

§6.2  Matlab的六大常见符号运算

§6.2.1 因式分解

symsx

f=x^6+1;

s=factor(f)

结果为：

s=(x^2+1)*(x^4-x^2+1)

§6.2.2 计算极限

求极限：

 

  syms h n x

  L=limit('(log(x+h)-log(x))/h',h,0)  %%单引号可省略掉

 M=limit('(1-x/n)^n',n,inf) 

结果为：

 L =1/x

 M =exp(-x) 

§6.2.3 计算导数

．

  syms a x;   y=sin(a*x);

 A=diff(y,x)

 B=diff(y,a) 

 C=diff(y,x,2)

结果为：

  A = cos(a*x)*a

 B = cos(a*x)*x

 C = -sin(a*x)*a^2  

§6.2.4 计算不定积分、定积分、反常积分

，

， 

．

    syms x

    f=(x^2+1)/(x^2-2*x+2)^2;

    g=cos(x)/(sin(x)+cos(x));

    h=exp(-x^2);

    I=int(f)

    J=int(g,0,pi/2)

K=int(h,0,inf) 

结果为：

  I =1/4*(2*x-6)/(x^2-2*x+2)+3/2*atan(x-1)

  J =1/4*pi

  K =1/2*pi^(1/2) 

 

§6.2.5 符号求和

求级数   的和S, 以及前十项的部分和S1．

    syms n

    S=symsum(1/n^2, 1, inf)

S1=symsum(1/n^2,1,10) 

结果为：

  S =1/6*pi^2

  S1 =1968329/1270080

重要说明：当求函数项级数   的和S2时，可用命令：

syms n x

S2=symsum(x/n^2, n, 1, inf)   

S2 =1/6*x*pi^2

 

两点说明：

（1）注意观察S2与S1的细微区别！

（2）当通项公式的Matlab表达式较长时，表达式要加上单引号．后面的练习中会遇到此问题．

§6.2.6 解代数方程和常微分方程

利用符号表达式解代数方程所需要的函数为solve(f)，即解符号方程式f．

例如：求一元二次方程a*x^2+b*x+c=0的根．

f=sym('a*x^2+b*x+c')  或 f='a*x^2+b*x+c'

solve(f)

ans=

[1/2/a*(-b+(b^2-4*c*a)^(1/2))]

[1/2/a*(-b-(b^2-4*c*a)^(1/2))]

solve(f,a)

ans=

-(b*x+c)/x^2

 

利用符号表达式可求解微分方程的解析解，所需要的函数为dsolve(f)，使用格式：

dsolve('equation1', ' equation2', …)

其中：equation为方程或条件．写方程或条件时，用 Dy 表示y 关于自变量的一阶导数，用D2y 表示 y 关于自变量的二阶导数，依此类推．

1． 求微分方程 的通解．

syms x y   %定义x,y为符号

dsolve('Dy=x',  'x')

ans =

1/2*x^2+C1

试比较：

若写成：

syms x y   %定义x,y为符号

dsolve('Dy=x')

结果将是什么？是否正确？为什么？

 

2． 求微分方程 的特解．

syms x y

dsolve('D2y=x+Dy','y(0)=1', 'Dy(0)=0', 'x')

ans =

-1/2*x^2+exp(x)-x

试比较：

若写成：

syms x y

dsolve('D2y=x+Dy','y(0)=1', 'Dy(0)=0')

 

结果将是什么？是否正确？为什么？

 

3． 求微分方程组 的通解．

syms x y

[x,y]=dsolve('Dx=y+x,Dy=2*x') 

x =

1/3*C1*exp(-t)+2/3*C1*exp(2*t)+1/3*C2*exp(2*t)-1/3*C2*exp(-t)

y =

2/3*C1*exp(2*t)-2/3*C1*exp(-t)+2/3*C2*exp(-t)+1/3*C2*exp(2*t)

 

试比较：

若写成：

（1）dsolve('Dx=y+x, Dy=2*x')

结果将是：

ans =

    x: [1x1 sym]

y: [1x1 sym]

试解释此结果的含义．

 

若写成：

（2）[x,y]=dsolve('Dx=y+x, Dy=2x')

结果将是：

x=

exp(t)*C1+C2*exp(t)-C2-2-2*t

y =

             C2+2*t

是否正确？为什么？
综合练习

 

1．按顺序进行如下的操作：

（1）产生一个5阶魔术方阵A；并计算A'与A^-1(即inv(A))；

（2）求A的特征值；

（3）计算A的各列的总和与平均值；

（4）计算A的各行的总和与平均值；

（5）若b=[1 2 3 45] '，求方程组 Ax=b的解；

（6）验证你的结论的正确性．

2．产生行向量S =[1.0, 1.2, 1.4, …,  20]，并计算S * S' 与  S' * S，你有何“发现”？

3．设A= ；B= ；求C=A * B – B * A，你有何“发现”？

4．若设矩阵A= ；B= ；求C=A * B – B * A，你又有何“发现”？

5．如何建立如下的矩阵（命令方式和程序方式）？

 

（１） ；　　（２） ；

 

 

（３） ；　　（４）

（５）                （6）  

（７）

6．绘制下列曲线的图形（散点图与折线图）：

7．绘制下列曲面的图形： （提示：曲面由两部分构成）

8．在同一个图形上作下列两个函数的图象：

（1） ；    （2）

9．假如你有一组实测数据，例如：

x=[53  56     60    67.5    75    90  110];

y=[109 120.5  130   141.1  157.5  180  185];

求其回归直线，画回归直线图形并计算最小误差平方和．

10．假如你有一组实测数据，例如：

x=[75 86  95  108  112  116 135  151  155 160  163  167 171  178 185];

y=[10 12  15  17   20   22   35   41   48  50   51   54  59   66   75];

求其回归直线，画回归直线图形并计算最小误差平方和．

11．随机产生500个0到100的整数FS作为学生的考试分数．

（1）       画出FS的简单直方图；

（2）       画出每个分数段（0～10、10～20、…，90～100）的统计频数直方图；

12．求下列各结果：

（1）用Matlab因式分解： ．

（2）用Matlab求极限： ．

（3）用Matlab求积分： ．

（4）用Matlab求幂级数： 的和函数（化简结果）．

13．非线性回归尝试

下表是到1994年的游泳世界纪录，试估计时间y与距离x的关系．

距离x（米）	50	100	200	400	800	1500
时间y（秒）	21.81	48.42	106.69	225	466.60	863.48

 

说明：用线性回归方法将得到： ，但当 时， ，这是非常荒唐的结果！显然，一个基本要求是当 时 ．试尝试使用非线性回归模型: .

 

14. (三维)符号作图尝试

命令	作用	解释
ezplot3	3-D parametric curve plotter	3D参数曲线图形
ezcontour	use contour plotter	等高线图
ezcontourf	filled contour plotter	填充等高线图
ezmesh	3-D mesh plotter	3D mesh曲面图形
ezmeshc	combination mesh/contour plotter	mesh曲面/等高线图
ezsurf	3-D colored surface plotter	3D surf曲面图形
ezsurfc	combination surf/contour plotter	surf曲面/等高线图

 

请尝试以下的命令：

ezplot3('sin(t)', ' cos(t)','t', [0,6*pi])

ezcontour('x*exp(-x^2 - y^2)')

ezcontourf('x*exp(-x^2 - y^2)')

ezmesh('(s-sin(s))*cos(t)','(1-cos(s))*sin(t)','s',[-2*pi,2*pi])

ezmeshc('(s-sin(s))*cos(t)','(1-cos(s))*sin(t)','s',[-2*pi,2*pi])

ezsurf('x*exp(-x^2 - y^2)')

ezsurfc('x*exp(-x^2 - y^2)')

 

 

绪  论

Matlab是“Matrix Laboratory”的缩写，意为“矩阵实验室”，是当今美国很流行的科学计算软件．信息技术、计算机技术发展到今天，科学计算在各个领域得到了广泛的应用．在许多诸如控制论、时间序列分析、系统仿真、图像信号处理等方面产生了大量的矩阵及其相应的计算问题．自己去编写大量的繁复的计算程序，不仅会消耗大量的时间和精力，减缓工作进程，而且往往质量不高．美国Mathwork软件公司推出的Matlab软件就是为了给人们提供一个方便的数值计算平台而设计的．

Matlab是一个交互式的系统，它的基本运算单元是不需指定维数的矩阵，按照IEEE的数值计算标准（能正确处理无穷数Inf(Infinity)、无定义数NaN(not-a-number)及其运算）进行计算．系统提供了大量的矩阵及其它运算函数，可以方便地进行一些很复杂的计算，而且运算效率极高．Matlab命令和数学中的符号、公式非常接近，可读性强，容易掌握，还可利用它所提供的编程语言进行编程完成特定的工作．除基本部分外，Matlab还根据各专门领域中的特殊需要提供了许多可选的工具箱，如应用于自动控制领域的Control System工具箱和神经网络中Neural Network工具箱等．

第一节　Matlab的安装及使用

§1.1 Matlab的安装

版本Matlab 6.x for Windows．因为它使用方便，界面美观，我们选择它作为主要讲解版本．Matlab还有许多附加的部分，最常见的部分称为Simulink，是一个用作系统仿真的软件包，它可以让您定义各种部件，定义各自对某种信号的反应方式及与其它部件的连接方式．最后选择输入信号，系统会仿真运行整个模拟系统，并给出统计数据．Simulink有时是作为Matlab的一部分提供的，称为Matlab with Simulink版本．Matlab还有许多工具箱，它们是根据各个特殊领域的需要，用Matlab自身的语言编写的程序集，使用起来非常方便．您可以视工作性质和需要购买相应的工具箱．常见的工具箱有：

Signal Process	信号处理	System Identification	系统辨识
Optimization	优化	Neural Network	神经网络
Control System	自动控制	Spline	样条
Symbolic Math	符号代数	Image Process	图像处理
Nonlinear Control	非线性控制	Statistics	统计

§1.2 Matlab基本用法

从Windows中双击Matlab图标，会出现Matlab命令窗口（Command Window），在一段提示信息后，出现系统提示符“>>”．您可以在提示符后键入各种命令，通过上下箭头可以调出以前打入的命令，用滚动条可以查看以前的命令及其输出信息．

如果对一条命令的用法有疑问的话，可以用Help菜单中的相应选项查询有关信息，也可以用help命令在命令行上查询，您可以试一下help、help help和help eig（求特征值的函数）命令．

下面我们先从输入简单的矩阵开始掌握Matlab的功能．

§1.2.1输入简单的矩阵

输入一个小矩阵的最简单方法是用直接排列的形式．矩阵用方括号括起，元素之间用空格或逗号分隔，矩阵行与行之间用分号分开．例如输入：

A=[1  2  3 ;  4  5  6 ;  7  8  0]

系统会回答

A =

1      2     3

4      5     6

7      8     0

表示系统已经接收并处理了命令，在当前工作区内建立了矩阵A．

大的矩阵可以分行输入，用回车键代替分号，如：

A=[ 1     2     3

4    5     6

7    8     0]

结果和上式一样，也是

A =

1      2     3

4      5     6

7      8     0

§1.2.2矩阵元素

Matlab的矩阵元素可以是任何数值表达式．如：

x=[ -1.3   sqrt(3)   (1+2+3)*4/5]

结果：

x =

 -1.3000   1.7321    4.8000

在括号中加注下标，可取出单独的矩阵元素．如：

x(5)=abs(x(1))

结果

x =

 -1.3000   1.7321    4.8000    0     1.3000

注：结果中自动产生了向量的第5个元素，中间未定义的元素自动初始为零．

大的矩阵可把小的矩阵作为其元素来完成，如：

A=[A; [10  11  12]]

结果

A =

1      2     3

4      5     6

7      8     0

10    11   12

小矩阵可用“:”从大矩阵中抽取出来，如：

A=A(1:3,:);

即从A中取前三行和所有的列，重新组成原来的A. (详细介绍参见第二节的相关内容）

§1.2.3语句和变量

Matlab的表述语句、变量的类型说明由Matlab系统解释和判断．Matlab语句通常形式为：

变量=表达式

或者使用其简单形式为：

表达式

表达式由操作符或其它特殊字符、函数和变量名组成．表达式的结果为一个矩阵，显示在屏幕上，同时保存在变量中以留用．如果变量名和“=”省略，则具有ans名(意思指回答)的变量将自动建立．例如：

键入1900/81

结果为：

ans =

23.4568

需注意的问题有以下几点：

l 语句结束键入回车键，若语句的最后一个字符是分号，即“;”，则表明不输出当前命令的结果．

l 如果表达式很长，一行放不下，可以键入“ …”（三个点，但前面必须有个空格，目的是避免将形如“数2 …”理解为“数2.”与“..”的连接，从而导致错误），然后回车．

l 变量和函数名由字母加数字组成，但最多不能超过63个字符，否则系统只承认前63个字符．

l Matlab变量字母区分大小写，如A和a不是同一个变量，函数名一般使用小写字母，如inv(A)不能写成INV(A)，否则系统认为未定义函数．

§1.2.4 who和系统预定义变量

输入who命令可检查工作空间中建立的变量，键入：

who

系统输出为：

Your variables are:

A  ans  x

这里表明三个变量已由前面的例子产生了．

但列表中列出的并不是系统全部的变量，系统还有以下内部变量：

eps、pi、Inf、NaN

变量eps在决定诸如矩阵的奇异性时，可作为一个容许差，容许差的初值为1.0到1.0以后计算机所能表示的下一个最大浮点数，IEEE在各种计算机、工作站和个人计算机上使用这个算法．用户可将此值置为任何其它值(包括0值)．

变量pi是p.

Inf表示无穷大．如果您想计算1/0

S=1/0

结果会是

Warning：Divide by zero

S=Inf

具有IEEE规则的机器，被零除后，并不引出出错条件或终止程序的运行，而产生一个警告信息和一个特殊值在计算方程中列出来．

变量NaN表示它是个不定值．由Inf/Inf或0/0运算产生．

要了解当前变量的信息请键入whos，屏幕将显示：

  Name           Size               Bytes             Class

  A                4x3                96                  double array

  S                 1x1                8                  double array

  ans                     1x1                8                  double array

  x                 1x5                40                  double array

Grand total is 19 elementsusing 152 bytes

从size及bytes项目可以看出，每一个矩阵实元素需8个字节的内存．4×3的矩阵使用96个字节，全部变量的使用内存总数为152个字节．自由空间的大小决定了系统变量的多少，如计算机上有虚拟内存的话，其可定义的变量个数会大大增加．

§1.2.5数和算术表达式

Matlab中数的表示方法和一般的编程语言没有区别．如：

3                    -99                       0.0001

9.63972           1.6021E-20           6.02252e23

在计算中使用IEEE浮点算法其舍入误差是eps．浮点数表示范围是10^-308~10³⁰⁸．

数学运算符有：

+            加

-             减

*            乘

/             右除

\             左除

^            幂

这里1/4和4\1有相同的值都等于0.25（注意比较：1\4=4）．只有在矩阵的除法时左除和右除才有区别．

§1.2.6复数与矩阵

在Matlab中输入复数首先应该建立复数单位：

i=sqrt(-1)

及j=sqrt(-1)

之后复数可由下面语句给出：

Z=3+4i    (注意： 在4与i之间不要留有任何空间!)

输入复数矩阵有两个方便的方法，如：

A=[1  2; 3  4] + i*[5 6; 7  8]

和 A=[1+5i  2+6i;  3+7i  4+8i]

两式具有相等的结果．但当复数作为矩阵的元素输入时，不要留有任何空间，如1+5i，如在“+”号左右留有空格，就会被认为是两个分开的数．

不过实际使用复数时并没有这么麻烦，系统有一个名为startup.m的Matlab命令文件，建立复数单位的语句也放在其中．当Matlab启动时，此文件自动执行，i和j将自动建立．

§1.2.7输出格式

任何Matlab语句执行结果都可在屏幕上显示，同时赋给指定的变量，没有指定变量时赋给ans．数字显示格式可由format命令来控制（Windows系统下的Matlab系统的数字显示格式可以由Option菜单中的NumericalFormat菜单改变）．format仅影响矩阵的显示，不影响矩阵的计算与存贮．(Matlab以双精度执行所有的运算)

首先，如果矩阵元素是整数则矩阵显示就没有小数，如x=[-1 0 1]，结果为：

        x=

               -1    0     1

如果矩阵元素不是整数则输出形式有：（用命令：format 格式进行切换）

格式	中文解释	说明
format	短格式 （缺省格式）	Default. Same as SHORT
format short	短格式 （缺省格式）	Scaled fixed point format with 5 digits （只显示五位十进制数）
format long	长格式	Scaled fixed point format with 15 digits
format short e	短格式e方式	Floating point format with 5 digits
format long e	长格式e方式	Floating point format with 15 digits
format short g	短格式g方式	Best of fixed or floating point format with 5 digits
format long g	长格式g方式	Best of fixed or floating point format with 15 digits
format hex	16进制格式	Hexadecimal format
format +	+格式	The symbols +, - and blank are printed                         for positive, negative and zero elements.                        Imaginary parts are ignored
format bank	银行格式	Fixed format for dollars and cents
format rat	有理数格式	Approximation by ratio of small integers
format compact	压缩格式	Suppress extra line-feeds
format loose	自由格式	Puts the extra line-feeds back in

例如：

x=[4/3  1.2345e-6]

在不同的输出格式下的结果为：

短格式                   1.3333                                       0.0000

短格式e方式         1.3333e+000                               1.234e-006

长格式                   1.333333333333333                    0.000001234500000

长格式e方式         1.333333333333333e-000            1.23450000000000e-006

有理数格式            4/3                                             1/810045

16进制格式           3ff5555555555555                      3eb4b6231abfd271

+格式                    +                                               +

对于短格式，如果矩阵的最大元素比数999999999大，或者比数0.0001小，则在打印时，将加入一个普通的长度因数．如y=1.e20*x，意为x被10²⁰乘，结果为：

y=

         1.0e+020*

1.3333             0.0000

“+”格式是显示大矩阵的一种紧凑方法，“+”，“-”和空格显示正数、负数和零元素．

最后format compact命令压缩显示的矩阵，以允许更多的信息显示在屏幕上．

§1.2.8 Help求助命令和联机帮助

Help求助命令很有用，它对Matlab大部分命令提供了联机求助信息．您可以从Help菜单中选择相应的菜单，打开求助信息窗口查询某条命令，也可以直接用help命令．

键入help

得到help列表文件，键入“help 指定项目”，如：

键入help eig

则提供特征值函数的使用信息．

键入help [

显示如何使用方括号等．

键入help help

显示如何利用help本身的功能．

还有，键入lookfor  <关键字>：可以从m文件的help中查找有关的关键字．

§1.2.9 退出和存入工作空间

退出Matlab可键入quit或exit或选择相应的菜单．中止Matlab运行会引起工作空间中变量的丢失，因此在退出前，应键入save命令，保存工作空间中的变量以便以后使用．

键入 save

则将所有变量作为文件存入磁盘Matlab.mat中，下次Matlab启动时，

键入load

将变量从Matlab.mat中重新调出．

save和load后边可以跟文件名或指定的变量名，如仅有save时，则只能存入Matlab.mat中．如save temp命令，则将当前系统中的变量存入temp.mat中去，命令格式为：

savetemp x 仅仅存入x变量．

savetemp X Y Z 则存入X、Y、Z变量．

load  temp可重新从temp.mat文件中提出变量，load也可读ASCII数据文件．详细语法见联机帮助．

第二节　向量与矩阵运算

Matlab能处理数、向量和矩阵．但一个数事实上是一个1×1的矩阵，1个n维向量也不过是一个1×n或n×1的矩阵．从这个角度上来讲，Matlab处理的所有的数据都是矩阵．Matlab的矩阵处理能力是非常灵活、强大的．以下我们将从矩阵的产生、基本运算、矩阵函数等几个方面来说明．

§2.1向量及矩阵的生成

除了我们在上节介绍的直接列出矩阵元素的输入方法，矩阵还可以通过几种不同的方式输入到Matlab中．

§2.1.1 通过语句和函数产生

1. 向量的产生

除了直接列出向量元素（即所谓的“穷举法”）外，最常用的用来产生相同增量的向量的方法是利用“:”算符（即所谓的“描述法”）．在Matlab中，它是一个很重要的字符．如：

z=1:5

z =

     1    2     3    4     5

即产生一个1~5的单位增量是1的行向量，此为默认情况.

用“:”号也可以产生单位增量不等于1的行向量，语法是把增量放在起始量和结尾量的中间．如：

x=0:pi/4:pi

即产生一个由0~pi的行向量，单位增量是pi/4=3.1416/4=0.7854．

x =

 0     0.7854           1.5708           2.3562           3.1416

也可以产生单位增量为负数的行向量．如：

y=6:-1:1

y =

     6    5     4     3    2     1

2. 矩阵的产生

Matlab提供了一批产生矩阵的函数：

zeros	产生一个零矩阵	diag	产生一个对角矩阵
ones	生成全1矩阵	tril	取一个矩阵的下三角
eye	生成单位矩阵	triu	取一个矩阵的上三角
magic	生成魔术方阵	pascal	生成PASCAL矩阵

例如：

ones(3)

ans =

 1    1     1

 1    1     1

 1    1     1

eye(3)

ans =

 1    0     0

 0    1     0

 0    0     1

除了以上产生标准矩阵的函数外，Matlab还提供了产生随机（向量）矩阵的函数rand和randn，及产生均匀级数的函数linspace、产生对数级数的函数logspace和产生网格的函数meshgrid等等．详细使用请查阅随机文档．

“ : ”冒号可以用来产生简易的表格，为了产生纵向表格形式，首先用冒号“ : ”产生行向量，再进行转置，计算函数值的列，然后形成有二列的矩阵．例如命令：

x=(0.0:0.2:3.0)';

y=exp(-x).*sin(x);

[x y]

产生结果为：

ans =

0             0

0.2000      0.1627

0.4000      0.2610

0.6000      0.3099

0.8000      0.3223

1.0000      0.3096

1.2000      0.2807

1.4000      0.2430

1.6000      0.2018

1.8000      0.1610

2.0000      0.1231

2.2000      0.0896

2.4000      0.0613

2.6000      0.0383

2.8000      0.0204

3.0000      0.0070

§2.1.2 通过后缀为.m的命令文件产生

如有文件data.m，其中包括正文：

A=[  1 2 3

4 5 6

7 8 0]

则用data命令执行data.m，可以产生名为A的矩阵．

§2.2 矩阵操作

在Matlab中可以对矩阵进行任意操作，包括改变它的形式，取出子矩阵，扩充矩阵，旋转矩阵等．其中最重要的操作符为“:”， 它的作用是取出选定的行与列．

例如：

A(:,:) 代表A的所有元素；试比较A(:),将A按列的方向拉成长长的1列（向量）；

A(:,J) 代表A的第J列；

A(J:K) 代表 A(J), A(J+1), …, A(K)，如同A(:)的第J到第K个元素；

A(:,J:K) 代表A(:,J), A(:,J+1), …, A(:,K)，如此类推.

对矩阵可以进行各种各样的旋转、变形、扩充：

 

矩阵的转置用符号“ ' ”表示：

如A=[1  2  3;  4  5  6 ;  7  8  0]

那么：计算B=A'

B =

 1    4     7

 2    5     8

 3    6     0

符号“ ' ”为矩阵的转置，如果Z为复矩阵，则Z'为它的复数共轭转置，非共轭转置使用Z.' 或conj(Z')求得．

reshape改变矩阵的形状，这是什么意思呢？可举一个例子来说明．

A=[A;[10 11 12]]

A =

 1    2     3

 4    5     6

 7    8     0

 10 11   12

则 reshape(A,2,6)

ans =

 1    7     2     8    3     0

 4    10   5     11   6     12

可见，reshape 是将矩阵元素以列为单位进行重组，原来4×3的矩阵变为了2×6的矩阵．

第三节　矩阵的基本运算

§3.1 加和减

如矩阵A和B的维数相同，则A+B与A-B表示矩阵A与B的和与差．如果矩阵A和B的维数不匹配，Matlab会给出相应的错误提示信息．如：

A=                                B=

1     2     3                   1     4     7

4     5     6                   2     5     8

7     8     0                   3     6     0

C =A+B返回：

C =

    2     6    10

    6    10    14

   10    14     0

如果运算对象是个标量（即1×1矩阵），可和其它矩阵进行加减运算．例如：

x=    -1                  y=x-1=    -2

0                                 -1

2                                 1

§3.2矩阵乘法

Matlab中的矩阵乘法有通常意义上的矩阵乘法，也有Kronecker乘法，以下分别介绍．

§3.2.1 矩阵的普通乘法

矩阵乘法用“ * ”符号表示，当A矩阵列数与B矩阵的行数相等时，二者可以进行乘法运算，否则是错误的．计算方法和线性代数中所介绍的完全相同．

如：A=[1  2 ; 3  4]; B=[5  6 ; 7  8];  C=A*B，

结果为

C= × = =

即Matlab返回：

C =

    19   22

    43   50

如果A或B是标量，则A*B返回标量A（或B）乘上矩阵B（或A）的每一个元素所得的矩阵．

§3.3 矩阵除法

在Matlab中有两种矩阵除法符号：“＼”即左除和“／”即右除．如果A矩阵是非奇异方阵，则A\B是A的逆矩阵乘B，即inv(A)*B；而B/A是B乘A的逆矩阵，即B*inv(A)．具体计算时可不用逆矩阵而直接计算．

通常：

x=A\B就是A*x=B的解；

x=B/A就是x*A=B的解.

当B与A矩阵行数相等可进行左除．如果A是方阵，用高斯消元法分解因数．解方程：A*x(:, j)=B(:, j)，式中的(:, j)表示B矩阵的第j列，返回的结果x具有与B矩阵相同的阶数，如果A是奇异矩阵将给出警告信息．

如果A矩阵不是方阵，可由以列为基准的Householder正交分解法分解，这种分解法可以解决在最小二乘法中的欠定方程或超定方程，结果是m×n的x矩阵．m是A矩阵的列数，n是B矩阵的列数．每个矩阵的列向量最多有k个非零元素，k 是A的有效秩．

右除B/A可由B/A=(A'\B')'左除来实现．

§3.4矩阵乘方

A^P意思是A的P次方．如果A是一个方阵，P是一个大于1的整数，则A^P表示A的P次幂，即A自乘P次．如果P不是整数，计算涉及到特征值和特征向量的问题，如已经求得：[V,D]=eig(A)，则：

A^P=V*D.^P/V（注：这里的.^表示数组乘方，或点乘方，参见后面的有关介绍）

如果B是方阵， a是标量，a^B就是一个按特征值与特征向量的升幂排列的B次方程阵． 如果a和B都是矩阵，则a^B是错误的．

§3.5 矩阵的超越函数

在Matlab中解释exp(A)和sqrt(A)时曾涉及到级数运算，此运算定义在A的单个元素上． Matlab可以计算矩阵的超越函数，如矩阵指数、矩阵对数等．

一个超越函数可以作为矩阵函数来解释，例如将“m”加在函数名的后边而成expm(A)和sqrtm(A)，当Matlab运行时，有下列三种函数定义：

expm           矩阵指数

logm        矩阵对数

sqrtm           矩阵开方

所列各项可以加在多种m文件中或使用funm．请见应用库中sqrtm.m，1ogm.m，funm.m文件和命令手册．

§3.6数组运算

数组运算由线性代数的矩阵运算符“*”、“/”、“\”、“^”前加一点来表示，即为“.*”、“./”、“.\”、“.^”．注意没有“.+”、“.-”运算．

§3.6.1数组的加和减

对于数组的加和减运算与矩阵运算相同，所以“+”、“-”既可被矩阵接受又可被数组接受．

§3.6.2数组的乘和除

数组的乘用符号.*表示，如果A与B矩阵具有相同阶数，则A.*B表示A和B单个元素之间的对应相乘．例如x=[1  2  3]; y=[ 4  5  6];

计算z=x.*y

结果z=4  10  18

数组的左除（.\）与数组的右除（./），由读者自行举例加以体会．

§3.6.3 数组乘方

数组乘方用符号.^表示．

例如：键入：

x=[ 1  2  3]

y=[ 4  5  6]

则z=x.^y=[1^4  2^5  3^6]=[1  32  729]

(1) 如指数是个标量，例如x.^2，x同上，则：

z=x.^2=[1^2  2^2  3^2]=[ 1  4  9]

(2) 如底是标量，例如2 .^[x y] ，x、y同上，则：

z=2 .^[x y]=[2^1  2^2  2^3  2^4 2^5  2^6]=[2  4  8  16  32  64]

从此例可以看出Matlab算法的微妙特性，虽然看上去与其它乘方没什么不同，但在2和“．”之间的空格很重要，如果不这样做，解释程序会把“．”看成是2的小数点． Matlab看到符号“^”时，就会当做矩阵的幂来运算，这种情况就会出错，因为指数矩阵不是方阵．

§3.7 矩阵函数

Matlab的数学能力大部分是从它的矩阵函数派生出来的，其中一部分装入Matlab本身处理中，它从外部的Matlab建立的M文件库中得到，还有一些由个别的用户为其自己的特殊的用途加进去的．其它功能函数在求助程序或命令手册中都可找到．手册中备有为Matlab提供数学基础的LINPACK和EISPACK软件包，提供了下面四种情况的分解函数或变换函数：

(1)三角分解；(2)正交变换；(3) 特征值变换；(4)奇异值分解．

§3.7.1三角分解

最基本的分解为“LU”分解，矩阵分解为两个基本三角矩阵形成的方阵，三角矩阵有上三角矩阵和下三角矩阵．计算算法用高斯变量消去法．

从lu函数中可以得到分解出的上三角与下三角矩阵，函数inv得到矩阵的逆矩阵，det得到矩阵的行列式．解线性方程组的结果由方阵的“\”和“/”矩阵除法来得到．

例如：

A=[ 1      2     3

4      5     6

7      8     0]

LU分解，用Matlab的多重赋值语句

[L,U]=lu(A)

得出

L =

0.1429	1.0000	0
0.5714	0.5000	1.0000
1.0000	0	0

U =

7.0000	8.0000	0
0	0.8571	3.0000
0	0	4.5000

注：L是下三角矩阵的置换，U是上三角矩阵的正交变换，分解作如下运算，检测计算结果只需计算L*U即可．

求逆由下式给出： x=inv(A)

x =

-1.7778	0.8889	-0.1111
1.5556	-0.7778	0.2222
-0.1111	0.2222	-0.1111

从LU分解得到的行列式的值是精确的，d=det(U)*det(L)的值可由下式给出：

d=det(A)

d =

27

直接由三角分解计算行列式：d=det(L)*det(U)

d =

27.0000

为什么两种d的显示格式不一样呢? 当Matlab做det(A)运算时，所有A的元素都是整数，所以结果为整数．但是用LU分解计算d时，L、U的元素是实数，所以Matlab产生的d也是实数．

例如：线性联立方程取 b=[ 1

                                           3

                                           5]

解Ax=b方程，用Matlab矩阵除得到

x=A\b

结果x=

0.3333

0.3333

0.0000

由于A=L*U，所以x也可以有以下两个式子计算：y=L\b，x=U\y．得到相同的x值，中间值y为：

y =

5.0000

0.2857

0.0000

Matlab中与此相关的函数还有rcond、chol和rref．其基本算法与LU分解密切相关．chol函数对正定矩阵进行Cholesky分解，产生一个上三角矩阵，以使R'*R=X．rref用具有部分主元的高斯－约当消去法产生矩阵A的化简梯形形式．虽然计算量很少，但它是很有趣的理论线性代数．为了教学的要求，也包括在Matlab中．

§3.7.2正交变换

“QR”分解用于矩阵的正交－三角分解．它将矩阵分解为实正交矩阵或复酉矩阵与上三角矩阵的积，对方阵和长方阵都很有用．

例如A=[  1     2     3

4     5     6

7     8     9

10    11    12]

是一个降秩矩阵，中间列是其它二列的平均，我们对它进行QR分解：

[Q,R]=qr(A)

Q =

-0.0776	-0.8331	0.5444	0.0605
-0.3105	-0.4512	-0.7709	0.3251
-0.5433	-0.0694	-0.0913	-0.8317
-0.7762	0.3124	0.3178	0.4461

R =

-12.8841	-14.5916	-16.2992
0	-1.0413	-2.0826
0	0	0.0000
0	0	0

可以验证Q*R就是原来的A矩阵．由R的下三角都给出0，并且R(3,3)=0.0000，说明矩阵R与原来矩阵A都不是满秩的．

下面尝试利用QR分解来求超定和降秩的线性方程组的解．

例如：

b=[ 1

3

5

7]

讨论线性方程组Ax=b，我们可以知道方程组是超定的，采用最小二乘法的最好结果是计算x=A\b．

结果为：

Warning: Rank deficient, rank = 2 tol =   1.4594e-014

x =

    0.5000

         0

    0.1667

我们得到了缺秩的警告．用QR分解法计算此方程组分二个步骤：

y=Q'*b

x=R\y

求出的y值为

y =

-9.1586

-0.3471

0.0000

0.0000

x的结果为

Warning: Rank deficient, rank = 2 tol =   1.4594e-014

x =

    0.5000

         0

    0.1667

用A*x来验证计算结果，我们会发现在允许的误差范围内结果等于b．这告诉我们虽然联立方程Ax=b是超定和降秩的，但两种求解方法的结果是一致的．显然x向量的解有无穷多个，而“QR”分解仅仅找出了其中之一．

§3.7.3奇异值分解

在Matlab中三重赋值语句

[U,S,V]=svd(A)

在奇异值分解中产生三个因数：

A=U*S*V '

U矩阵和V矩阵是正交矩阵，S矩阵是对角矩阵，svd(A)函数恰好返回S的对角元素，而且就是A的奇异值（其定义为：矩阵A'*A的特征值的算术平方根）．注意到A矩阵可以不是方的矩阵．

奇异值分解可被其它几种函数使用，包括广义逆矩阵pinv(A)、秩rank(A)、欧几里德矩阵范数norm(A,2)和条件数cond(A)．

§3.7.4 特征值分解

如果A是n×n矩阵，若l满足Ax=lx，则称l为A的特征值，x为相应的特征向量．

函数eig(A)返回特征值列向量，如果A是实对称的，特征值为实数．特征值也可能为复数，例如：

A=[  0     1

-1     0]

eig(A)

产生结果

ans =

0 + 1.0000i

0 - 1.0000i

如果还要求求出特征向量，则可以用eig(A)函数的第二个返回值得到：

[x,D]=eig(A)

D的对角元素是特征值．x的列是相应的特征向量，以使A*x=x*D．

计算特征值的中间结果有两种形式：

Hessenberg形式为hess(A)，Schur形式为schur(A)．

schur形式用来计算矩阵的超越函数，诸如sqrtm(A)和logm(A)．

如果A和B是方阵，函数eig(A,B)返回一个包含一般特征值的向量来解方程

Ax=lBx

双赋值获得特征向量

[X,D]=eig(A,B)

产生特征值为对角矩阵D．满秩矩阵X的列相应于特征向量，使A*X=B*X*D，中间结果由qz(A,B)提供．

§3.7.5秩

Matlab计算矩阵A的秩的函数为rank(A)，与秩的计算相关的函数还有：rref(A)、orth(A)、null(A)和广义逆矩阵pinv(A)等．

利用rref(A)，A的秩为非0行的个数．rref方法是几个定秩算法中最快的一个，但结果上并不可靠和完善．pinv(A)是基于奇异值的算法．该算法消耗时间多，但比较可靠．其它函数的详细用法可利用Help求助．

第四节　Matlab中的图形

§4.1　二维作图

绘图命令plot绘制x-y坐标图；polor命令绘制极坐标图．

§4.1.1 基本形式

如果y是一个向量，那么plot(y)绘制一个y中元素的线性图．假设我们希望画出

y=[0., 0.48,  0.84,  1., 0.91,  6.14 ]

则用命令：plot(y)

它相当于命令：plot(x, y)，其中x=[1,2,…,n]或x=[1;2;…;n]，即向量y的下标编号, n为向量y的长度

Matlab会产生一个图形窗口，显示如下图形，请注意：坐标x和y 是由计算机自动绘出的．

图4.1.1.1plot([0.,0.48,0.84,1.,0.91,6.14])

 

上面的图形没有加上x轴和y轴的标注，也没有标题．用xlabel，ylabel，title命令可以加上．

如果x，y是同样长度的向量，plot(x,y)命令可画出相应的x元素与y元素的x-y坐标图．例：

x=0:0.05:4*pi;  y=sin(x);  plot(x,y)

grid on, title(' y=sin( x ) 曲线图' )

xlabel(' x = 0 : 0.05 : 4Pi ')

结果见下图．

图4.1.1.2 y=sin(x)的图形

 

title	图形标题
xlabel	x坐标轴标注
ylabel	y坐标轴标注
text	标注数据点
grid	给图形加上网格
hold	保持图形窗口的图形

表4.1.1.1 Matlab图形命令

§4.1.2 多重线

在一个单线图上，绘制多重线有三种办法.

第一种方法是利用plot的多变量方式绘制：

plot(x1,y1,x2,y2,...,xn,yn)

x1,y1,x2,y2,...,xn,yn是成对的向量，每一对x, y在图上产生如上方式的单线．多变量方式绘图是允许不同长度的向量显示在同一图形上．

第二种方法也是利用plot绘制，但加上hold on/off命令的配合：

plot(x1,y1)

hold on

plot(x2,y2)

hold off

第三种方法还是利用plot绘制，但代入矩阵：

如果plot用于两个变量plot(x,y)，并且x，y是矩阵，则有以下情况：

（1）如果y是矩阵，x是向量，plot(x,y)用不同的画线形式绘出y的行或列及相应的x向量，y的行或列的方向与x向量元素的值选择是相同的．

（2）如果x是矩阵，y是向量，则除了x向量的线族及相应的y向量外，以上的规则也适用．

（3）如果x，y是同样大小的矩阵，plot(x,y)绘制x的列及y相应的列．

还有其它一些情况，请参见Matlab的帮助系统．

§4.1.3 线型和颜色的控制

如果不指定划线方式和颜色，Matlab会自动为您选择点的表示方式及颜色．您也可以用不同的符号指定不同的曲线绘制方式．例如：

plot(x,y,'*')                    用'*'作为点绘制的图形．

plot(x1,y1,':',x2,y2,'+')    用':'画第一条线，用'+'画第二条线．

线型、点标记和颜色的取值有以下几种：

线型		点标记		颜色
-	实线	.	点	y	黄
:	虚线	o	小圆圈	m	棕色
-.	点划线	x	叉子符	c	青色
--	间断线	+	加号	r	红色
		*	星号	g	绿色
		s	方格	b	蓝色
		d	菱形	w	白色
		^	朝上三角	k	黑色
		v	朝下三角
		>	朝右三角
		<	朝左三角
		p	五角星
		h	六角星

表4.1.3.1线型和颜色控制符

如果你的计算机系统不支持彩色显示，Matlab将把颜色符号解释为线型符号，用不同的线型表示不同的颜色．颜色与线型也可以一起给出，即同时指定曲线的颜色和线型．

例如：     t=-3.14:0.2:3.14;

x=sin(t);  y=cos(t);

plot(t,x, '+r',t,y, '-b')

图4.1.3.1不同线型、颜色的sin,cos图形

 

§4.1.4对数图、极坐标图及条形图

 

 

例6：极坐标绘图

t=0:0.01:2*pi;

polar(t,sin(6*t))   %% (图4.1.4.3)

    

             图4.1.4.3 极坐标绘图         图4.1.4.4正态分布的统计直方图与其正态分布拟合曲线

 

例7：正态分布图

我们可以用命令normrnd生成符合正态分布的随机数．

normrnd(u,v,m,n)

其中，u表示生成随机数的期望，v代表随机数的方差．

运行：

a=normrnd(10,2,10000,1);

histfit(a)   %% (图4.1.4.4)

我们可以得到正态分布的统计直方图与其正态分布拟合曲线．

 

例8：比较正态分布（图4.1.4.5（1））与平均分布（图4.1.4.5（2））的分布图：

yn=randn(30000,1);    %% 正态分布

x=min(yn): 0.2 : max(yn);

subplot(121)

hist(yn,x)

yu=rand(30000,1);     %% 平均分布

subplot(122)

hist(yu,25)

4.1.4.5(1)                          4.1.4.5(2)

 

图4.1.4.5 正态分布与平均分布的分布图

§4.1.5 子图

在绘图过程中，经常要把几个图形在同一个图形窗口中表现出来，而不是简单地叠加（例如上面的例8）．这就用到函数subplot．其调用格式如下：

subplot(m,n,p)

subplot函数把一个图形窗口分割成m×n个子区域，用户可以通过参数p调用个各子绘图区域进行操作．子绘图区域的编号为按行从左至右编号．

 

例9：绘制子图

x=0:0.1*pi:2*pi;

subplot(2,2,1)

plot(x,sin(x),'-*');

title('sin(x)');

subplot(2,2,2)

plot(x,cos(x),'--o');

title('cos(x)');

subplot(2,2,3)

plot(x,sin(2*x),'-.*');

title('sin(2x)');

subplot(2,2,4);

plot(x,cos(3*x),':d')

title('cos(3x)')

得到图形如下：

图4.1.5.1子图

§4.1.6 填充图

利用二维绘图函数patch，我们可绘制填充图．绘制填充图的另一个函数为fill．

下面的例子绘出了函数humps（一个Matlab演示函数）在指定区域内的函数图形．

例10：用函数patch绘制填充图

fplot('humps',[0,2],'b')

hold on

patch([0.5 0.5:0.02:1 1],[0humps(0.5:0.02:1) 0],'r');

hold off

title('A region under an interesting function.')

grid

图4.1.6.1填充图

 

我们还可以用函数fill来绘制类似的填充图．

例11：用函数fill绘制填充图

x=0:pi/60:2*pi;

y=sin(x);

x1=0:pi/60:1;

y1=sin(x1);

plot(x,y,'r');

hold on

fill([x1 1],[y10],'g')

图4.1.6.2填充图

§4.2 三维作图

§4.2.1 mesh(Z)语句

mesh(Z)语句可以给出矩阵Z元素的三维消隐图，网络表面由Z坐标点定义，与前面叙述的x-y平面的线格相同，图形由邻近的点连接而成．它可用来显示用其它方式难以输出的包含大量数据的大型矩阵，也可用来绘制Z变量函数．

显示两变量的函数Z=f(x,y)，第一步需产生特定的行和列的x-y矩阵．然后计算函数在各网格点上的值．最后用mesh函数输出．

下面我们绘制sin(r)/r函数的图形．建立图形用以下方法：

x=-8:.5:8;

y=x';

x=ones(size(y))*x;

y=y*ones(size(y))';

R=sqrt(x.^2+y.^2)+eps;

z=sin(R)./R;

mesh(z)       %%  试运行mesh(x,y,z)，看看与mesh(z)有什么不同之处？

各语句的意义是：首先建立行向量x，列向量y；然后按向量的长度建立1-矩阵；用向量乘以产生的1-矩阵，生成网格矩阵，它们的值对应于x-y坐标平面；接下来计算各网格点的半径；最后计算函数值矩阵Z．用mesh函数即可以得到图形．

图4.2.1三维消隐图

第一条语句x的赋值为定义域，在其上估计函数；第三条语句建立一个重复行的x矩阵，第四条语句产生y的响应，第五条语句产生矩阵R（其元素为各网格点到原点的距离）．用mesh方法结果如上．

另外，上述命令系列中的前4行可用以下一条命令替代：

[x, y]=meshgrid(-8:0.5:8)

§4.2.2与mesh相关的几个函数

(1) meshc与函数mesh的调用方式相同，只是该函数在mesh的基础上又增加了绘制相应等高线的功能．下面来看一个meshc的例子：

[x,y]=meshgrid([-4:.5:4]);

z=sqrt(x.^2+y.^2);

meshc(z)       %% 试运行 meshc(x,y,z)，看看与meshc(z)有什么不同之处？

我们可以得到图形：

图4.2.2.1  meshc 图

 

地面上的圆圈就是上面图形的等高线．

(2) 函数meshz与mesh的调用方式也相同，不同的是该函数在mesh函数的作用之上增加了屏蔽作用，即增加了边界面屏蔽．例如：

[x,y]=meshgrid([-4:.5:4]);

z=sqrt(x.^2+y.^2);

meshz(z)       %% 试运行 meshz(x,y,z)，看看与meshz(z)有什么不同之处？

我们得到图形：

图4.2.2.2  meshz 图

 

§4.2.3 其它的几个三维绘图函数

(1) 在Matlab中有一个专门绘制圆球体的函数sphere，其调用格式如下：

[x,y,z]=sphere(n) 

此函数生成三个(n+1)×(n+1)阶的矩阵，再利用函数surf(x,y,z)可生成单位球面．

[x,y,z]=sphere   此形式使用了默认值n=20

sphere(n)   只绘制球面图，不返回值．

运行下面程序：

sphere(30);

axis square;

我们得到球体图形：

图4.2.3.1  球面图

 

若只输入sphere画图，则是默认了n=20的情况．

 

(2) surf函数也是Matlab中常用的三维绘图函数．其调用格式如下：

surf(x,y,z,c)

输入参数的设置与mesh相同，不同的是mesh函数绘制的是一网格图，而surf绘制的是着色的三维表面．Matlab语言对表面进行着色的方法是，在得到相应网格后，对每一网格依据该网格所代表的节点的色值（由变量c控制），来定义这一网格的颜色．若不输入c，则默认为c=z．

我们看下面的例子：

%绘制地球表面的气温分布示意图．

[a,b,c]=sphere(40);

t=abs(c);    %求绝对值

surf(a,b,c,t);

axis equal

colormap('hot')

我们可以得到图形如下：

图4.2.3.2  等温线示意图

 

 

§4.2.4图形的控制与修饰

(1) 坐标轴的控制函数axis，调用格式如下：

axis([xmin,xmax,ymin,ymax,zmin,zmax])

用此命令可以控制坐标轴的范围．

与axis相关的几条常用命令还有：

axis auto                自动模式，使得图形的坐标范围满足图中一切图元素

axis equal                     严格控制各坐标的分度使其相等

axis square             使绘图区为正方形

axis on                  恢复对坐标轴的一切设置

axis off                 取消对坐标轴的一切设置

axis manual           以当前的坐标限制图形的绘制

 

（2）grid on  在图形中绘制坐标网格．

grid off  取消坐标网格．

 

（3）xlabel,ylabel, zlabel分别为x轴, y轴, z轴添加标注．title为图形添加标题．

以上函数的调用格式大同小异，我们以xlabel为例进行介绍：

xlabel('标注文本'，'属性1'，'属性值1'，'属性2'，'属性值2'，…)

这里的属性是标注文本的属性，包括字体大小、字体名、字体粗细等．

 

例如：

[x, y]=meshgrid(-4:.2:4);

R=sqrt(x.^2+y.^2);

z=-cos(R);

mesh(x,y,z)

xlabel('x\in[-4,4]','fontweight','bold');

ylabel('y\in[-4,4]','fontweight','bold');

zlabel('z=-cos(sqrt(x^2+y^2))','fontweight','bold');

title('旋转曲面','fontsize',15,'fontweight','bold','fontname','隶书');

 

 

图4.2.4.1 添加标注

 

以上各种绘图方法的详细用法，请看联机信息．

§4.3 统计回归图

对平面上n个点：

在平面直线族{ 为实数}中寻求一条直线 ，使得散点到与散点相对应的在直线上的点之间的纵坐标的误差的平方和最小，用微积分的方法可得：

 

所求得的这条直线： 称为回归直线．

例：已知如下点列，求其回归直线，并计算最小误差平方和．

x	0.1	0.11	.12	.13	.14	.15	.16	.17	.18	.2	.21	.23
y	42	43.5	45	45.5	45	47.5	49	53	50	55	55	60

参考的程序如下：

x=[0.10.11 .12 .13 .14 .15 .16 .17 .18 .2.21 .23];

y=[4243.5 45 45.5 45 47.5 49 53 50 55 55 60];

n=length(x);

xb=mean(x);  

yb=mean(y);

x2b=sum(x.^2)/n;

xyb=x*y'/n;

b=(xb*yb-xyb)/(xb^2-x2b);

a=yb-b*xb;

y1=a+b.*x;

plot(x,y,'*',x,y1);

serror=sum((y-y1).^2)

 

 

图4.3.1回归直线

第五节　Matlab编程

§5.1关系运算

§5.1.1 比较运算

比较两个同阶矩阵有下面六种相关操作符：

相关操作符
<	小于
<=	小于等于
>	大于
>=	大于等于
==	等于
~=	不等于

表5.1.1.1相关操作符

比较两个元素的大小，结果是“1”表明为真，结果是“0”表明为假．

例如2+2~=4

结果是“0”，表明为假．

例如一个6阶魔术方阵，矩阵元素计算满足各种条件：

A=magic(6)

ans =

    35    1     6    26   19    24

     3   32     7    21   23    25

    31    9     2    22   27    20

     8   28    33    17   10    15

    30    5    34    12   14    16

     4   36    29    13   18    11

阶数为n的魔术方阵，即n×n矩阵，是由1～n²的整数组成(n=6)．仔细观察这个矩阵，我们会发现任何行和、任何列和都相等．另外，每个3×3子行列式的对角线元素和，都可被3整除．为了显示这一特性，键入：

p=(rem(A,3)==0)

p =

     0    0     1     0    0     1

     1    0     0     1    0     0

     0    1     0     0    1     0

     0    0     1     0    0     1

     1    0     0     1    0     0

     0    1     0     0    1     0

为了再仔细地观察这个模式，可以用format+格式画出矩阵的压缩格式．此格式用“+”代表正元素，“-”代表负元素，空格代表0．

format+

p =

  +  +

+  + 

 +  +

  +  +

+  + 

 +  +

find 函数在关系运算中很有用，它可以在0-1矩阵中找非零元素的下标．

若y是一个向量，例如：y=[1   3   2   4   3.5   2.9]，则find(y<3.0)，将指出y的分量在哪些位置上小于3.0．

ans=     1     3    6

即：向量y的第1、3、6位置上的元素小于3.0．

当输入x==NaN时结果为NaN，因为按照IEEE算法规定任何具有NaN的操作，结果都是NaN．调试NaN很有用，例如测试x，输入isnan(x)函数，如果x元素是不定值则得1，否则得0．isfinite(x)更有用，如-¥<x<¥时则得1．

§5.1.2逻辑运算

&	与
|	或
~	非

表5.1.2.1逻辑运算符

 

“＆”和“|”操作符可比较两个标量或两个同阶矩阵．对于矩阵来说必须符合规则，如果A和B都是0-1矩阵，则A＆B或A|B也都是0-1矩阵，这个0-1矩阵的元素是A和B对应元素之间逻辑运算的结果，逻辑操作符认定任何非零元素都为真，给出“1”，任何零元素都为假，给出“0”．

非(或逻辑非)是一元操作符，即~A：当A是非零时结果为“0”；当A为“0”时，结果为“1”．因此下列两种表示：

p |(~p)             结果为1．

p ＆ (~p)          结果为0．

any和all函数在连接操作时很有用，设x是0-1向量，如果x中任意有一元素非零时，any(x)返回“1”，否则返回“0”；all(x)函数当x的所有元素非零时，返回“1”，否则也返回“0”．这些函数在if语句中经常被用到．如：

ifall(A<5)

do something

end

§5.2控制流

Matlab与其它计算机语言一样，也有控制流语句．控制流语句可使原本简单地在命令行中运行的一系列命令或函数，组合成为一个整体——程序，从而提高工作效率．

§5.2.1 for 循环

Matlab与其它计算机语言一样有do或for循环，完成一个语句或一组语句在一定时间内反复运行的功能．例如：

for  i = 1:n ,  x( i )=0,  end

x的第一个元素赋0值，如果n＜1，结构上合法，但内部语句不运行，如果x不存在或比n元素小，额外的空间将会自动分配．

多重循环写成锯齿形是为了增加可读性．例如：

m=9;n=9;

for i= 1:m

for j=1:n

               A( i, j ) = 1/( i + j - 1);

end

end

A

程序的说明：

（1）事实上，上述程序给出了Hilbert矩阵的构造过程，可参见函数hilb(n)．

（2）语句内部使用分号，表示计算过程不输出中间结果．

（3）循环后的A命令表示显示矩阵A的结果．

（4）每个for语句必须以end语句结束，否则是错误的．

for循环的通用形式为：

for v=expression

statements

end

其中expression表达式是一个矩阵，因为Matlab中都是矩阵，矩阵的列被一个接一个的赋值到变量v，然后statements语句运行．

通常expression是一些m:n或m:k:n仅有一行的矩阵，并且它的列是个简单的标量．但如注意到expression可以为矩阵，即v可以为向量，对某些问题的处理将大大简化．

§5.2.2  while 循环

Matlab中的while 循环语句为一个语句或一组语句在一个逻辑条件的控制下重复未知的次数．

它的一般形式为：

while  expression

statements

end

当expression的所有运算为非零值时，statements语句组将被执行．如果判断条件是向量或矩阵的话，可能需要all或any函数作为判断条件．

例如计算expm(A)，在A并不是太大时，直接计算expm(A)是可行的．

expm(A)=I+A+A^2/2!+A^3/3!+…   注意：这里的I表示单位矩阵．

程序为：

E = 0*A; F = E + eye(size(E)); N = 1;

while norm(F,1) > 0,

F = A*F/N;

E = E + F;

N = N + 1;

end

§5.2.3 if和break语句

下面介绍if语句的二个例子．

(1) 一个计算如何被分成三个部分，用符号校验：

if n<0

A=negative(n)

elseif mod(n,2)==0

A=even(n)

else

A=odd(n)

end

其中的三个函数negative(n)、even(n)、odd(n)是自编的输出函数．参见下面的函数文件．

 

(2) 这个例子涉及数论中一个很有趣的问题，取任何的正整数，如果是偶数，用2除；如果是奇数，用3乘，并加上1，反复这个过程，直到你的整数成为1．这个极有趣不可解的问题是：有使这个过程不中止的整数吗?

%classic "3n+1"problem from number theory

while 1

n=input('Enter n, negative quits:');

if n<=0 break,end

        while n>1

               if rem(n,2) == 0  %% 是连续的2个等号

                      n=n/2

               else

                      n=3*n+1

               end;

        end

end

这个过程能永远进行吗?

程序的说明：

（1）本程序用到了if语句与while语句，过程比较复杂；

（2）使用input函数，可使程序在执行过程中，从键盘输入一个数（矩阵）；

（3）break语句提供了程序跳出死循环的途径．

§5.3 M文件、命令文件及函数文件

§5.3.1 M文件

Matlab通常使用命令驱动方式，当单行命令输入时，Matlab立即处理并显示结果，同时将运行说明或命令存入文件．

Matlab语句的磁盘文件称作M文件，因为这些文件名的未尾是.M形式，例如一个文件名为bessel.m，提供bessel函数语句．

一个M文件包含一系列的Matlab语句，一个M文件可以循环地调用它自己．

M文件有两种类型：

第一类型的M文件称为命令文件，它是一系列命令、语句的简单组合．

第二类型的M文件称为函数文件，它提供了Matlab的外部函数．用户为解决一个特定问题而编写的大量的外部函数可放在Matlab工具箱中，这样的一组外部函数形成一个专用的软件包．

这两种形式的M文件，无论是命令文件，还是函数文件，都是普通的ASCII文本文件，可选择编辑或字处理文件来建立．

§5.3.2命令文件

当一个命令文件被调用时，Matlab运行文件中出现的命令而不是交互地等待键盘输入，命令文件的语句在工作空间中运算全局数据，对于进行分析解决问题及做设计中所需的一长串繁杂的命令和解释是很有用的．

例如：一个自编的命令文件fibo.m，用于计算Fibonnaci数列

% An M-file to calculate Fibonnaci numbers

f=[1, 1 ]; i = 1;

while f(i)+f(i+1)<1000

f(i+2)=f(i)+f(i+1);

i=i+1;

end

plot(f)

在Matlab命令窗口中键入fibo命令，并回车执行，将计算出所有小于1000的Fibonnaci数，并绘出图形．

要注意的是：文件执行后，f和i变量仍然留在工作空间．

§5.3.3函数文件

如果M文件的第一行包含function，这个文件就是函数文件，它与命令文件不同，所定义变量和运算都在文件内部，而不在工作空间．函数被调用完毕后，所定义变量和运算将全部释放．函数文件对扩展Matlab函数非常有用．

例如：一个自编的函数文件mean.m，用于求向量的（或矩阵按列的）平均值

function  y=mean(x)

% MEAN Average or mean value，For Vectors,

% MEAN (x) returns the mean value

% For matrix MEAN (x) is a row vector

% containing the mean value of each column

[m,n]=size(x);

if m==1

m=n;

end

y=sum(x)/m;

磁盘文件中定义的新函数称为mean函数，它与Matlab函数一样使用，例如z为从1到99的实数向量：

z=1:99;

计算均值：mean(z)

ans=

50

mean.m程序的说明：

（1）第一行的内容：函数名，输入变量，输出变量，没有这行这个文件就是命令文件，而不是函数文件．

（2）％：表明％右边的行是说明性的内容注释．前一小部分行来确定M文件的注释，并在键入helpmean后显示出来．显示内容为连续的若干个％右边的文字．

（3）变量m,n和y是mean的局部变量，在mean运行结束后，它们将不在工作空间z中存在．如果在调用函数之前有同名变量，先前存在的变量及其当前值将不会改变．

 

再例如：一个计算标准差的函数文件stat.m

function [mean,stdev]=stat(x)

[m,n]=size(x);

if m==1

m=n

end

mean=sum(x)/m;

stdev=sqrt(sum(x.^2)/m-mean.^2);

stat表明返回多输出变量是可能的．

 

又如：使用多输入变量计算矩阵秩函数

function r=rank(x,tol)

% rank of a matrix

s=svd(x);

if(nargin==1)

tol=max(size(x))*s(1)*eps;

end

r=sum(s>tol);

这个变量说明利用永久变量nargin确定输入变量的个数，变量nargout虽然这里没有使用，但它包含有输出变量的个数．

一些有用的说明：

当M函数文件第一次在Matlab运行时，它被编译并放入内存，以后使用时不用重新编译即可得到．

what命令：显示磁盘当前目录中的M文件，

dir命令：列出所有文件．

一般而言，输入一个名字到Matlab，例如键入whoopie命令，Matlab用以下步骤解释：

(1) 看whoopie是否为变量．

(2) 检验whoopie是否为在线函数．

(3) 检验whoopie文件的当前目录．

(4) 将whoopie看成Matlab的PATH中的一个文件，在Matlab PATH目录中搜索．

如果whoopie存在，Matlab首先将其作为变量而不是作为函数．

§5.4 字符串、输入及输出

§5.4.1 echo、input、pause、keyboard

一般来说，当一个M文件运行时，文件的命令不在屏幕上显示，而echo命令则使M文件运行时，命令在屏幕上显示，这对于调试、演示相当有用．

input功能：输入Input('How many apples')给用户一个提示串，等待，然后显示用户通过键盘输入的大量表达式．可以用input命令建立驱动M文件的菜单．

与input功能相同，但功能更强的keyboard命令将计算机作为一个命令文件来调用，放入M文件中，此特性对调试或正在运行期间修改变量很有用．

pause命令：使用户暂停运行一个程序，当再按任一键时恢复执行，pause(n)等待n秒钟后再继续执行．

§5.4.2串和宏串

字符串用单个引号输入到Matlab中，例如：

s='Hello'

结果显示为：

s =

Hello

字符存在向量中，每个元素就是一个字符，如：

size(s)

ans =

1       5

表明S为一个1×5的矩阵，有五个元素．字符以ASCII值存入，abs函数或double函数将显示以下值(即Hello的ASCII值)

abs(s)

ans =

   72   101  108   108   111

getstr函数，使向量作为字符显示，而不显示ASCII值．

disp可在变量中显示字符．sprintf, num2str和int2str可以将数字转换成串．

字符变量通过括号连成大串．例如：

s='hello';

s=[s,'world']

s =

hello world

eval是与字符变量—起工作的函数，执行简单字符宏调用．eval( t )执行包含在t内的字符．如果t是任何Matlab表达式或语句的源字符，则字符串被解释执行．例如：

t='eye(2)', eval(t)

结果为：

ans=

1     0

0     1

又例如，给矩阵元素赋值

t='1/(i+j-1)';

for i=1:n

for j=1:n

               a( i,j)=eval(t);

end

end

这儿有一个例子，介绍如何一起使用eval与load命令，装入十个具有顺序文件名的文件中的数据：

fname='mydata';

for i = 1:10

eval([ 'load ', fname, int2str( i )])

end

§5.4.3外部程序

Matlab与外部独立程序的通讯方式可以是多种多样的，下面介绍其中的一个办法：

(1) Matlab中将变量存入磁盘

(2) 运行外部程序（读数据文件，进行处理），将结果写到磁盘上

(3) 将处理后的文件装回到工作空间中

例如：用外部程序gareqn找garfield方程的结果：

function y=garfield(a,b,q,r)

save gardata a,b,q,r

! gareqn

load gardata

使用FORTRAN或其它语言写gareqn程序，使其可以读gardata.mat，进行处理，将结果存入文件中．

这个程序可将计算机的“连接码”提供给Matlab，在许多系统中它将新的目标码连接到程序中比物理联接要方便得多．

§5.4.4输入输出数据

可使用各种方法将其它程序和外部世界的数据送入Matlab，同样可把Matlab数据输送到外部世界，使你的程序以Matlab使用的文件形式直接计算数据．

最好的方法取决于多少数据，数据是否可读，什么形式等：

(1) 清晰的元素表输入：

如果你有少量数据，比如说小于10~15个元素，使用方括号[]输入．

(2) 使用文本编辑建立命令文件，将数据列为清晰的元素表输入．如果数据不是可读形式，又不得不以一种方法键入，可以重复运行M文件，重复修改数据．

(3) 如果数据以ASCII形式存贮，并有固定长度，行尾有回车符，各数间有空格的文件称为flatfile(ASCII的flat file可由普通文本编辑来编辑)，flat file通过load命令直接读进Matlab，结果存入名为文件名的变量中去．

(4) 将数据文件译成Matlab文件形式，使用load命令，translate程序由Matlab中的应用程序库支持，translate程序将ASCII文件、二进制文件、FORTRAN非格式文件和DIF文件转换为Matlab使用的特定的MAT文件，当磁盘文件中存有大量数据时，这个方法输入最好．

Matlab数据输出到外部世界的方法：

(1) 小矩阵时：使用diary命令建立日志文件，在文件中列出变量，用文本编辑处理日志文件，日志的输出包括运行中的Matlab命令．

(2) 使用save命令存入变量，退出Matlab，用translate程序将MAT文件转换成任一种其它文件形式．

第六节　Matlab符号运算

 

Matlab本身并没有符号计算功能，1993年通过购买Maple的使用权后，开始具备符号运算的功能．符号运算的类型很多，几乎涉及数学的所有分支．

§6.1 Matlab符号运算的工作流程

§6.1.1工作过程

 

§6.1.2 核心工具

sym函数  VS    syms 语句

sym函数：构造符号变量和表达式:  a=sym('a')

（Construct symbolic numbers,variables and objects）

syms语句：构造符号对象的简捷方式（Short-cut forconstructing symbolic objects）

§6.1.3 符号变量确定原则

（1）除了i 和j之外，字母位置最接近x的字母；若距离相等，则取ASCII码大的；

（2）若没有除了i 与j以外的字母，则视x为默认的符号变量；

（3）可利用函数findsym(string,N)来询问在众多符号中，哪N个为符号变量．例如：键入findsym(3*a*b+y^2,1)，即可得到答案y．更多的例子见下表：

 

符号表达式	默认符号变量
a*x^2+b*x+c	x
1/(4+cos(t))	t
4*x/y	x
2*a+b	b
2*i	x

§6.2  Matlab的六大常见符号运算

§6.2.1 因式分解

symsx

f=x^6+1;

s=factor(f)

结果为：

s=(x^2+1)*(x^4-x^2+1)

§6.2.2 计算极限

求极限：

 

  syms h n x

  L=limit('(log(x+h)-log(x))/h',h,0)  %%单引号可省略掉

 M=limit('(1-x/n)^n',n,inf) 

结果为：

 L =1/x

 M =exp(-x) 

§6.2.3 计算导数

．

  syms a x;   y=sin(a*x);

 A=diff(y,x)

 B=diff(y,a) 

 C=diff(y,x,2)

结果为：

  A = cos(a*x)*a

 B = cos(a*x)*x

 C = -sin(a*x)*a^2  

§6.2.4 计算不定积分、定积分、反常积分

，

， 

．

    syms x

    f=(x^2+1)/(x^2-2*x+2)^2;

    g=cos(x)/(sin(x)+cos(x));

    h=exp(-x^2);

    I=int(f)

    J=int(g,0,pi/2)

K=int(h,0,inf) 

结果为：

  I =1/4*(2*x-6)/(x^2-2*x+2)+3/2*atan(x-1)

  J =1/4*pi

  K =1/2*pi^(1/2) 

 

§6.2.5 符号求和

求级数   的和S, 以及前十项的部分和S1．

    syms n

    S=symsum(1/n^2, 1, inf)

S1=symsum(1/n^2,1,10) 

结果为：

  S =1/6*pi^2

  S1 =1968329/1270080

重要说明：当求函数项级数   的和S2时，可用命令：

syms n x

S2=symsum(x/n^2, n, 1, inf)   

S2 =1/6*x*pi^2

 

两点说明：

（1）注意观察S2与S1的细微区别！

（2）当通项公式的Matlab表达式较长时，表达式要加上单引号．后面的练习中会遇到此问题．

§6.2.6 解代数方程和常微分方程

利用符号表达式解代数方程所需要的函数为solve(f)，即解符号方程式f．

例如：求一元二次方程a*x^2+b*x+c=0的根．

f=sym('a*x^2+b*x+c')  或 f='a*x^2+b*x+c'

solve(f)

ans=

[1/2/a*(-b+(b^2-4*c*a)^(1/2))]

[1/2/a*(-b-(b^2-4*c*a)^(1/2))]

solve(f,a)

ans=

-(b*x+c)/x^2

 

利用符号表达式可求解微分方程的解析解，所需要的函数为dsolve(f)，使用格式：

dsolve('equation1', ' equation2', …)

其中：equation为方程或条件．写方程或条件时，用 Dy 表示y 关于自变量的一阶导数，用D2y 表示 y 关于自变量的二阶导数，依此类推．

1． 求微分方程 的通解．

syms x y   %定义x,y为符号

dsolve('Dy=x',  'x')

ans =

1/2*x^2+C1

试比较：

若写成：

syms x y   %定义x,y为符号

dsolve('Dy=x')

结果将是什么？是否正确？为什么？

 

2． 求微分方程 的特解．

syms x y

dsolve('D2y=x+Dy','y(0)=1', 'Dy(0)=0', 'x')

ans =

-1/2*x^2+exp(x)-x

试比较：

若写成：

syms x y

dsolve('D2y=x+Dy','y(0)=1', 'Dy(0)=0')

 

结果将是什么？是否正确？为什么？

 

3． 求微分方程组 的通解．

syms x y

[x,y]=dsolve('Dx=y+x,Dy=2*x') 

x =

1/3*C1*exp(-t)+2/3*C1*exp(2*t)+1/3*C2*exp(2*t)-1/3*C2*exp(-t)

y =

2/3*C1*exp(2*t)-2/3*C1*exp(-t)+2/3*C2*exp(-t)+1/3*C2*exp(2*t)

 

试比较：

若写成：

（1）dsolve('Dx=y+x, Dy=2*x')

结果将是：

ans =

    x: [1x1 sym]

y: [1x1 sym]

试解释此结果的含义．

 

若写成：

（2）[x,y]=dsolve('Dx=y+x, Dy=2x')

结果将是：

x=

exp(t)*C1+C2*exp(t)-C2-2-2*t

y =

             C2+2*t

是否正确？为什么？
综合练习

 

1．按顺序进行如下的操作：

（1）产生一个5阶魔术方阵A；并计算A'与A^-1(即inv(A))；

（2）求A的特征值；

（3）计算A的各列的总和与平均值；

（4）计算A的各行的总和与平均值；

（5）若b=[1 2 3 45] '，求方程组 Ax=b的解；

（6）验证你的结论的正确性．

2．产生行向量S =[1.0, 1.2, 1.4, …,  20]，并计算S * S' 与  S' * S，你有何“发现”？

3．设A= ；B= ；求C=A * B – B * A，你有何“发现”？

4．若设矩阵A= ；B= ；求C=A * B – B * A，你又有何“发现”？

5．如何建立如下的矩阵（命令方式和程序方式）？

 

（１） ；　　（２） ；

 

 

（３） ；　　（４）

（５）                （6）  

（７）

6．绘制下列曲线的图形（散点图与折线图）：

7．绘制下列曲面的图形： （提示：曲面由两部分构成）

8．在同一个图形上作下列两个函数的图象：

（1） ；    （2）

9．假如你有一组实测数据，例如：

x=[53  56     60    67.5    75    90  110];

y=[109 120.5  130   141.1  157.5  180  185];

求其回归直线，画回归直线图形并计算最小误差平方和．

10．假如你有一组实测数据，例如：

x=[75 86  95  108  112  116 135  151  155 160  163  167 171  178 185];

y=[10 12  15  17   20   22   35   41   48  50   51   54  59   66   75];

求其回归直线，画回归直线图形并计算最小误差平方和．

11．随机产生500个0到100的整数FS作为学生的考试分数．

（1）       画出FS的简单直方图；

（2）       画出每个分数段（0～10、10～20、…，90～100）的统计频数直方图；

12．求下列各结果：

（1）用Matlab因式分解： ．

（2）用Matlab求极限： ．

（3）用Matlab求积分： ．

（4）用Matlab求幂级数： 的和函数（化简结果）．

13．非线性回归尝试

下表是到1994年的游泳世界纪录，试估计时间y与距离x的关系．

距离x（米）	50	100	200	400	800	1500
时间y（秒）	21.81	48.42	106.69	225	466.60	863.48

 

说明：用线性回归方法将得到： ，但当 时， ，这是非常荒唐的结果！显然，一个基本要求是当 时 ．试尝试使用非线性回归模型: .

 

14. (三维)符号作图尝试

命令	作用	解释
ezplot3	3-D parametric curve plotter	3D参数曲线图形
ezcontour	use contour plotter	等高线图
ezcontourf	filled contour plotter	填充等高线图
ezmesh	3-D mesh plotter	3D mesh曲面图形
ezmeshc	combination mesh/contour plotter	mesh曲面/等高线图
ezsurf	3-D colored surface plotter	3D surf曲面图形
ezsurfc	combination surf/contour plotter	surf曲面/等高线图

 

请尝试以下的命令：

ezplot3('sin(t)', ' cos(t)','t', [0,6*pi])

ezcontour('x*exp(-x^2 - y^2)')

ezcontourf('x*exp(-x^2 - y^2)')

ezmesh('(s-sin(s))*cos(t)','(1-cos(s))*sin(t)','s',[-2*pi,2*pi])

ezmeshc('(s-sin(s))*cos(t)','(1-cos(s))*sin(t)','s',[-2*pi,2*pi])

ezsurf('x*exp(-x^2 - y^2)')

ezsurfc('x*exp(-x^2 - y^2)')