[工学]matlab教案.doc

1 目录目录 第 1 章 绪论.2 1.1 科学计算简介.2 1.2 Matlab 概述2 第 2 章 Matlab 运算基础.5 2.1 变量与赋值.5 2.2 矩阵.6 2.3 表达式.9 2.4 数学函数.9 第 3 章 Matlab 程序设计.10 3.1 m 文件10 3.2 数据的输入输出.10 3.3 程序结构.11 3.4 函数文件.12 第 4 章 图形与声音.14 4.1 二维图形.14 4.2 三维图形.14 4.3 图形窗口控制.15 4.4 图形控制.16 4.5 动画.16 4.6 声音.16 第 5 章 线性代数.18 5.1 矩阵.18 5.2 向量空间.19 5.3 线性方程组.21 5.4 特征值与特征向量.25 5.5 二次型及其标准型*27 第 6 章 数据处理与多项式.30 6.1 基本统计处理.30 6.2 多项式.31 6.3 数据插值.32 6.4 曲线拟合.35 6.5 离散傅立叶变换.36 第 7 章 数值积分与微分方程.38 7.1 数值积分.38 7.2 数值微分.38 7.3 常微分方程的数值解.41 7.4 非线性方程(组)求解45 7.5 函数优化.47 第 8 章 符号计算.48 8.1 符号计算基础.48 8.2 微积分.49 8.3 线性代数.51 8.4 方程求解.53 2 第 1 章 绪论 1.1 科学计算简介 科学计算，即对科学和工程中的数学问题进行数值计算。数值计算的过程主要包括建立数学模型、 建立求解的计算方法、计算机实现三个阶段。数值计算的特点是计算方法比较复杂，方法种类多种多样， 如数值微分、数值积分、常/偏微分方程、线性代数方程、有限元等。数值计算所关心的焦点是计算精度 （误差影响） 。 科学计算可分为两类：一类是纯数值的计算，例如求函数的值，方程的数值解;另一类计算是符号计 算，又称代数运算，这是一种智能化的计算，处理的是符号。符号可以代表整数，有理数，实数和复数， 也可以代表多项式，函数，还可以代表数学结构如集合，群的表示等等。我们在数学的教学和研究中用 笔和纸进行的数学运算多为符号运算。 主要的数学软件有：Mathematica、MATLAB、Maple 和 MathCAD。 尽管计算机代数系统在代替人进行繁琐的符号运算上有着无比的优越性，但计算机毕竟是机器，它 只能执行人们给它的指令。此外，虽然计算机代数系统包含大量的数学知识，但这仅仅是数学的一小部 分，目前有许多数学领域计算机代数系统还未能涉及。 1.2 Matlab 概述 1980 年前后，当时的新墨西哥大学计算机系主任 Cleve Moler 教授在讲授线性代数课程时，发现了 用其他高级语言编程极为不便，便构思并开发了 MATLAB (MATrix LABoratory，即矩阵实验室), 这一软 件利用了当时数值线性代数领域最高水平的 EISPACK 和 LINPACK 两大软件包中可靠的子程序，用 Fortran 语言编写了集命令翻译、科学计算于一身的一套交互式软件系统。只要给出一条命令，立即就可 以得出该命令的结果，而无需像 C 和 Fortran 语言那样，首先编写源程序，然后对之进行编译、连接，最 终形成可执行文件。这无疑会给使用者带来了极大的方便。此后，Cleve Moler 等人成立了 MathWorks 的 公司，用 C 语言作了完全的改写，于 1984 年推出了第一个 MATLAB 的商业版本。 其后又增添了丰富 多彩的图形图像处理、多媒体功能、符号运算和它与其他流行软件的接口功能，使得 MATLAB 的功能 越来越强大，应用范围越来越广。到目前为止其最高版本 7.0 版已经推出。 MATLAB 具有强大的数学运算能力、方便实用的绘图功能及语言的高度集成性，它在其他科学与工 程领域的应用也是越来越广，并且有着更广阔的应用前景和无穷无尽的潜能。 MATLAB 可以将使用者 从繁琐、无谓的底层编程中解放出来，把有限的宝贵时间更多地花在解决问题中，这样无疑会提高工作 效率。目前，MATLAB 已经成为国际上最流行的科学与工程计算的软件工具，现在的 MATLAB 已经不 仅仅是一个“矩阵实验室”了，它已经成为了一种具有广泛应用前景的全新的计算机高级编程语言了， 在国内外高校和研究部门正扮演着重要的角色。可以预见，在科学运算、自动控制与科学绘图领域 MATLAB 语言将长期保持其独一无二的地位。 一、 MATLAB 的主要特点 有高性能数值计算的高级算法特别适合矩阵代数领域； 有大量事先定义的数学函数并且有很强的用户自定义函数的能力； 有强大的绘图功能以及具有教育、科学和艺术学的图解和可视化的二维、三维图； 适合个人应用的强有力的面向矩阵(向量)的高级程序设计语言； 3 与其它语言编写的程序结合和输入输出格式化数据的能力； 有在多个应用领域解决难题的工具箱。 本教程基于 MATLAB6.5 版。 二、 操作环境 matlab 主界面：菜单、工具栏、工作区（命令窗口） 、workplace brower、commmand history 等。 Figtrue 窗口： m 文件编辑窗口： 有些工具箱也提供了操作界面，可以通过命令或选择 start 的 toolboxes 菜单进入。 三、 引例 2 + 6 4 ans = 4 ans/2 ans = 2 也可以定义变量 a = 5 a = 5 b = 6 b = 6 4 c = b/a c = 1.2000 可以使用内建的变量和函数 pi ans = 3.1416 sin(ans/4) ans = 0.7071 四、 帮助 在 MATLAB 系统中相关的线上(on-line)求助方式有： 1. 利用 help 或 helpwin 指令：在命令窗口键入 help ，例如 help sqrt, help(sqrt), helpwin sqrt。 2. 利用 lookfor 指令：键入的 lookfor (key-word)，列出所有相关的题材，例如 lookfor cosine, lookfor sine。 3. 利用指令视窗的功能选单中的 Help，从中选取 Table of Contents（目录）或是 Index（索引） 。 4. 利用 doc 指令：doc；或 doc function 5 第 2 章 Matlab 运算基础 2.1 变量与赋值 一、 变量 matlab 变量主要是矩阵，矩阵的元素可以是数值或字符串。其中数值可以是实数或复数，虚数单位用 i 或 j 表示。如 a=3+9i（最好不要写成 3+9*i，因为 i、j 是可以作为普通变量被赋值的） 特别地，如果矩阵只有一行或一列，就称为向量(vector),，如果只有一行一列，就称为标量(scalar)。 其他的变量类型有：多维数组、结构、对象等。 变量的命名：matlab 变量无需声明即可直接使用，变量名以字母开头，接字母、数字或下划线，最多 63 个字符（多余的被忽略） ，区分大小写。 系统变量： ans：value of an expression when that expression is not assigned to a variable eps： fl oating point precision 2.2204e-016 (2-52) pi： , (3.141492 . . .) realmax： 最大正浮点数 1.7977e+308 (21024) realmin：最小正浮点数 2.2251e-308 (2-1022) Inf： ，正无穷大, a number larger than realmax, the result of evaluating 1/0. NaN： not a number, the result of evaluating 0/0 i, j：虚数单位 不要对系统变量赋值，否则使用内建函数时会出错。例外：i, j 常作为循环变量使用。 二、 赋值语句 1）变量=表达式 2）直接输入表达式：结果被赋给系统变量 ans 注意： 语句后加分号,表示仅执行，不显示结果 续行符为（三个句点） 注释：以%开头 例子：p14，计算表达式的值，并将结果赋给变量 a。 2sin(85 ) 153i a=2*sin(85*pi/180)/(1+sqrt(5)+3i) a = 0.3311 - 0.3070i 变量的查看：who, whos。who 列出工作空间中的所有变量名，whos 同时给出变量的维数和性质。变 量也列在 workplace browser 中。 数据的格式：format 格式符，格式符参见 p18。如 format short。也可以通过 file-prefrences 菜单， 在 prefrences 对话框的 command window 项中设置。 6 变量的清除：clear，清除工作空间中的所有变量，等同 edit-clear workplace 菜单。 清除命令窗口：clc，等同 edit-clear conmmand window 菜单。 2.2 矩阵 一、 矩阵的建立 1） 直接输入 将矩阵元素用方括号括起来，同一行各元素间用空格或逗号隔开；各行之间用分号或回车隔开。 例：A=1 2 3 ; 4 5 6; 7 8 9 或 A=1, 2, 3 4, 5, 6 7, 8, 9 2）利用函数建立 很多函数都会返回一个或几个矩阵作为计算结果。以下函数则专门用于建立特殊矩阵（向量） 。 linspace：建立等差数列，x = linspace(startValue,endValue,nelements)，若省略 nelements，则为 100。 例： v = linspace(0,9,4) v = 0 3 6 9 logspace：建立等比数列。 y = logspace(a,b,n)：在 10a 与 10b 之间产生 n 个数，若省略 n，则为 50 y = logspace(a,pi,n)：在 10a 与 pi 之间产生 n 个数，若省略 n，则为 50。 例： y = logspace(0,2,5) y = 1.0000 3.1623 10.0000 31.6228 100.0000 特殊矩阵的建立（参见 p83）： 空矩阵：，产生行数、列数均为 0 的矩阵。 对角线矩阵：diag(V)用向量 V 建立一个对角线矩阵，diag(A)从矩阵 A 中提取主对角线向量。 diag(A,k)提取主对角线上方第 k 条对角线。 单位矩阵：eye(n) 幺矩阵：ones(m), ones(m,n)产生元素值均为 1 的矩阵 零矩阵：zeros(m), zeros(m,n)产生元素值均为 0 的矩阵 随机矩阵：元素值为随机值。 rand(n)、rand(m.n)：分别建立 n*n、m*n 阶矩阵，矩阵元素为 01 之间均匀分布的数；randn(n)、randn(m,n)：01 之间正态分布。 3）通过 mat 文件保存和建立矩阵 save filename var1 var2 .：把当前工作空间 (workpalce) 中的指定变量保存到文件 filename（扩展名 为.mat）中，如果不指定变量在，则保存所有变量。 load filename var1 var2 .：把文件 filename.mat 中的变量载入到当前工作空间中。如果不指定变量在， 则载入所有变量。 也通过 workplace browser 的工具栏按钮可以完成以上操作。 4) 通过 workplace browser 编辑矩阵 双击 workplace browser 中的变量，或选中变量后单击工具栏的 open 按钮，可以打开 Array Editor 编 辑矩阵。 7 二、 矩阵的基本操作 1、矩阵元素的操作：如 A(3,2)=200，如果给出的行数或列数大于原来矩阵的范围，则会自动扩展原来 的矩阵，并将扩展后未赋值的元素置 0。 2、冒号操作：冒号是一个重要的运算符，利用冒号可以产生向量、拆分矩阵。 产生向量的一般格式是：e1:e2:e3。其中 e1 为初始值；e2 为步长，若省略则为 1；e3 为终止值。 例： t=0:2:10 t = 0 2 4 6 8 10 拆分矩阵：A(i, :) 提取 A 矩阵的第 i 行元素。A(:, j) 提取 A 矩阵的第 j 列元素 A(i:i+m, k:k+m)：提取 A 矩阵的第 ii+m 行元素、第 kk+m 列元素。 例：A=1 2 3 4 5; 6 7 8 9 10; 11 12 13 14 15; 16 17 18 19 20; A(3, :) ans = 11 12 13 14 15 A(2:3, 4:5) ans = 9 10 1415 A(2:3, 1:2:5) ans = 6 8 10 11 13 15 ans(:) ans = 6 11 8 13 10 15 3、矩阵的合并 大矩阵可由方括号中的小矩阵建立起来。 例：A=1 2 3; 4 5 6; 7 8 9; C=A eye(size(A) C = 1 2 3 1 0 0 4 5 6 0 1 0 8 7 8 9 0 0 1 4、删除矩阵元素 利用空矩阵，可将矩阵的某些行和列删除。如 A(:,2,4)=将删除 A 的 2 列和第 4 列元素。 例：A=1 2 3 4 5; 6 7 8 9 10; 11 12 13 14 15; A(:,2,4)= %删除 A 的 2 列和第 4 列元素 A = 1 3 5 6 8 10 11 13 15 5、矩阵的加减法 运算规则：对应元素相加、减。A+B，A-B 6、数乘：k*A 7、 矩阵的乘法 运算规则：按线性代数中矩阵乘法运算进行，即放在前面的矩阵的各行元素，分别与放在后面的矩 阵的各列元素对应相乘并相加。即 n ijikkj k 1 ca *b 例：A=1 2 3; 4 5 6; B=1 2; 3 0; 7 4; A*B ans = 28 14 61 32 8、数组乘法：将矩阵的对应元素相乘。x.*y 例：x=1 2 ; 3 4 ; y=5 6; 7 8; z=x.*y z = 5 12 21 32 9、矩阵的除法 有两种除法运算：（左除）和/（右除） 。AB 相当于 inv(A)*B，B/A 相当于 B*inv(A)。对于标量来 说，/就是通常所用的除法，且 A/B 和 BA 是等价的。 矩阵的除法在讲线性方程组求解再详细讲解。 10、数组除法：将矩阵的对应元素相除。x.y 和 y./x 是等价的。 11、矩阵的乘方： 例：A=1 2; 3 4; A2 ans = 7 10 15 22 12、数组的乘方：求矩阵的每个元素的幂，. 例：A=1 2; 3 4; A.2 ans = 1 4 9 9 16 12、转置 把 A 的所有列换成相应的列，得到的矩阵称为 A 的转置矩阵。运算符是。 （对称矩阵 A=A；反对称 矩阵 A=-A） 例：A=1 2; 3 4; A ans = 1 3 2 4 2.3 表达式 算术表达式：+, -, *, /, , （乘方） 关系表达式：, =, =, =，若结果成立，关系表达式结果为 1，否则为 0。 逻辑表达式：1 5 11 A = 3 6 1 5 2. pause 函数 pause(延迟秒数)，省略秒时，需要按任意键后才继续执行 3. disp 函数 disp(输出项)，输出项可以是字符串或矩阵。 3.3 程序结构 一、 选择结构 1. if 语句 if 语句根据逻辑表达式的值，决定是否执行相应的语句。其格式为： if expression1 commands evaluated if expression1 is True elseif expression2 commands evaluated if expression2 is True elseif expression3 commands evaluated if expression3 is True else commands evaluated if no other expression is True end 如果在表达式中的所有元素为真(非零)，那么就执行 if 和 end 语言之间的commands。在表达式包 含有几个逻辑子表达式时，即使前一个子表达式决定了表达式的最后逻辑状态，仍要计算所有的子表达 式。 2. switch 语句 switch 语句根据表达式取值的不同，分别执行不同的语句。其格式为： switch expression (scalar or string) case value1 statements % Executes if expression is value1 case value2 statements % Executes if expression is value2 otherwise statements % Executes if expression does not match any case 二、 循环 1. for 语句 For 循环允许一组命令以固定的和预定的次数重复。For 循环的一般形式是： for x = array 12 commands end 在 for 和 end 语句之间的commands按数组中的每一列执行一次。在每一次迭代中，x 被指定为数组 的下一列，即在第 n 次循环中，x=array(:, n)。 2. while 语句 与 For 循环以固定次数求一组命令的值相反，While 循环以不定的次数求一组语句的值。While 循环 的一般形式是： while expression commands end 只要在表达式里的所有元素为真，就执行 while 和 end 语句之间的commands。通常，表达式的求 值给出一个标量值，但数组值也同样有效。在数组情况下，所得到数组的所有元素必须都为真。 3. continue 语句：跳过后面的语句，进入下一次迭代。 4. break 语句：跳出循环 3.4 函数文件 主函数与子函数：除内联函数外的所有函数都需要在 m 文件中定义。 每个 M 文件中可以定义一个 或多个函数，最先出现在该 m 文件中的为主函数，其余的为子函数。主函数的范围比子函数要广，可以 在 M 文件外部调用（在命令窗口或者是其他的 M 文件中） ，而子函数则只在主函数和该 M 文件的其他子 函数中可见。 私有函数：私有函数是主函数文件的一种。私有函数放在以 private 命名的子目录下，他们只是对其 父目录中的函数是可见的。因此他们可以采用与其他目录下函数相同的名字，这在当你想创建自己特定 的函数的新版本而想在另外目录保存原来版本的函数的时候很有用。 全局变量和局部变量：函数文件与命令文件在通信方面是不同的，函数内部的变量是局部变量，不 出现在 MATLAB 工作空间。但在函数文件中可以用 global 定义全局变量，全局变量的作用域是整个 matlab 工作空间。例：global a，则 a 为全局变量。全局变量应尽量少使用。 一、 函数文件格式 function 输出形参表=函数名(输入参数表) %注释说明部分 函数体 二、 函数文件的规则和属性 1. 函数名和文件名必须相同。例如，函数 fliplr 存储在名为 fliplr.m 文件中。 2. 在函数 M 文件中，到第一个非注释行为止的注释行是帮助文本。当需要帮助时，返回该文本。 例如， help fliplr 返回 fliplr 的注释部分。第一行帮助行，名为 H1 行，是由 lookfor 命令搜索的行。 3. 函数可以有零个或更多个输入参量。函数可以有零个或更多个输出参量。 4. 从函数 M 文件内可以调用命令文件。在这种情况下，命令文件查看函数工作空间，不查看 MATLAB 工作空间。 5. 函数可以递归调用。即 M 文件函数能调用它们本身。 例： p42，求 n!。 在编制要递归调用的函数时，必须确保会终止，否则 MATLAB 会陷入死循环。 13 6. 当函数 M 文件到达 M 文件终点，或者碰到返回命令 return，就结束执行和返回。return 命令提供 了一种结束一个函数的简单方法，而不必到达文件的终点。 7. MATLAB 函数 error 在命令窗口显示一个字符串，放弃函数执行，把控制权返回给键盘。这个函 数对提示函数使用不当很有用，如在以下文件片段中： if length(val)1 error( VAL must be a scalar. ) end 这里，如果变量 val 不是一个标量，error 显示消息字符串，把控制权返回给命令窗口和键盘。 三、 内联函数 内联函数用于定义简单的函数，他可以直接在命令窗口中定义。 内联函数的定义为：g = inline(expr, arg1, arg2, .)。 例：g = inline(sin(alpha*x),x,alpha); g(0.5,pi) ans = 1 四、 函数调用 格式：输出参数表=函数名输入参数表 函数调用的注意事项： 1. 函数调用时，各实参出现的顺序要与函数定义时形参的顺序一致。函数可以按少于函数文件中所 规定的输入和输出变量进行调用，但不能用多于函数 M 文件中所规定的输入和输出变量数目。 2. 当函数有一个以上输出变量时，输出变量包含在括号内。例如，V,D = eig(A)。不要把这个句法 与等号右边的V,D 相混淆。右边的V,D 是由数组 V 和 D 所组成。 3. 利用系统变量 nargin 和 nargout，可以得到调用时的输入参数和输出参数的数目，从而决定函数如 何处理。nargin 为输入参量个数；nargout 为输出参量个数。 例如，考虑 MATLAB 函数 linspace ： function y = linspace(d1, d2, n) if nargin = 2 n = 100; end y = d1+(0:n-2)*(d2-d1)/(n-1) d2 ; 如果用户只用两个输入参量调用 linspace ，例如 linspace(0,10) ，linspace 产生 100 个数据点。相反， 如果输入参量的个数是 3，例如 linspace(0,10,50)，第三个参量决定数据点的个数。 3. 函数句柄 有些函数的输入参数也是函数（如 fplot） ，这时可以通过函数句柄来传入。 函数句柄的定义方法是：在函数名前加一个。 例：fhandle = sin; fplot(fhandle ,0,2*pi) 可以使用格式：y1,y2,. = feval(fhandle,x1,.,xn)来调用函数句柄 fhandle 所指向的函数。例： x=feval(sin,pi/2) x = 1 习题：教材 p49-51：2、3、4、5、6 题 14 第 4 章 图形与声音 4.1 二维图形 一、 plot 函数 plot(x,y, 属性值)： plot(x,y1,属性值，x,y2,属性值) h=plot() 以 x 为横坐标，y 为纵坐标，按照坐标(xj ,yj)的有序排列绘制曲线。属性值用于定义图形属性，如线 型、颜色等（见 p61） ，属性值可省略。 例：x = -pi:.1:pi; y1 = sin(x); y2 = cos(x); plot(x,y1,b-,x,y2,r) 二、 fplot 函数 fplot 绘制函数曲线，它可以自适应地对函数采样。 格式：fplot(函数, 变量取值范围, 相对误差, 选项) 绘制函数的可以是标准函数、表达式，或用户自定义的函数。函数形式必须为 y=f(x)。变量取值范围 为xmin xmax，或xmin xmax ymin ymax，默认的相对误差为 2e-3。 例：fplot(sin(x),cos(x), 0, 2*pi, 1e-3, *) fplot(sin,0,2*pi) %使用函数句柄 fplot(inline(sin(x),0,2*pi) %使用内联函数 三、 特殊坐标图形 双对数坐标系：loglog(x,y)。两个坐标轴均用以 10 底的对数刻度标定。 单对数坐标系：semilogx(x,y), x 轴用以 10 底的对数刻度标定。semilogy(x,y)，y 轴用以 10 底的对数 刻度标定 极坐标：polar(theta,r)，向量 theta 的元素代表弧度参数，向量 r 代表从极点开始的长度。 四、 其他图形函数 stairs(x,y)：阶梯图形 bar(x,y)：条状图形 fill(x,y,c)：填充图形。绘制并填充二维多边图形，x 和 y 为多边形顶点向量，c 规定填充颜色。 errorbar ( x , y , e)：误差条形图。绘制向量 y 对 x 的误差条形图，误差条对称地分布在 y 的上方和下 方，长度为 e。 4.2 三维图形 一、 plot3 函数 plot3 ( x, y, z )：用(xi, yi, zi)所定义的点绘制三维曲线。向量 x、y 和 z 必须为等长度的。该命令与 15 plot 类似。 例：绘制三维螺旋线 t = 0:pi/50:10*pi; plot3(sin(t),cos(t),t) 二、 mesh 函数 用法：h=mesh(X,Y,Z,C) 功能：生成由 X，Y 和 Z 指定的网线面。若 X 与 Y 均为向量，length(X)=n，length(Y)=m，而m，n =size(Z)，则空间中的点 (X(j),Y(I),Z(i,j) 为所画曲面网线的交点，分别地，X 对应于 Z 的列，Y 对应于 Z 的行。 C 是一个矩阵，用于指定三维网格点的颜色。若省略 C，则颜色由 Z 指定。 h 为返回的图形对象句柄，可省略。 例：x,y=meshgrid(0:0.15:2*pi,0:0.15:2*pi); z=sin(y).*cos(x); mesh(x,y,z) (a) mesh 图 (b) surf 图 (c) waterfall 图 三、 surf 函数 用法：h=surf(X,Y,Z,C) 功能：绘制三维曲面图。将网格线之间的补面(patch)用颜色填充。 例：接上个例子，将 mesh(x,y,z)改为 surf(x,y,z)。 四、 waterfall 函数 用法：h = waterfall(X,Y,Z) 功能：绘制瀑布图。 例：接上个例子，将 mesh(x,y,z)改为 waterfall(x,y,z)。 4.3 图形窗口控制 figure：创建一个新图形窗口 figure(h)：激活图形窗 h close(h)：关闭图形窗 h close all：关闭所有图形窗口 hold on：设置图形保持状态 hold off：关闭图形保持 子图：subplot(m,n,p)将图形窗口分割成 m 行 n 列，并设置 p 所指定的子窗口为当前窗口。子窗口按 行由左至右，由上至下进行编号。 16 4.4 图形控制 一、 图形标记 title：Titles for 2-D and 3-D plots xlabel, ylabel, zlabel：labels for 2-D and 3-D plots text(x,y,string)：adds the string in quotes to the location specified by the point (x,y). gtext(string)：Place text on 2-D graph using mouse legend：Graph legend for lines and patches。添加图例 texlabel：Produce the TeX format from character string 二、 坐标轴设置 axis 函数设置坐标轴的，其格式较多，参见 p63。 如 axis(xmin xmax ymin ymax)，设置坐标轴的最大值和最小值。 axis off，关闭坐标轴。 三、 观察点 view (az, el) 设置观察点角度。标量 az 是方位角，即观察点在 xy 平面上投影得到的角度。el 为仰角， 表示观察点在 xy 平面上下方的角度。az 和 el 都用度来衡量。 以上内容都可以在 Figure 窗口中，通过菜单或工具栏按钮进行操作。 4.5 动画 制作动画的一般方法为：利用循环把多个图形保存到一个列向量中，然后再播放出来。 for i= 1:k plotting commands M(i) = getframe; %拷贝一个图形窗口到一个列向量中。 end movie(M,n,fps) %播放动画，n 为播放次数，缺省为 1；fps 为帧数/秒，缺省为 12。 例：播放直径不断变化的球 x,y,z=sphere; for i = 1:30 surf(i*x, i*y, i*z); M(i) = getframe; end movie(M,3) % Play the movie 3 times 4.6 声音 sound：将向量 y 传送给扬声器。向量确定了最大振幅。支持大多数 unix 平台。 格式：sound(y,Fs)，Fs 为采样频率，缺省为 8192 Hz wavplay：Play sound on PC-based audio output device 格式：wavplay(y,Fs)，Fs 为采样频率，缺省为 11025 Hz 17 wavrecord：Record sound using PC-based audio input device 格式：y = wavrecord(n,Fs)，Fs 为采样频率，缺省为 11025 Hz wavread：Read Microsoft WAVE (.wav) sound file 格式：y,Fs,bits = wavread(filename,N)，Fs 为采样频率，bits 为采样位数（8 或 16） ，N 为读取的采 样点数。 wavwrite：Write Microsoft WAVE (.wav) sound file 格式：wavwrite(y,Fs,N,filename) 例 1：x=sin(linspace(0, 10000, 10000); % 一个纯正弦波。 sound(x) ; % 产生声音。 例 2：Fs = 11025; y = wavrecord(2*Fs,Fs); wavplay(y,Fs); 习题：教材 p81-42：1、2、3、4 题 18 第 5 章 线性代数 5.1 矩阵 一、 矩阵的行列式 函数：d = det(A)，求得方阵 A 的行列式 d。 若 det(A)=0，则称 A 为奇异矩阵，否则称为非奇异矩阵。 算法：将 A 作 LU 分解。 二、 矩阵的条件数 条件数（2-范数）的定义为：。矩阵 A 的条件数是一个不小于 1 的实数。 1 ( ) |cond AAA 函数 c = cond(A) %求 A 的 2-范数的条件数，即 A 的最大奇异值和最小奇异值的商。 说明 线性方程组 AX=b 的条件数用来衡量关于数据中的扰动，也就是 A 或 b 对解 X 的灵敏度。一 个差条件的方程组的条件数很大。 三、 矩阵的秩 矩阵 A 中最高阶非零子式的阶数称 r 为 A 的秩(rank)。 矩阵的秩代表矩阵中独立方程式个数。如果一矩阵的秩数和其矩阵的列数相等，则此矩阵为满秩矩 阵，否则称为降满秩矩阵。 matlab 函数：k = rank(A) 例：A=2 -1 -1 1 1 1 -2 -1 4 -6 2 -2 3 6 -9 7; rank(A) ans = 3 四、 矩阵的初等变换 下面三种变换称为矩阵的初等行变换： 1） 互换两行； 2） 以一个非零数乘以某一行； 3） 把某一行的 k 倍加到另一行上。 若将定义中的“行”换成“列” ，则称之为初等列变换，初等行变换和初等列变换统称为初等变换。 若矩阵 A 经有限次初等变换变成矩阵 B，则称 A 与 B 等价，记作 AB 。 矩阵 A 作一次初等行变换后得到的矩阵 B 等于以一个相应的初等阵 P 左乘 A，即 B=PA 矩阵 A 作一次初等列变换后得到的矩阵 C 等于以一个相应的初等阵 P 右乘 A，即 C=AP 初等变换不改变矩阵的秩。 对秩为 r 的矩阵 A，必可经有限次初等行变换，使 A 化为阶梯形；经有限次初等变换，使 A 化为标 19 准形，即，其中 Er 为 r 阶单位矩阵。 r E0 00 matlab 函数：R=rref(A) 功能：用初等行变换将矩阵 A 变为行简化阶梯形（Reduced row echolon form）R。其每行第一个非 0 元素均为 1。 例：接上例 R=rref(A) R = 1 0 -1 0 0 1 -1 0 0 0 0 1 0 0 0 0 从而可以看出 A 的秩为 3。 五、 矩阵求逆 若方阵 A、B 满足 AB=BA=E，则称 A、B 是可逆的，B 是 A 的逆矩阵，记为 B=A-1。 一些重要结论： 可逆矩阵 A 的逆矩阵是唯一的； n 阶矩阵 A 可逆的充分必要条件是：det(A)R,jb=rref(A) R = 1 0 -1 0 0 1 -1 0 20 0 0 0 1 0 0 0 0 jb = 1 2 4 jb 表明第 1、2、4 列向量，即 A(:,jb) 是一个极大无关组（主元在 1、2、4 列上） ， 。 二、 向量空间的基与维 若向量空间 V 中的任意一个向量，都可以由一组线性无关的向量 a1, a2, , ar 线性表出，则称向量组 a1, a2, , ar 为 V 的一个基，其向量的个数 r 称为 V 的维数。 向量空间的基就是向量空间这个向量组的极大无关组，维数就是这个向量组的秩。 三、 向量的内积与范数 1. 向量的内积（点积） 若 a=a1, a2, an，b=b1,b2, ,bn，则 a, b 的内积为： =a1b1+a2b2+anbn = ab matlab 函数：C = dot(A,B)。若 A,B 为行向量，相当于 A*B，或 sum (A.*B) 2. 向量的范数（长度） 向量的 2 范数：。若|a|=1，则称 a 为单位向量。 2 1 n i i Vv 2 = = () 对 3 维向量，上式表示空间坐标系下，起点为原点，终点为（v1, v2, v3）的有向线段的长度。 Matlab 计算范数：n=norm(x)，返回 x 的 2 范数 n。 3. 向量的夹角 向量 a, b 的夹角为：cos= / (|a| |b| )，若 =0，则称 a, b 正交。 在空间坐标系中，由于= (|a| |b| ) cos，所以 a, b 的内积为 a 的长度与 b 在 a 上的投影长度值 的乘积，物理上力所作的功就是这一概念。 Matlab 求向量夹角：theta=subspace(A,B)，A, B 为列向量 例：求力 G，N 的合力 F 的大小和方向。 G=0 -100; N=50,50; F=G+N %F=50 -50 Fh=norm(F) %F 的大小，70.7017 theta=subspace(F,G)*180/pi %F 与 G 的夹角，45 度。或 theta=acos(F*G/(norm(G)*norm(F)*180/pi 4. 向量叉乘 两向量的叉乘是一个过两相交向量的交点，且垂直于两向量所在平面的向量。 Matlab 中函数：C = cross(A,B) %若 A、B 为向量，则返回 A 与 B 的叉乘，即 C=AB，A、B 必 须是 3 个元素的向量；若 A、B 为矩阵，则返回一个 3n 矩阵，其中的列是 A 与 B 对应列的叉积， A、B 都是 3n 矩阵。 例：计算垂直于向量(1, 2, 3)和(4, 5, 6)的向量。 a=1 2 3; b=4 5 6; c=cross(a,b) 结果显示：c= -3 6 -3，可得垂直于向量(1, 2, 3)和(4, 5, 6)的向量为(-3, 6, -3)。 21 四、 标准正交基与正交矩阵 若向量 a1, a2, , ar 两两正交，则称为正交向量组。若将 a1, a2, , ar 单位化，则得到标准正交向量 组。 向量空间 V 中由标准正交向量组构成的基称为标准正交基（规范正交基） 。 若方阵满足 AA=AA=E，则称 A 为正交矩阵。 A 为正交矩阵的充要条件是它的 n 个列向量是两两正交的单位向量组。 若 A 是正交矩阵，则：1）A 的行列式等于 1 或-1。 2）A-1=A 正交变换：若 T 为正交矩阵，则线性变换 y=Tx 的正交变换。 正交变换保持向量的长度不变。 Matlab 中函数：Q = orth(A)，求矩阵 A 代表的列向量组的标准正交基。B*B = eye(rank(A) 。 例：接 5.1.2 的例 B=orth(A) B = -0.0584 0.3754 0.1697 -0.1036 0.0946 -0.9806 0.3666 0.8682 -0.0018 -0.9228 0.3105 0.0986 5.3 线性方程组 一、 齐次线性方程组的解 齐次线性方程组 AX=0 有非 0 解的充要条件是：A 的列向量是线性相关的，即 A 的秩 rA=1 2 2 1; 2 1 -2 -2; 1 -1 -4 -3; format rat %指定有理式格式输出 B=null(A,r) %求解空间的有理基 B = 2 5/3 %rank(A)=2，所以 B 有两列 -2 -4/3 1 0 0 1 %从而方程组的通解为： X = 2*k1+5/3*k2 22 -2*k1-4/3*k2 k1 k2 二、 非齐次线性方程组的唯一解或特解 非齐次线性方程组 AX=b 有解的充要条件是：列向量 b 可由 A 的列向量线性表出，即 A 的秩 r1=增 广矩阵A b的秩 r2。方程组有唯一解的充要条件是：b 可由 A 的列向量线性表出，且表示法唯一，即 A 的秩 r1=增广矩阵的秩 r2=A 的列数 n。若 r1=r2A=1 1 -3 ; 3 -1 2; 1 5 -9 ; B=1 4 0; r=rank(A) %求秩 X=AB %求解 X =1.2500 -0.2500 -0.0000 也可用 rref 求解： C=A,B; %构成增广矩阵 R=rref(C) R = 1.0000 0 0 1.2500 0 1.0000 0 -0.2500 0 0 1.0000 0 R 的最后一列元素就是所求之解。 2利用矩阵的分解求方程组的解 1）LU 分解： LU 分解又称 Gauss 消去分解，可把任意方阵分解为下三角矩阵的基本变换形式（行交换）和上三角 23 矩阵的乘积。即 A=LU，L 为下三角阵，U 为上三角阵。则方程 A*X=b 变成 L*U*X=b，所以 X=U(Lb)， 这样可以大大提高运算速度。 格式：L，U=lu (A) 例：解上题 L,U=lu(A) X=U(LB) L = 0.3333 0.2500 1.0000 1.0000 0 0 0.3333 1.0000 0 U = 3.0000 -1.0000 2.0000 0 5.3333 -9.6667 0 0 -1.2500 X =1.2500 -0.2500 -0.0000 2）Cholesky 分解 若 A 为对称正定矩阵，则 Cholesky 分解可将矩阵 A 分解成：A=RR，其中 R 为上三角阵，R为下三 角阵。则方程 A*X=b 变成 RRX=b，所以 X=R(Rb) 。 格式：R = chol(A) 3）QR 分解 对于任何长方矩阵 A，都可以进行 QR 分解，其中 Q 为正交矩阵，R 为上三角矩阵的初等变换形式， 即：A=QR，则方程 AX=b 变成 QRX=b，所以 X=R(Qb) 格式：Q, R=qr(A) 这三种分解在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。 3. 线性方程组的迭代解法 原理：将 Ax=b 改写成一个等价的方程组 x=Bx+g，从而建立迭代式：xk+1=Bxk+g 给出初始向量 x0 后，按公式进行迭代，就可以得到一个向量序列xk，若 xk 收敛于确定的向量 x， 则 x 就是方程组的解。 雅可比迭代法：x k+1=D-1(E+F)x k+D-1b 高斯-赛德尔迭代法：x k+1=(D-E)-1Fx k+(D-E)-1b 线性方程组的 LQ 解法 x = symmlq(A,b) 双共轭梯度法解方程组 x = bicg(A,b) 稳定双共轭梯度方法解方程组 x =bicgstab(A,b) 复共轭梯度平方法解方程组 x = cgs(A,b) 共轭梯度的 LSQR 方法 x = lsqr(A,b) 广义最小残差法 x = gmres(A,b) 最小残差法解方程组 x = minres(A,b) 预处理共轭梯度方法 x = pcg(A,b) 准最小残差法解方程组 x = qmr(A,b) 例 调用 MATLAB 数据文件 west0479。 lo