matlab教案

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 微积分 ............................................................................................................................................. 50 8.3 线性代数 ......................................................................................................................................... 52 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： ?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 注意：

? 语句后加分号,表示仅执行，不显示结果 ? 续行符为…（三个句点） ? 注释：以%开头

2sin(85?)例子：p14，计算表达式的值，并将结果赋给变量a。

1?5?3i>> 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)：在10^a与10^b之间产生n个数，若省略n，则为50 y = logspace(a,pi,n)：在10^a与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阶矩阵，矩阵元素为0~1之

间均匀分布的数；randn(n)、randn(m,n)：0~1之间正态分布。 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矩阵的第i—i+m行元素、第k—k+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

14 15 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、 矩阵的乘法

运算规则：按线性代数中矩阵乘法运算进行，即放在前面的矩阵的各行元素，分别与放在后面的矩阵的各列元素对应相乘并相加。即cij?n?ak?1ik*bkj

例：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、矩阵的除法

有两种除法运算：\\（左除）和/（右除）。A\\B相当于inv(A)*B，B/A相当于B*inv(A)。对于标量来说，/就是通常所用的除法，且A/B和B\\A是等价的。

矩阵的除法在讲线性方程组求解再详细讲解。

10、数组除法：将矩阵的对应元素相除。x.\\y和y./x是等价的。 11、矩阵的乘方：^ 例：A=[1 2; 3 4];

A^2

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，否则为0。 运算优先级：算术运算优先级最高，逻辑运算最低。

2.4 数学函数

2.4.1 常用数学函数

参见p25-26

2.4.2 矩阵的超越函数

一、 矩阵指数

Y = expm(X)，etM???k?0?tM?kk!

二、 矩阵对数

Y = logm(X)

通常有：logm(expm(X)) = X = expm(logm(X))

三、 矩阵平方根

Y = sqrtm(X)

习题：

教材p31-32：2 -(1) (2) (3)、3、4题

10 第3章 Matlab程序设计

3.1 m文件

matlab有两种工作方式，一种是交互式的命令行工作方式，一种是m文件的程序工作方式。

m文件是后缀名为.m的文件，m文件有两类：命令文件和函数文件，命令文件没有输入参数，也不返回输出参数。函数文件可以有输入参数，也可返回输出参数。

函数文件和命令文件的执行如同普通的matlab命令。当输入文件名时(如果有自变量就一起附带参数，就执行文件中的语句。

一、 m文件的建立与编辑

所有的M文件都是普通的ASCII文件，可用任何文本编辑器建立和编辑，但一般用matlab自带的编辑器，在其中可以设置程序断点，进行调试。

m文件的建立：File->New->m-file m文件的建立：File->Open…

二、 命令文件

命令文件又称脚本文件（script file），它是将多条需要运行的命令编辑保存到一个文件中，然后在命令窗口中输入文件名，就可运行文件中包含的命令序列。matlab将首先在当前工作目录下寻找此文件，如果它不在当前目录下，那么在该路径下的所有目录中搜索。搜索路径可以通过File->set path…菜单指定。

例：p34

3.2 数据的输入输出

matlab提供了一些输入输出函数，便于用户和计算机之间进行交互。 1. ipput函数

A=input(提示信息,选项)，选项可省略，为's'时允许用户输入字符串。 如：>> A=input('please input A:') please input A:[3 6;1 5] A =

11 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 {commands} end

12 在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会陷入死循环。

6. 当函数M文件到达M文件终点，或者碰到返回命令return，就结束执行和返回。return命令提供了一种结束一个函数的简单方法，而不必到达文件的终点。

13 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必须为等长度的。该命令与plot

15 类似。

例：绘制三维螺旋线 t = 0:pi/50:10*pi; plot3(sin(t),cos(t),t)

二、 mesh函数

用法：h=mesh(X,Y,Z,C)

Y和Z指定的网线面。length(X)=n，length(Y)=m，n]=size(Z)，功能：生成由X，若X与Y均为向量，而[m，

则空间中的点 (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-范数）的定义为：cond(A)?||A||||A?1||。矩阵A的条件数是一个不小于1的实数。 函数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等价，记作A<->B 。

矩阵A作一次初等行变换后得到的矩阵B等于以一个相应的初等阵P左乘A，即B=PA 矩阵A作一次初等列变换后得到的矩阵C等于以一个相应的初等阵P右乘A，即C=AP 初等变换不改变矩阵的秩。

对秩为r的矩阵A，必可经有限次初等行变换，使A化为阶梯形；经有限次初等变换，使A化为标准

19 ?Er形，即??00?，其中Er为r阶单位矩阵。 ?0?

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。

五、 矩阵求逆

-1

若方阵A、B满足AB=BA=E，则称A、B是可逆的，B是A的逆矩阵，记为B=A。 一些重要结论：

? 可逆矩阵A的逆矩阵是唯一的；

? n阶矩阵A可逆的充分必要条件是：det(A)<>0，或rank(A)=n（满秩）；否则为不可逆矩阵。 ? 若A，B可逆，则A-1、AB、A'均可逆且(A-1)-1=A，(AB)-1=B-1A-1，(A')-1=(A-1)'。 逆矩阵的求法

-1

公式法：A=A*/det(A)。A*称为A的伴随矩阵，它是A的每个元素换成其对应的代数余子式，然后再转置得到。

初等行变换法：任何满秩矩阵都能经过一系列初等行变换化为单位矩阵E，同时用一系列同样的初等

-1-1

行变换作用到E上，最终I就化成了A。即：(P1P2…Pm)(C|E)=(E|A)

matlab函数：B=inv(A)，求方阵A的逆矩阵B。若A为奇异阵或近似奇异阵，将给出警告信息。 例；接上例，B=inv(A)，会提示Warning: Matrix is singular to working precision.

5.2 向量空间

一、 向量的线性相关性

对一组n个向量v1, v2, …, vn，若存在不全为0的数k1, k2, …, kn，使成立： k1v1 + k2v2 + … + knVn = 0， 称这n个向量是线性相关的，否则是线性无关的。 若向量组中r个向量线性无关，而任意r+1个向量线性相关，则称向量组的秩为r。 若向量组的秩为r，则向量组中任意r个线性无关的向量都称为向量组的极大无关组。

矩阵的秩等于其行向量组的秩，也等于其列向量组的秩。因此可用rank()函数求向量组的秩。

以列向量构建矩阵A，把A用初等行变换化为阶梯形矩阵B，则根据B的主元可求出极大无关组。 例：5.1.2的例子中， >>[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范数：V2==(v)。若||a||=1，则称a为单位向量。 ?i2i?1n对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=A×B，A、B必须是3个元素的向量；若A、B为矩阵，则返回一个3×n矩阵，其中的列是A与B对应列的叉积，A、B都是3×n矩阵。

例：计算垂直于向量(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'=A'A=E，则称A为正交矩阵。

A为正交矩阵的充要条件是它的n个列向量是两两正交的单位向量组。

-1

若A是正交矩阵，则：1）A的行列式等于1或-1。 2）A=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的秩r

解集合的极大无关组称为解的一个基础解系，基础解系包含n-r个解向量。 如果AX=0有基础解系X1, X2, …Xn-r，则其通解可表示为： k1X1 + k2X2+ … + k n-r X n-r，这里k1, k2, … k n-r为任意常数。

在Matlab中，函数null用来求出齐次线性方程组AX=0的一组基础解系。

格式 z = null(A) % z的列向量为方程组的正交规范基，满足Z??Z?I。利用svd()求得。

z = null(A,'r') % z的列向量是方程AX=0的有理基。利用rref()求得，便于教学。

?x1?2x2?2x3?x4?0?例 求解方程组的通解：?2x1?x2?2x3?2x4?0

??x1?x2?4x3?3x4?0解：

>>A=[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] [ -2*k1-4/3*k2] [ k1]

22 [ k2]

二、 非齐次线性方程组的唯一解或特解

非齐次线性方程组AX=b有解的充要条件是：列向量b可由A的列向量线性表出，即A的秩r1=增广矩阵[A b]的秩r2。方程组有唯一解的充要条件是：b可由A的列向量线性表出，且表示法唯一，即A的秩r1=增广矩阵的秩r2=A的列数n。若r1=r2

线性方程组的无穷解 = 对应齐次方程组的通解 + 非齐次方程组的一个特解。

线性方程组的唯一解或特解的数值解法有两类：一类是直接法，主要用于解低阶稠密矩阵；另一类是迭代法，用于解大型稀疏矩阵。

一般不用矩阵求逆的方法解线性方程组，因为计算量较大，而且精度较低。（若AX=B，且A可逆，

-1

则方程的解为：X= AB 。） 直接法是在不考虑舍入误差的情况下，通过有限次四则运算求得准确解。克莱姆法则就是一种直接法，即x1=D1/D，x2=D2/D，…，xn=Dn/D，D为系数行列式，Dj为把系数行列式中第j列用方程组右端的常数列代替后得到的行列式。但当方程组阶数较高时，克莱姆法则运算量很大，无法实用。

常用而有效的直接法是消元法，这种解法先把方程组化成同解的上三角形方程组（称为消元过程），然后按方程相反顺序解上三角形方程组（称为回代过程），得出原方程的解。

1．利用矩阵除法求线性方程组的特解 方程：AX=b 解法：X=A\\b

?x1?x2?3x3?1?例 求方程组?3x1?x2?2x3?4的一个特解。

?x?5x?9x?023?1解：

>>A=[1 1 -3 ; 3 -1 2; 1 5 -9 ]; B=[1 4 0]'; r=rank(A) %求秩 X=A\\B %求解 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消去分解，可把任意方阵分解为下三角矩阵的基本变换形式（行交换）和上三角矩阵的乘积。即A=LU，L为下三角阵，U为上三角阵。则方程A*X=b 变成L*U*X=b，所以X=U\\(L\\b)，这样可以大大提高运算速度。

23 格式：[L，U]=lu (A) 例：解上题 >>[L,U]=lu(A) X=U\\(L\\B)

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=R'R，其中R为上三角阵，R'为下三角阵。则方程A*X=b 变成R'RX=b，所以X=R\\(R'\\b) 。

格式：R = chol(A)

3）QR分解

对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形式，即：A=QR，则方程 AX=b 变成QRX=b，所以X=R\\(Q\\b)

格式：[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。 >> load west0479

>> A=west0479; %将数据取为系数矩阵A。

>> b=sum (A,2); %将A的各行求和，构成一列向量。 >> X=A\\b; %用“\\”求AX=b的解。

24 >> norm(b-A*X)/norm(b) %计算解的相对误差。 ans =

1.2454e-017

>> [x,flag,relres,iter,resvec] = bicg(A,b) %用bicg函数求解。 x = （全为0，由于太长，不显示出来） flag =

1 %表示在默认迭代次数（20次）内不收敛。

relres = %相对残差relres = norm(b-A*x)/norm(b) =norm(b)/norm(b) = 1。 1

iter = %表明解法不当，使得初始估计值0向量比后来所有迭代值都好。 0

三、 求非齐次线性方程组的通解

非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解。 因此，步骤为：

第一步：判断AX=b是否有解，若有解则进行第二步 第二步：求AX=b的一个特解 第三步：求AX=0的通解

第四步：AX=b的通解 = AX=0的通解 + AX=b的一个特解。

?x1?2x2?3x3?x4?1?例 求解方程组?3x1?x2?5x3?3x4?2

?2x?x?2x?2x?3234?1解：在Matlab中建立M文件如下：

A=[1 -2 3 -1; 3 -1 5 -3; 2 -1 2 -2]; b=[1 2 3]'; B=[A b]; n=4;

R_A=rank(A) R_B=rank(B) format rat

if R_A==R_B&R_A==n %判断有唯一解 X=A\\b

elseif R_A==R_B&R_A

C=null(A,'r') %求AX=0的基础解系 else X='equition no solve' %判断无解 end

运行后得到：

R_A = 3 R_B = 3 Xt = 2 -1 -1 0

25 C = 1 0 0 1

从而通解为：X=Xt+k*C

四、 超定方程组

对于一个系数矩阵是m×n的线性方程组A x = b来说，如果m>n，也就是说方程的个数多于未知数，则称为超定方程组。通常这些方程组是矛盾的，所以方程组没有精确解。在拟合实验数据的曲线时，常会遇到这个问题。

关键是要找到一个向量x使它对m个方程的总误差最小。有几种方法可以求得，但是最常用的方法是最小二乘法。

最小二乘解可以用除法运算符\\或命令nnls和lscov来解得。

MATLAB在用左除解超定方程组时都用QR因式分解。如果方程组的系数矩阵不满秩，就没有唯一的最小二乘解，将不定值设为零。然后给出一个警告信息。

5.4 特征值与二次型

一、 特征值与特征向量的基本概念

设A为n阶矩阵，若存在复数λ0和非零n维列向量X使得AX=λX （1）

X为A的属于 则称数λ为A的一个特征值(eigenvalue)，（或对应于）特征值λ的特征向量(eigenvector)。

（1）式可表达为：λX—AX=（λE—A）X=0 （2）

由此可见，特征值λ是使齐次线性方程组（λE—A）X=0 有非零解X≠0的λ值。 方程组（2）有非零解的充要条件是：|λE－A|=0（或秩（λE－A））。

|λE－A|=0称为A的特征方程，|λE－A|称为A的特征多项式。A的特征值就是特征方程的根 广义特征值问题：求方程Ax??Bx解的问题。

若A的特征值为λ1、λ2、λn，则有： λ1+λ2+…+λn=trace(A) λ1λ2…λn=det(A)

二、 特征值与特征向量的求法

雅可比法是用来求实对称矩阵全部特征值与特征向量的方法，其基本思想是把对称阵A经一系列正交相似变换，化为一个对角阵，对角阵的对角元就是A的特征值，而由正交变换的矩阵的乘积可求得特征向量。

QR方法是一种求一般矩阵的全部特征值的有效方法。其原理是：利用矩阵A的QR分解，得到A的一系列相似矩阵。设A=A1=Q1R1，令A2=R1Q1=Q2R2，…，A k=R k-1Q k-1=Q kR k。

由于A k-1=Q k-1R k-1，所以有A k=R k-1Q k-1=Q' k-1A k-1Q k-1，即A k相似于A k-1，因此A1，A2，…，Ak

均相似，它们具有相同的特征值。当k—>∞时，Ak的对角线下（不包括对角元）元素收敛于0，根据相似关系，Ak的对角元收敛于A的特征值。从而可以求出A的特征值，根据特征值由（λE—A）X=0可以解出相应的特征向量。

格式 E= eig(A) %求矩阵A的特征值E，以向量形式存放E。

[V,D] = eig(A) %计算A的特征值对角阵D和特征向量V，使AV=VD成立。特征值在D的对角线上，对应的特征向量为V的列。V的列为单位向量。

26 [V,D] = eig(A,B) %计算广义特征值向量阵V和广义特征值阵D，满足AV=BVD。 例：>> A=[1 -2 2; -2 -2 4; 2 4 -2]; [V,D]=eig(A)

V = 0.3333 0.9339 -0.1293 0.6667 -0.3304 -0.6681 -0.6667 0.1365 -0.7327 D = -7.0000 0 0 0 2.0000 0 0 0 2.0000

即：特征值-7对应特征向量(0.3333 0.6667 -0.6667)T

TT

特征值2对应特征向量(0.9339 -0.3304 0.1365)和(-0.1293 -0.6681 -0.7327)

三、 提高特征值的计算精度

函数 balance

格式 [T,B] = balance(A) %求相似变换矩阵T和平衡矩阵B，满足B?T?1AT。 B = balance(A) %求平衡矩阵B

四、 实对称矩阵的相似对角化

-1

相似矩阵：设A、B为n阶方阵，若存在n阶可逆矩阵P，使PAP=B，A相似于B，记作A~B。称P为A到B的相似变换矩阵。

若A、B相似，则： A、 B的特征值相同；

Bm= (P-1AP)m = P-1AmP，因此若A为对角阵，则很容易通过 P-1AmP求得Bm，这在很多工程问题中有重要作用。

实对称矩阵（A'=A）一定可以相似对角化。[V,D]=eig(A)求得的V就是相似变换矩阵，D就是A的相

-1

似对角矩阵。VAV=D

例：

>>A=[3 2 4; 2 0 2; 4 2 3]; [V,D]=eig(A); inv(V)*A*V ans =

-1.0000 0.0000 -0.0000 0.0000 -1.0000 0.0000 0.0000 -0.0000 8.0000

求线性常微分方程组的非0解

dx1??x2?dt?dx2??x3 ?dt??dx3?dt?12x1?4x2?3x3?A=[0 1 0; 0 0 1; -12 4 3]; [V,D]=eig(A)

V=[V(1,:)./V(1,:);V(2,:)./V(1,:);V(3,:)./V(1,:)]

27 运行结果为

V = 0.10483 -0.21822 -0.21822 0.31449 -0.43644 0.43644 0.94346 -0.87287 -0.87287 D = 3 0 0 0 2 0 0 0 -2 V = 1 1 1 3 2 -2 9 4 4

?x1??1??1??1?????3t??2t???2t即解为：?x1??C1?3?e?c2?2?e?C3??2?e，其中C1, C2, C3为任意常数。

?????x1???9???4???4??

五、 二次型及其标准型

n个变数x1、x2…xn的二次齐次多项式f(x1,x2,…xn)称为二次型。一般将二次型写成如下形式：

nf(x1,x2,...xn)=?aijxixj?X'AX

i,j?1?x1??a11?a1n????????，x???? 式中aij=aji, A????a??x??m1?amn??n?A是n阶[实]对称阵，并称这个实对称阵为二次型的矩阵。显然，二次型与实对称矩阵是一一对应的。

若在上式中对一切i≠j均有aij=0，即二次型的矩阵A为对角阵时，称这样的二次型具有标准形。 若存在n阶可逆矩阵P，使C'AC=B，则称方阵A合同于B，记作A?B。

作线性变换x=Py，则f=X'AX=(Py)'Apy=y'(P'AP)y。因此二次型f能否化为标准形等价于实对称矩阵A能否与对角阵合同。

二次型的标准型不是唯一的，但其正系数与负系数的个数是不变的。

对于实二次型f=x'Ax，若对任意x≠0，都有f>0，则称f为正定二次型，并称矩阵A是正定的，记作A>0。若对任意x≠0，都有f<0，则称f为负定二次型，并称矩阵A是负定的，记作A<0。

实二次型f正定的充要条件可以用以下几种方法表述 ? 二次型的标准型的n个系数全为正。 ? 二次型的对称矩阵A的特征值全为正。

? 霍尔维茨定理：对称矩阵A的各阶顺序主子式全大于零，则其为正定；A的奇数阶主子式为负，

偶数阶主子式为正，则其为负定。

例：求一个正交变换X=PY，把二次型f?x1?2x2?2x3?4x1x2?4x1x3?8x2x3化成标准形。 A=[1 -2 2; -2 -2 4; 2 4 -2]; [V,D]=eig(A)

V = 0.33333 0.93391 -0.12925 0.66667 -0.33042 -0.66812 -0.66667 0.13654 -0.73274

222 28 D = -7 0 0 0 2 0 0 0 2

因此通过正交变换X=VY，化为标准型f??7y1?2y2?2y3 验证

X=[13; 21; 37]; Y=inv(V)*X; X'*A*X Y'*D*Y ans =

3597 ans =

3597

222

29 习题

?1234???1. 分别用rank和rref函数求矩阵?2341?的秩。

??3574???5?10??21?????2. 解矩阵方程：??231?X??20??2X

???2?16???35??3. 已知向量组a1=[1,4,2,1], a2=[-2,1,5,1], a3=[-1,2,4,1], a4=[-2,1,-1,1], a5=[2,3,0,1/3], 求该向量组的

秩，并求一个极大无关组。

4. 已知R4上的两个向量u=[4, 1, 2 , 3]'，v=[-2, 3, 1, 4]'，求u, v的长度，以及u在v上的投影向量。

?400???5. 将矩阵A??031?正交规范化。

?013????5x1?6x2?x1?5x2?6x3?x2?5x3?6x46. 求方程组??x3?5x4?6x5?x4?5x5??1?0?0的解。 ?0?1?x1?x2?3x3?x4?1?7. 求方程组?3x1?x2?3x3?4x4?4的通解。

??x1?5x2?9x3?8x4?0??211???8. 求矩阵A??020?的特征值和特征向量。

??413???9. 将矩阵A=[-1 0 2; 0 1 2; 2 2 0]相似对角化。

10. 将二次型f(x1,x2,x3)=x1-x2?x3?2x2x3?6x1x2化为标准型。

222d2xdx11. 解微分方程2?2?3x?0 （提示：令x=x1, dx/dt=x2）

dtdt

30 第6章 数据处理与多项式

6.1 基本统计处理

? 求最大值c=max(x), [c,I]=max(x) 返回各列中的最大数，I为最大数所在的列。 ? 求最小值min(x), [c,I]=min(x) 例：>>v = [3 -2 4 -1 5 0];

max(v)

ans = 5

>>[vmin, kmin] = min(v) vmin = -2 kmin = 2

? 求中位数median(x)：若数据序列为奇数个数时，中位数就是中间那个数；个数为偶数时，为中间

两个数的平均值。

例：>> x1=[9 –2 5 7 12]; median(x); ans=7

>> x2=[9 –2 5 6 7 12]; median(x); ans=6.5

? 升序排序sort(x) ? 求和sum(x)

? 求平均值mean(x)

例：>> data1 = 2*rand(1,500) + 2; mean(data1)

ans =2.9940

? 求积prod(x)

例：>> nfact = prod(1:4) %求24！ nfact =24

? 求累计和cumsum(x) ? 求累计积cumprod(x) 例：>>cumsum(1:5) ans =

[1 3 6 10 15] >>cumprod(1:5) ans =

1 2 6 24 120

? 求标准方差std(x)：s = std(X,flag) ，flag=0（缺省值）时采用(1)式，flag=1时采用(2)式

31

>> x=randn(1,500); %标准正态分布 std(x)

ans = 0.9764

6.2 多项式

一、 多项式的表示

多项式用向量形式表达，例如对多项式p(x)=x3-2x-5，用p = [1 0 -2 -5]表示。

p=poly(x)，若已知多项式的全部根x1, x2…xn，即(x- x1) (x- x2)…(x- xn)=0，则可用poly函数建立多项式：

x为根组成的向量，p为多项式系数组成的向量。如:

>>p=poly([0 0 0])

p = 1 0 0 0 %p(x)=x3

二、 多项式求根

x=roots(c)

3

例：p = [1 0 -2 -5]; %p(x)=x-2x-5 x=roots(p)

ans = 2.0946 -1.0473 + 1.1359i -1.0473 - 1.1359i

三、 求多项式的值

Y=polyval(A,x)

例：>>p = [1 0 -2 -5]; polyval(p,[5 7 9])

ans = 110 324 706

即x=5, 7, 9时，多项式的值分别为110, 324, 706。

四、 多项式的四则运算

多项式的加、减：直接对系数向量进行加、减（如果多项式的次数不同，需要把低次的多项式系数不足的高次项用0补足）。

多项式乘法：w = conv(u,v)

多项式除法：[q,r] = deconv(u,v)， q为商，r为余项，即u= conv(v,q)+r 例：

>>u=[1 8 0 0 -10]; v=[2 -1 3];

w = conv(u,v) %求x4+8x3-10与2x3-x3+3的乘积

w = 2 15 -5 24 -20 10 -30 %乘积为 [q,r] = deconv(v,u)

q = 0.5000 4.2500 1.3750

32 r = 0 0 0 -11.3750 -14.1250

6.3 数据插值

插值法是实用的数值方法，是函数逼近的重要方法。在生产和科学实验中，自变量x与因变量y的函数y = f(x)的关系式有时不能直接写出表达式，而只能得到函数在若干个点的函数值或导数值。当要求知道观测点之外的函数值时，需要估计函数值在该点的值。

如何根据观测点的值，构造一个比较简单的函数y=φ(x)，使函数在观测点的值等于已知的数值或导数值。用简单函数y=φ(x)在点x处的值来估计未知函数y=f(x)在x点的值。寻找这样的函数φ(x)，办法是很多的。φ(x)可以是一个代数多项式，或是三角多项式，也可以是有理分式；φ(x)可以是任意光滑（任意阶导数连续）的函数或是分段函数。函数类的不同，自然地有不同的逼近效果。

如拉格朗日n次插值多项式可写为：

Ln(x)??yili(x)??i?1i?1nn(x?x0)...(x?xi?1)(x?xi?1)...(x?xn)yi

(xi?x0)...(xi?xi?1)(xi?xi?1)...(xi?xn)特别地，一次插值（线性插值）多项式为：

L1(x)?x?x0y?y0x?x1y0?y1 或 L1(x)?y0?1(x?x0)

x0?x1x1?x0x1?x0

命令1 interp1

功能 一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。各个参量之间的关系示意图为图2-14。

Y：原始数据点 Yi：插值点 f(x) x：原始数据点 xi：插值点

图 6-1 数据点与插值点关系示意图

格式 yi = interp1(x,Y,xi) %返回插值向量yi，每一元素对应于参量xi，同时由向量x与Y的

yi是阶数为length(xi)*size(Y,2)内插值决定。参量x指定数据Y的点。若Y为一矩阵，则按Y的每列计算。

的输出矩阵。

yi = interp1(Y,xi) %假定x=1:N，其中N为向量Y的长度，或者为矩阵Y的行数。

yi = interp1(x,Y,xi,method) %用指定的算法计算插值：’nearest’：最近邻点插值；’linear’：线性插值

（缺省方式）；’spline’：三次样条函数插值。’pchip’：分段三次Hermite插值。’cubic’：与’pchip’操作相同；

yi = interp1(x,Y,xi,method,'extrap') %对于超出x范围的xi中的分量将执行特殊的外插值法extrap。 yi = interp1(x,Y,xi,method,extrapval) %确定超出x范围的xi中的分量的外插值extrapval，其值通常取NaN或0。

例：已知1900年至1990年中每隔10年的美国人口数（万） t = 1900:10:1990;

33 p = [7599 9197 10571 12320 13166 15069 17932 20321 22650 24963]; interp1(t,p,1975) %求1975年的人口数 ans = 21486

x = 1900:1:2000; %绘制分段线性插值曲线 y = interp1(t,p,x,'spline'); plot(t,p,'o',x,y)

图 6-2 一维数据插值

命令2 interp2

功能 二维数据内插值

格式 ZI = interp2(X,Y,Z,XI,YI) %返回矩阵ZI，其元素包含对应于参量XI与YI（可以是向量、或同型矩阵）的元素，即Zi(i,j)←[Xi(i,j),yi(i,j)]。用户可以输入行向量Xi和列向量Yi，此时输出向量Zi与矩阵meshgrid(xi,yi)是同型的。同时取决于由输入矩阵X、Y与Z确定的二维函数Z=f(X,Y)。参量X与Y必须是单调的，且相同的划分格式，就像由命令meshgrid生成的一样。若Xi与Yi中有在X与Y范围之外的点，则相应地返回nan（Not a Number）。

ZI = interp2(Z,XI,YI) %缺省地，X=1:n、Y=1:m，其中[m,n]=size(Z)。再按第一种情形进行计算。 ZI = interp2(Z,n) %作n次递归计算，在Z的每两个元素之间插入它们的二维插值，这样，Z的阶数将不断增加。interp2(Z)等价于interp2(z,1)。

ZI = interp2(X,Y,Z,XI,YI,method) %用指定的算法method计算二维插值： ’linear’：双线性插值算法（缺省算法）；’nearest’：最临近插值；’spline’：三次样条插值；’cubic’：双三次插值。

例：平板温度

length=1:5; % index for length of plate (i.e.,the x-dimension) with=1:3; % index for with of plate (i,e,,the y-dimension)

temps=[82 81 80 82 84; 79 63 61 65 81; 84 84 82 85 86]; %求长2.5、宽1.2处的温度

interp2(length, with, temps, 2.5, 1.2) ans = 76.8

%绘制分段线性插值曲面

[l,w]=meshgrid(1:0.2:5,1:0.2:3); %或l=[1:0.2:5]; w=[1:0.2:3]' choose higher resolution t=interp2(length, with, temps, l, w, 'cubic') ; % linear interpolation mesh(l, w, t)

xlabel(' length of Plate '), ylabel(' with of Plate ') zlabel(' Degrees Celsius ')

34

图 6-3 二维数据内插值

命令3 interp3

功能 三维数据插值（查表）

格式 VI = interp3(X,Y,Z,V,XI,YI,ZI) %找出由参量X,Y,Z决定的三元函数V=V(X,Y,Z)在点（XI,YI,ZI）的值。参量XI,YI,ZI是同型阵列或向量。若向量参量XI,YI,ZI是不同长度，不同方向（行或列）的向量，这时输出参量VI与Y1,Y2,Y3为同型矩阵。其中Y1,Y2,Y3为用命令meshgrid(XI,YI,ZI)生成的同型阵列。若插值点(XI,YI,ZI)中有位于点(X,Y,Z)之外的点，则相应地返回特殊变量值NaN。

VI = interp3(V,XI,YI,ZI) %缺省地，X=1:N，Y=1:M，Z=1:P，其中，[M,N,P]=size(V)，再按上面的情形计算。

VI = interp3(V,n) %作n次递归计算，在V的每两个元素之间插入它们的三维插值。这样，V的阶数将不断增加。interp3(V)等价于interp3(V,1)。

VI = interp3(?,method) %用指定的算法method作插值计算：‘linear’：线性插值（缺省算法）；‘cubic’：三次插值；‘spline’：三次样条插值；‘nearest’：最邻近插值。

说明 在所有的算法中，都要求X,Y,Z是单调且有相同的格点形式。当X,Y,Z是等距且单调时，用算法’*linear’，’*cubic’，’*nearest’，可得到快速插值。

>>[x,y,z,v] = flow(20);

[xx,yy,zz] = meshgrid(.1:.25:10, -3:.25:3, -3:.25:3); vv = interp3(x,y,z,v,xx,yy,zz);

slice(xx,yy,zz,vv,[6 9.5],[1 2],[-2 .2]); shading interp; colormap cool 插值图形如图。

图 6-4 三维插值图

命令4 spline

功能 三次样条数据插值

格式 yy = spline(x,y,xx) %该命令用三次样条插值计算出由向量x与y确定的一元函数y=f(x)在点xx处的值。若参量y是一矩阵，则以y的每一列和x配对，再分别计算由它们确定的函数在点xx处的值。

35 则yy是一阶数为length(xx)*size(y,2)的矩阵。

pp = spline(x,y) %返回由向量x与y确定的分段样条多项式的系数矩阵pp，它可用于命令ppval、unmkpp的计算。

例：对离散地分布在y=exp(x)sin(x)函数曲线上的数据点进行样条插值计算： >>x = [0 2 4 5 8 12 12.8 17.2 19.9 20]; y = exp(x).*sin(x); xx = 0:.25:20;

yy = spline(x,y,xx); plot(x,y,'o',xx,yy) 插值图形结果如图。

图 6-5 三次样条插值

6.4 曲线拟合

从几何意义上将，拟合是给定了空间中的一些点，找到一个已知形式未知参数的连续曲面来最大限度地逼近这些点；而插值是找到一个(或几个分片光滑的)连续曲面来穿过这些点。

通常采用最小二乘拟合，以最小二乘多项式拟合为例说明。

设有变量x和y的一组数据(xi,yi)，i=1,2…m。最小二乘多项式拟合就是要求取多项式P(x)=a0+a1x+…+anx的系数a0,a1…an (n

[P(xi)?yi]最小。

mi?1m?S?0,k?0,1,...n。即使S取最小值的条件是：

?aksj??X,tj??Xilyi,l?0,1...n

lii?1i?1mm?sj?0nj?ka?jt,kk?0,1...n。 式中

根据这个线性方程组，就可求出P(x)的系数，从而得到最小二乘拟合多项式P(x)。

MATLAB的polyfit函数提供了从一阶到高阶多项式的回归法，其语法为 coef=polyfit(x,y,n)，n=1就是一阶 的线性回归法。其中x,y为输入数据组n为多项式的阶数，从polyfit

coef(1)=a0 , coef(2)=a1,...,coef(n+1)=an 。函数得到的输出值就是多项式的各项系数，注意上式对n 阶的多 项

式会有 n+1 项的系数。

[p,s] = polyfit(x,y,n)，s为一个包含误差估计的结构，其中就有上述的S值。 由 polyfit 算出多项式的各个系数后，即可以 polyval 计算多项式的函数值。 例：t = 1900:10:1990;

p = [7599 9197 10571 12320 13166 15069 17932 20321 22650 24963]; pt=polyfit(t,p,2);

36 ti=1900:5:1990;

pi=polyval(pt,ti); % pt=0.9272 -3414.1 3.1474e+006 注意format为short g plot(t,p,'o',ti,pi)

图 6-6 曲线拟合

此外，通过Figure窗口的tools->basic fitting菜单，打开basic fitting对话框，在其中可以进行曲线拟合。

图 6-7 basic fitting对话框

6.5 离散傅立叶变换

傅立叶变换用于频域分析，即以频率为横坐标，研究信号的幅度、相位在各频率上的分布情况。 傅立叶变换的定义为：F(?)?????f(t)e?j?tdt

F(ω)为复数，可表示为：F(?)?|F(?)|ej?(?)

幅度频谱：|F(ω)|与ω的关系图。相位频谱：φ(ω) 与ω的关系图。

为利用计算机计算傅立叶变换，需要对无限长连续信号，进行时域采样、时域截断、频域采样，周期延拓，从而得到离散傅立叶变换DFT：

X(k)?DFT[x(n)]??x(n)en?0N?1?j2?nkN,0?k?N?1

忽略由频谱混叠和泄漏引起的误差，离散傅立叶变换与傅立叶变换具有下列关系：

37 1DFT[x(t)|t?NTs]?X(f)|f?kf0，0?f?(N?1)f0Ts式中，f0=1/NTs，N为采样点数，Ts为采样周期。

DFT方法计算量太大，限制了应用。1965年美国学者提出了一种快速计算DFT的算法, 称为快速傅立叶变换FFT。matlab中的函数为：

Y = fft(X,n)，n表示只计算x的前n点，省略时计算全部点。

例：信号x(t)=sin(2πf1) +2sin(2πf2)+N。式中，f1 =50，f2=120，N为正态分布的噪声信号。 t=0:0.001:0.6; `0点

x=sin(2*pi*50*t)+2*sin(2*pi*120*t)+2*randn(size(t)); subplot(2,1,1)

plot(x(1:256)) %绘制信号的时域波形 y=fft(x,512); %快速傅立叶变换

f=(0:256)/(0.001*512); %频率轴，只绘制前面257点 subplot(2,1,2)

plot(f,abs(y(1:257))) %绘制信号的幅度频谱 信号的时域波形和幅度频谱如下。从时域波形只能看出信号的振幅范围，很难看出信号的频率特征；而从频谱图可以很容易看出信号中主要由频率分别为50与120的简谐信号构成，且后者的幅度是前者的两倍。

图 6-8 (a)信号的时域波形 (b) 信号的幅度频谱

38 第7章 微积分

7.1 数值微分

函数f(x)在x=a的微分可用差分近似表示。 前向差分：?fk??fk?1??fk 后向差分：?fk??fk??fk?1 中间差分：?fk??fk?h??fk?h

22在MATLAB用diff 函数来计算二相邻点的差值，语法为d=diff(y)。其中 y为一系列的数据向量，d为差分值输出。d=diff(y,n)计算n阶差分。

函数y=f(x)的导数则为 dy=diff(y)./diff(x)。 >> x=[1 4 9 16 25]; >> dx=diff(x)

dx = 3 5 7 9 >> dx2=diff(x,2)

dx2 = 2 2 2

现假设有一组流率与压降之关系的实验数据，如下表: x(流率,L/s) 0 10.80 16.03 22.91 32.56 36.76 39.88 43.68 y(压降,KPa) 0 0.299 0.576 1.036 1.781 2.432 2.846 3.304

现利用diff来进行後向差分,藉以求得压降对流率之变化率,并与中间差分之微分近似值结果做一比较. >> x=[0 10.80 16.03 22.91 32.56 36.76 39.88 43.68]; >> y=[0 0.299 0.576 1.036 1.781 2.432 2.846 3.304]; >> n=length(x);

>> d1=diff(y)./diff(x); %前、后向差分

>> d2=0.5*(d1(1:n-2)+d1(2:n-1)); %中间差分

>> plot(x(2:n),d1,'o--',x(1:n-1),d1,'*-.',x(2:n-1),d2,'s-') >> xlabel('x')

>> ylabel('微分近似值')

>> legend('后向差分','前向差分','中间差分')

7.2 数值积分

bn函数f(x)在区间[a,b]上的定积分为： S=S??af(x)dx?lim?f(xi)n??i?1b?a 2定积分的数值解法的基本思想是：将积分区间[a,b]分成n个子区间，再在每个子区间上求定积分的近

39 似值。即S???i?1nxi?1xif(x)dx。

图 7-1 数值积分

梯形求积：在每个子区间上用线性插值代替f(x)，从而得到S的数值解为：

n?1hn?111T=?[f(xk)?f(xk?1)]?h[f(a)?f(b)??Xk]，h=(b-a)/n，称为步长。

2k?022k?1h2梯形求积公式的误差为：R?(b?a)f''(?)，??[a,b]

12辛普生求积(Simpleson quadrature)：在每两个相邻的子区间上用二次插值代替f(x)，从而得到S的数值解为：

nn?1hh4T=[f(a)?f(b)?4?f(X2k?1)?2?f(X2k)]，其误差为：R?(b?a)f(4)(?)

3180k?1k?1实际应用时，往往采用逐步缩小步长的办法，直到两次的计算结果相差不大时为止，称为变步长积分

法。缩小步长最好每次减半，这样前次计算过的函数值不必重复计算。

一、 一元函数的数值积分

函数1 quad、quadl、quad8

功能 求函数的数值积分。quad采用变步长Simpleson积分法， quad8采用高阶方法，新版本中quad8用quadl取代，quadl采用变步长Gauss/Lobatto方法。

格式 q = quad(fun,a,b) %近似地从a到b计算函数fun的数值积分，误差为10-6。若给fun输入向量x，应返回向量y，即fun是一单值函数。

q = quad(fun,a,b,tol) %用指定的绝对误差tol代替缺省误差。tol越大，函数计算的次数越少，速度越快，但结果精度变小。

q = quad(fun,a,b,tol,trace,p1,p2,?) %将可选参数p1,p2,…等传递给函数fun(x,p1,p2,…)，再作数值积分。若tol=[]或trace=[]，则用缺省值进行计算。

[q,n] = quad(fun,a,b,?) %同时返回函数计算的次数n

… = quadl(fun,a,b,?) %用高精度进行计算，效率可能比quad更好。

例：求

??0sin(x)dx

quad('sin(x)',0,pi) ans = 2.0000

函数2 trapz

功能 梯形法数值积分

40 格式 T = trapz(Y) %用等距梯形法近似计算Y的积分。若Y是一向量，则trapz(Y)为Y的积分；若Y是一矩阵，则trapz(Y)为Y的每一列的积分；若Y是一多维阵列，则trapz(Y)沿着Y的第一个非单元集的方向进行计算。

T = trapz(X,Y) %用梯形法计算Y在X点上的积分。若X为一列向量，Y为矩阵，且size(Y,1) = length(X)，则trapz(X,Y)通过Y的第一个非单元集方向进行计算。

例： 求

??0sin(x)dx

>>X = 0:pi/100:pi; Y = sin(X); Z = trapz(X,Y) Z = 1.9998

二、 重积分的数值计算

函数1 dblquad

功能 矩形区域上的二重积分的数值计算

格式 q = dblquad(fun,xmin,xmax,ymin,ymax) %调用函数quad在区域[xmin,xmax, ymin,ymax]上计算二元函数z=f(x,y)的二重积分。输入向量x，标量y，则f(x,y)必须返回一用于积分的向量。

q = dblquad(fun,xmin,xmax,ymin,ymax,tol) %用指定的精度tol代替缺省精度10-6，再进行计算。 q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method) %用指定的算法method代替缺省算法quad。method的取值有@quadl或用户指定的、与命令quad与quadl有相同调用次序的函数句柄。

q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method,p1,p2,?) %将可选参数p1,p2,..等传递给函数fun(x,y,p1,p2,…)。若tol=[]，method=[]，则使用缺省精度和算法quad。

例

>>fun = inline(’y./sin(x)+x.*exp(y)’); >>Q = dblquad(fun,1,3,5,7) Q = 3.8319e+003

函数2 quad2dggen（在NIT工具箱中） 功能 任意区域上二元函数的数值积分。

格式 q = quad2dggen(fun,xlower,xupper,ymin,ymax) %在由[xlower,xupper, ymin,ymax]指定的区域上计算二元函数z=f(x,y)的二重积分。

例 计算单位圆域上的积分：I?x2?y2?111?y2?x22??e?x22sin(x2?y)dxdy

先把二重积分转化为二次积分的形式：I?dy?1??1?y2?esin(x2?y)dx

f = inline(’exp(-x.^2/2).*sin(x.^2+y)’,’x’,’y’);

xlower = inline(’-sqrt(1-y.^2)’,’y’); xupper = inline(’sqrt(1-y.^2)’,’y’); Q = quad2dggen(fun,xlower,xupper,-1,1,1e-4) 计算结果为：Q = 0.5368603818

函数3：triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol)：长方体区域的三重积分

函数4：quadndg('Fun',xlow,xhigh,tol)，在NIT工具箱中，n-D hyper-rectangular 区域的n重积分。

41

7.3 常微分方程的解

微分方程是研究函数变化规律的重要工具，有着非常广泛的应用。

若微分方程的自变量个数为1则称为常微分方程（ODE，Ordinary Differential Equation）；若自变量的个数大于1，则称为偏微分方程。一般的，n阶常微分方程具有形式：

dxdnxF(t,x,,...,n)?0

dtdtdxdnx若式子左端为x及,...,n的一次有理整式，则称为n阶线性微分方程。一般的，n阶线性常微分

dtdt方程具有形式：

dnxdn?1xdx?a(t)?...?a(t)?an(t)x?f(t) 1n?1nn?1dtdtdt高阶线性微分方程一般可化为一阶线性微分方程组，方法为：令xn=dxn-1/dt, xn-1=dxn-2/dt,?,x2=dx/dt, x1=x，则高阶的方程（组）可写成一阶微分方程组：X'=AX+f(t)。

若上式中f(t)=0，则称为齐线性方程组。n阶齐线性方程组X'=AX一定存在n个线性无关的解：x1(t), x2(t), ?xn(t)，x=C1x1(t)+ C2x2(t)+ ?+Cnxn(t)，x1(t), x2(t), ?xn(t)且其通解可表为：其中Ci为任意常数。称为方程组的基础解系，其构成的矩阵称为基解矩阵。

非齐线性方程组的通解 = 对应齐线性方程组的通解 + 非齐线性方程组的一个解。 若式中A为常数，则称为常系数线性微分方程组。 为了确定微分方程的特定解，需要给出这个解必须满足的条件，即定解条件。如果定解条件是初始条件，就称为初值问题（IVP, Initial Value Problem）。

一、 常系数齐线性微分方程组的代数解法

At

定理：矩阵e是常系数齐线性微分方程组X'=AX的基解矩阵。

?ea1t?k?(At)AteAt??，如果A是一个对角矩阵，则e???k!k?1??ea2t???，式中a1,a2?为A的对角?...?eant?At

线元素。若A是实的，则e也是实的

定理：如果A具有n个线性无关的特征向量v1, v2,?vn，它们对应的特征值分别为λ1,λ2,?λn（不

必各不相同），则矩阵φ(t)=[e1v1,e2vn,...,evn]是齐线性微分方程组X'=AX的一个基解矩阵。 关于这一结论，在第5章举例说明过。

At-1

如果A是实的，则通过e=Φ(t)Φ(0)，可以得到一个实的基解矩阵。

-1

如果A的特征向量并非是线性无关的，则不能通过相似变换VAV化成对角矩阵，而只能化成准对角矩阵或约当标准型（Jordan Canonical Form）矩阵，利用约当标准型矩阵可计算基解矩阵。

λtλtλnt

42 ?J1?约当标准型矩阵J?????kJ2?λj?????，其中J????...??Jk?n?n??1λj??1??......?1?λj??称为约当块，为nj阶矩

nj?nj阵，且

?nj?1j?n。λ1，λ2，…λk是A的特征值，nj是λj的重特征值数。

-1从而基解矩阵eAt=e(VJV)t?VeJtV?1，而

?eJ1t?eJt=????eJ2t??1t????Jtj?，其中e=?01??...??eJkt????t2...2!tn-1tj??nj?1?n-2?tj?λjte ...nj?2??...?1??matlab符号工具箱中的函数jordan(A)可以求出非奇异的矩阵V，使得inv(V)*A*V成为约当标准型，

函数jordan的形式为[V,J]=jordan(A)，V为相似变换矩阵，J为约当标准型矩阵。

例：>>A=[2 1; -1 4];

[V,J]=jordan(A) % find the Jordan form and eigenvectors V =

-1 1 -1 0 J =

3 1 0 3 iV=inv(V) iV =

0 -1 1 -1

t???11??1t?3t?0?1?3t?1?tAte=e?e从而基解矩??10??01??1?1???t1?t?

????????

也可利用符号工具箱求得： >> syms t

V*[1 t; 0 1]*exp(3*t)*inv(V) ans =

[ exp(3*t)*(-t+1), exp(3*t)-exp(3*t)*(-t+1)] [ -exp(3*t)*t, exp(3*t)+exp(3*t)*t]

43 二、 利用拉普拉斯变换解齐线性微分方程组

（略）

三、 常微分方程的数值解

一个一阶常微分方程式 (ordinary differential equation, ODE) 可以下式表示

?y'?f(t,y),t0?t?T ?y(t)?y0?0其中t为独立变数而 y是t的函数。数值解法就是要求方程的解在点t0

以数值方法求解上述的 ODE 的问题，可以依据泰勒级数对 y(x) 做展开。一个一阶的泰勒级数近似式为 y(x k+1)=y(x k+h)=y(x k)+hy'(x k)=y(x k)+hf(x k,y(x k))

这就是欧拉公式，根据这个公式就可以求得个x k处的解y(x k)。

图 7-2 欧拉公式

龙格-库塔 (Runge-Kutta) 方法是最通用的解 ODE 的方法，它可以依计算精确度的要求有低阶到高阶的各个计算式。

根据微分中值定理，有y'(xk??h)?y(xk?1)?y(xk)yk?1?yk? 式中0<θ<1

xk?1?xkh而y'?f(x,y)，y'(xk??h)?f(xk??h,y(xk??h))，故有yk?1?yk?hf(xk??h,y(xk??h)) 由于f (x k+θh, y(x k+θh))是(x k,x k+1)中满足微分中值定理的点的导数，因此f(x k+θh, y(x k+θh))称

为平均斜率，记作K，则y k+1=y k+hK。

如取x k处的斜率f (x k,y k)作为平均斜率k的近似值，则得欧拉公式。

如果通过在[x k, x k+1]中多取几个点，再以它们的某种平均值作为平均斜率，就可以提高精度。

?，xk??nh，再以???1?1??2?2????n?n，假设在[x k, x k+1]中取n个点：xk??1h，xk??2h，则：

yk?1?yk?h??m?m 式中?m?f(xk??mh,yk?h??mj?j)

m?1j?1nm?1这就是龙格—库塔法的一般公式。恰当的选取公式中的常数λ、α、β，就可以使精度尽量地高。λ、

α、β可以用待定系数法确定，确定的原则是将局部截断误差按泰勒级数展开，选取系数使它关于h的阶数尽可能高一些。

函数 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb 功能 常微分方程（ODE）组初值问题的数值解

格式 [T,Y] = solver(odefun,tspan,y0) %在区间tspan=[t0,tf]上，从t0到tf，用初始条件y0求解显式

44 微分方程y’=f(t,y)。对于标量t与列向量y，函数f=odefun(t,y)必须返回一f(t,y)的列向量f。解矩阵Y中的每一行对应于返回的时间列向量T中的一个时间点。要获得问题在其他指定时间点t0,t1,t2,…上的解，则令tspan=[t0,t1,t2,…,tf]（要求是单调的）。

[T,Y] = solver(odefun,tspan,y0,options) %用参数options（用命令odeset生成）设置的属性（代替了缺省的积分参数），再进行操作。常用的属性包括相对误差值RelTol（缺省值为1e-3）与绝对误差向量AbsTol（缺省值为每一元素为1e-6）。

[T,Y] =solver(odefun,tspan,y0,options,p1,p2…) 将参数p1,p2,p3,..等传递给函数odefun，再进行计算。若没有参数设置，则令options=[]。

参数说明：

solver为命令ode45、ode23,ode113,ode15s,ode23s,ode23t,ode23tb之一。因为没有一种算法可以有效地解决所有的ODE问题，为此，MATLAB提供了多种求解器Solver，对于不同的ODE问题，采用不同的Solver。 ode45 ode23 ode113 ode23t ode15s ode23s ode23tb 非刚性 非刚性 非刚性 适度刚性 刚性 刚性 刚性 一步算法；4，5阶Runge-Kutta方程；累计截断误差达(△x)3 一步算法；2，3阶Runge-Kutta方程；累计截断误差达(△x)3 多步法；Adams算法；高低精度均可到10-3～10-6 采用梯形算法 多步法；Gear’s反向数值微分；精度中等 一步法；2阶Rosebrock算法；低精度 梯形算法；低精度 大部分场合的首选算法 使用于精度较低的情形 计算时间比ode45短 适度刚性情形 若ode45失效时，可尝试使用 当精度较低时，计算时间比ode15s短 当精度较低时，计算时间比ode15s短 Odefun 为显式常微分方程y’=f(t,y)，或为包含一混合矩阵的方程M(t,y)*y’=f(t,y)。命令ode23只能求解常数混合矩阵的问题；命令ode23t与ode15s可以求解奇异矩阵的问题。

Tspan 积分区间（即求解区间）的向量tspan=[t0,tf]。要获得问题在其他指定时间点t0,t1,t2,…上的解，则令tspan=[t0,t1,t2,…,tf]（要求是单调的）。

Y0 包含初始条件的向量。

Options 用命令odeset设置的可选积分参数。 P1,p2,… 传递给函数odefun的可选参数。

1．求解具体ODE的基本过程：

（1）根据问题所属学科中的规律、定律、公式，用微分方程与初始条件进行描述。

F(y,y’,y’’,…,y(n),t) = 0

y(0)=y0,y’(0)=y1,…,y(n-1)(0)=yn-1

而y=[y;y(1);y(2);…,y(m-1)]，n与m可以不等

（2）运用数学中的变量替换：yn=y(n-1),yn-1=y(n-2),…,y2=y1=y，把高阶（大于2阶）的方程（组）

?y???y1(0)??y0??f1(t,y)??1???f(t,y)??y(0)??y???2y?????1? 2y0??2写成一阶微分方程组：y???，???????????????????????fn(t,y)??yn(0)??yn??yn??（3）根据（1）与（2）的结果，编写能计算导数的M-函数文件odefile。

（4）将文件odefile与初始条件传递给求解器Solver中的一个，运行后就可得到ODE的、在指定时间区间上的解列向量y（其中包含y及不同阶的导数）。

d2ydy例 求解描述振荡器的经典的Ver der Pol微分方程2??(1?y2)?1?0

dtdty(0)=1，y’(0)=0

令x1=y，x2=dy/dx，则

45 dx1/dt = x2

dx2/dt = μ(1-x2)-x1

编写函数文件verderpol.m： function xprime = verderpol(t,x) global MU

xprime = [x(2);MU*(1-x(1)^2)*x(2)-x(1)]; 再在命令窗口中执行： >>global MU >>MU = 7; >>Y0=[1;0]

>>[t,x] = ode45(‘verderpol’,0,40,Y0); >>x1=x(:,1);x2=x(:,2); >>plot(t,x1,t,x2)

图 7-3 Ver der Pol微分方程图

7.4 非线性方程(组)求解

一、 非线性方程的解

非线性方程的标准形式为f(x)=0

求方程实根的最简单的近似方法是二分法，但其求解时间长，一般不单独采用。求方程根的主要方法是迭代法。迭代法的原理如下：

将方程f(x)=0化为一个等价（同根）的方程：x=ψ(x)。

给定一个初值x0，代入右端可算得一个x1=ψ(x0)，再将x1代入右端，又可算得x2=ψ(x1)，如此继续下去，得到一个序列{xk}，其中x k+1=ψ(x k)，{xk}称为迭代序列，ψ(x)称为迭代函数。如果迭代序列是收敛的，且收敛于x*，则x*就是方程的根。实际计算时，当两次迭代的值相差很小时，就可取x k+1作为方程的根。

如何构造迭代函数，才能保证迭代序列收敛呢？可以将f(x)在xk处作泰勒展开，并取近似值，得 f(x)=f(x k)+f'(x k)(x-x k)，从而得到近似的方程f(x k)+f'(x k)(x-x k)=0。

所以xk?1?xk?f(xk)f(xk)?xk?(xk?xk?1)

f'(xk)f(xk)?f(xk?1)

函数 fzero

说明 该函数采用数值解求方程f(x)=0的根。

46 格式 x = fzero (fun,x0) %用fun定义表达式f(x)，x0为初始解。

[x,fval] = fzero(fun,x0) %fval=f(x) 例：求x3-2x-5=0的根 >> fun='x^3-2*x-5';

z=fzero(fun,2) %初始估计值为2 z =

2.0946

[x,fval] =fzero(fun,2) x =

2.0946 fval =

-8.8818e-016

二、 非线性方程组的解

非线性方程组的标准形式为：F(X) = 0，其中X为向量，F(X)为函数向量。 求解非线性方程组常用办法是迭代法，把方程组化为如下的等价形式：

li (x1,x2,…xn)=ψ(x1,x2,,…,xn)。

其中左端是线性式，建立迭代式：

(k?1)?1)?1)(k)(k)li(x1,x(k,...,x(k)??i(x1,x(k)2n2,...,xn),(i?1,2,...,n),k?0,1,2...

取初值(x1,x2,...,xn)后，按上式迭代，右端成为已知量，因而每步迭代需要求解一个线代数方程组，解这个线代数方程组若用迭代法，就形成一个双重迭代。

函数 fsolve

格式 x = fsolve(fun,x0) %用fun定义向量函数，定义方式为：

先定义方程函数function F = myfun (x)。F =[表达式1；表达式2；?表达式m] %保存为myfun.m，并用下面方式调用：x = fsolve(@myfun,x0)，x0为初始估计值。

x = fsolve(fun,x0,options)

[x,fval] = fsolve(?) %fval=F(x)，即函数值向量 例: 求下列系统的根

?x??2x1?x2?e1 ??x2???x1?2x2?e?x??2x1?x2?e1?0解：化为标准形式?，设初值点为x0=[-5 -5]。 ?x2???x1?2x2?e?0(0)(0(0)先建立方程函数文件，并保存为myfun.m：

function F = myfun(x)

F = [2*x(1) - x(2) - exp(-x(1)); -x(1) + 2*x(2) - exp(-x(2))]; 然后调用优化程序

x0 = [-5; -5]; % 初始点 [x,fval] = fsolve(@myfun,x0) x =

0.5671

47 0.5671 fval =

1.0e-006 * -0.4059 -0.4059

7.5 函数优化

最优化是求最优解，也就是在某个区间内有条件约束或者无条件约束地找到函数的最大值或者最小值。MATLAB使用数字方法求函数的最小值。使用迭代算法，也就是有些步骤要重复许多次。现在，假设要求函数f在某个区间内的最小值xmin。

迭代方法需要一个初始估计值x0。从x0开始找到一个更接近xmin的新值x1，这个值的好坏取决于使用的数学方法。直到找到有足够精度的近似值xi才停止迭代，也就是绝对值|xmin－xi|足够小。

这里提到了标准MATLAB系统的两个最优化命令，fmin命令可以求单变量函数的最小值；fmins命令可以求多变量函数的最小值，同时它还要求有一个初始向量。

在新版本中，fmin和fmins分别被fminbnd和fminsearch取代。

Matlab中没有求函数f的最大值的命令，但可以通过求其相反函数h=－f 的最小值间接求得。 函数格式：

x=fmin ( fcn , x1 , x2 , options) 求函数在区间(x1，x2 )内取最小值时的x值，采用黄金分割法和二次插fcn是目标函数名。值。如果没有局部最小值，则返回区间内的最小x值。向量options为控制参数，如options(1) = 1，显示中间结果；options(2)表示得到的结果x的精度，缺省为10－4，options可省略。

x=fmins (fcn, x0, options) 求函数fcn的最小值。由用户自己给出一个初始估计向量x0，优化参数options可省略。

3

例：在区间[ 0，2 ]内求函数f(x)=x-2x-5的最小值。 f=inline('x.^3-2*x-5'); x = fmin(f, 0, 2) x =

0.8165 f(x) ans =

-6.0887

例：求香蕉函数f(x)?100(x2?x1)?(1?x1) 的最小值。 >> banana= inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2');

x = fmins(banana, [-1.2, 1]) x =

1 1

banana ([1,1]) ans = 0

此外，MATLAB提供了最优化工具箱。

222

48 第8章 符号计算

符号数学工具箱是操作和解决符号表达式的符号数学工具箱(函数)集合，有复合、简化、微分、积分以及求解代数方程和微分方程的工具。另外还有一些用于线性代数的工具，求解逆、行列式、正则型式的精确结果，找出符号矩阵的特征值而无由数值计算引入的误差。

符号数学工具箱中的工具是建立在功能强大的称作Maple软件的基础上。它最初是由加拿大的滑铁卢（Waterloo）大学开发的。当要求MATLAB进行符号运算时，它就请求Maple去计算并将结果返回到MATLAB命令窗口。因此，在MATLAB中的符号运算是MATLAB处理数字的自然扩展。

8.1 符号计算基础

一、 定义符号变量

参与符号运算的对象可以是符号变量、符号表达式或符号矩阵。符号变量要先定义，后引用。 1. sym函数：每次定义一个符号变量。 x=sym('x') %创建变量x。

x=sym('x', 'real') %创建实变量x。x=sym('x', 'unreal')去掉x的实数属性。 a=sym('alpha') %创建变量a，打印为alpha。 2. syms函数

sysms arg1 arg2 … argN。一次可定义多个符号变量，使用起来更方便。注意各变量间用空格隔开。

二、 符号表达式

符号表达式的书写格式与数值表达式一样，由符号变量、函数、运算符组成，如： syms a b c x

f=a*x^2 + b*x + c

2

也可以用f = sym('a*x^2 + b*x + c') 把符号表达式ax+bx+c赋给f，但此时没有定义变量a b c x。

三、 生成符号函数

有两种方法创建符号函数，一是用符号表达式，如上面f=a*x^2 + b*x + c 定义了符号函数f，从而可以用diff等符号函数对f进行运算。如diff(f)将返回结果2*a*x+b。

一是用M文件，如定义抽样函数：

function z = sinc(x)

if isequal(x,sym(0)) %数值不能直接参与符号运算，先转化为符号元素。 z = 1; else

z = sin(x)/x; end

四、 默认符号变量

当字符表达式中含有多于一个的变量，但其中只有一个变量是独立变量时，如果不告诉MATLAB哪

49 一个变量是独立变量，MATLAB将基于以下规则选择：选择在字母顺序中最接近x的字母，如果有两个距x相等的字母，就选择在字母表中较后的那一个。

可以利用find(f,1)找出符号函数f中的默认符号变量。如： syms a b x f=sin(a*x+b);

findsym(f,1) %返回ans = x

五、 符号矩阵

符号矩阵的元素为符号表达式，可用函数sym来产生符号矩阵。 例：>> G=sym( ' [cos(t)， sin(t)；-sin(t)，cos(t)] ' ) G= [ cos(t)， sin(t)] [-sin(t)， cos(t)]

由于数值不能直接参与符号运算，因此先要将数值矩阵转化为符号矩阵。函数sym也可以把数值矩阵转换成符号形式。如：

G=sym( '[2/3 sin(2); 0.33 1/0.7] ' ) G = [ 2/3, sin(2)]

[ 0.33, 1/0.7]

六、 简化表达式

1. simplify：Simplify是功能强大、通用的工具。它利用各种类型代数恒等式，包括求和、积分和分数幂、三角、指数和log函数、Bessel函数、超几何函数和 函数，来简化表达式。

例：syms x

S = (x^2+5*x+6)/(x+2); R = simplify(S) R =x+3

2. expand：将表达式按幂由低到高展开成多项式形式（也可以是三角、指数、对数形式）。 例：expand((x-2)*(x-4)) ans =x^2-6*x+8

3. factor：因式分解，如factor(x^3-y^3) 返回 (x-y)*(x^2+x*y+y^2) 4. collect：合并同类项，例： collect((exp(x)+x)*(x+2))

ans =x^2+(exp(x)+2)*x+2*exp(x)

5. numden：如果表达式是一个有理分式(两个多项式之比)，或是可以展开为有理分式(包括哪些分母为1的分式)，可利用numden来提取分子或分母

6. horner(s)：嵌套， 例：

horner(x^3-6*x^2+11*x-6) ans =-6+(11+(-6+x)*x)*x

7. simple：simple试用上面几种不同的简化方法，然后选择在结果表达式中含有最少字符的那种形式。

七、 可变精度运算

符号工具箱支持可变精度运算，即支持符号计算并能以指定的精度返回结果。 digits(5) %设置精度 vpa('pi') %返回3.1416

50 八、 函数计算器

函数 funtool

格式 funtool %该命令将生成三个图形窗口，Figure No.1用于显示函数f的图形，Figure No.2用于显示函数g的图形，Figure No.3为一可视化的、可操作与显示一元函数的计算器界面。在该界面上由许多按钮，可以显示两个由用户输入的函数的计算结果：加、乘、微分等。funtool还有一函数存储器，允许用户将函数存入，以便后面调用。在开始时，funtool显示两个函数f(x) = x与g(x) = 1在区间[-2*pi, 2*pi]上的图形。Funtool同时在下面显示一控制面板，允许用户对函数f、g进行保存、更正、重新输入、联合与转换等操作。

图 8-1 函数计算器

8.2 微积分

一、 极限

函数 limit

格式 limit(F,x,a) %计算符号表达式F=F(x)当x→a时的极限值。

limit(F,a) %用命令findsym(F)确定F中的自变量，设为变量x，再计算F的极限值，当x→a时。 limit(F) %用命令findsym(F)确定F中的自变量，设为变量x，再计算F的极限值，当x→0时。 limit(F,x,a,'right')或limit(F,x,a,'left') %计算符号函数F的单侧极限：左极限x→a- 或右极限x→a+。 例：

syms x a t h;

Then limit(sin(x)/x) => 1 limit(1/x,x,0,'right') => inf limit((sin(x+h)-sin(x))/h,h,0) => cos(x)

二、 导数（包括偏导数）

函数 diff

格式 diff(S,'v')、diff(S,sym('v')) %对表达式S中指定符号变量v计算S的1阶导数。 diff(S) %对表达式S中的符号变量v计算S的1阶导数，其中v=findsym(S)。 diff(S,n) %对表达式S中的符号变量v计算S的n阶导数，其中v=findsym(S)。 diff(S,'v',n) %对表达式S中指定的符号变量v计算S的n阶导数。 例：

syms x t

diff(sin(x^2)) %结果为 2*cos(x^2)*x

本文来源：https://www.bwwdw.com/article/pk5o.html