下三角矩阵的逆简单算法矩阵单独一行乘一个数变吗

发布时间：2023-06-16 作者：admin 来源：文学

形容快乐的成语-擂沙汤圆

2023年3月3日发(作者：宫东风)

基础篇

本篇包含五个线性代数的基础实验，从矩阵运算到方程组的求解；从向量组线性相关性

分析到矩阵的对角化；从矩阵特征值和特征向量求解到二次型的标准化及正定性的分析，都

给出了MATLAB的解决方法。实验5利用MATLAB的绘图功能，对线性代数若干概念的几何意

义进行了分析讨论。

实验1矩阵的基本运算

1.1实验目的

1.掌握Matlab软件的矩阵赋值方法；

2.掌握Matlab软件的矩阵加法、数乘、转置和乘法运算；

3.掌握Matlab软件的矩阵幂运算及逆运算；

4.掌握Matlab软件的矩阵元素群运算；

5.通过Matlab软件进一步理解和认识矩阵的运算规则。

1.2实验指导

MATLAB是一种功能强大的科学及工程计算软件，它的名字由“矩阵实验室”（Matrix

Laboratoy）组成，它具有以矩阵为基础的数学计算和分析功能，并且具有丰富的可视化图

形表现功能及方便的程序设计能力。它的应用领域极为广泛。本实验学习用MATLAB软件

进行矩阵基本运算。

启动MATLAB后，将显示MATLAB操作界面，它包含多个窗口，其中命令窗口是最

常用的窗口，如图1.1所示。

图1.1MATLAB的操作桌面

本实验所有例题的MATLAB命令都是在命令窗口中键入的。在本实验中用到MATLAB

的运算符号及命令或函数列举如下：

1、运算符号

表1.1给出了本实验用到的MATLAB基本运算符号。

表1.1MATLAB的基本运算符号

运算符号＝＋－

*^‘.

说明赋值加减乘左除右除幂运算转置群运算

2、命令或函数

表1.2给出了与本实验相关的MATLAB命令或函数。若要进一步了解和学习某个命令

或函数的详细功能和用法时，MATLAB提供了一个help命令。

表1.2与本实验相关的MATLAB命令或函数

命令说明位置

helpinv

在命令窗口中显示函数inv的帮助信息

[]

创建矩阵例1.1

，矩阵行元素分割符号例1.1

；矩阵列元素分割符号例1.1

注释行例1.1

eye(n)

创建n阶单位矩阵例1.1

zeros(m,n)

创建m×n阶零矩阵例1.1

zeros(n)

创建n阶零方阵例1.1

ones(m,n)

创建m×n阶元素全为1的矩阵例1.1

rand(m,n)

创建m×n阶元素为从0到1的均匀分布的随机数矩阵例1.2

round(A)

对矩阵A中所有元素进行四舍五入运算例1.2

inv(A)

求矩阵A的逆例1.3

A^-1

用幂运算求矩阵A的逆例1.3

1.3实验内容

例1.1用MATLAB软件生成以下矩阵：

（1）















066

656

239

A（2）















100

010

001

B（3）















C（4）















1111

解：（1）在MATLAB命令窗口输入：

A=[9,3,2;6,5,6;6,6,0]%矩阵同行元素以逗号或空格分割

或：A=[932;656;660]%行与行之间必须用分号或回车分隔

或：A=[932

656

660]

结果都为：

932

656

660

（2）输入：

B=eye(3)

结果为：

100

010

001

（3）输入：

C=zeros(2)

结果为：

（4）输入：

D=ones(4)

结果为：

1111

Matlab对矩阵赋值有直接输入和命令生成两种方法，本例中矩阵A就是键盘直接输入

的；而矩阵B、C和D都是用Matlab命令而生成。

例1.2随机生成两个3阶方阵A和B，分别计算：

（1）A＋B；（2）A－B；（3）5A；（4）AB；（5）TA

解：输入：

A=round(rand(3)*10)%rand(3)：生成3阶元素为0－1的随机实数方阵

%round（）：对矩阵元素进行四舍五入运算

B=round(rand(3)*10)

结果为：

1023

5109

937

123

035

971

（1）输入：

A＋B

结果为：

ans=

1146

51314

18108

其中“ans”表示这次运算的结果。

（2）输入：

A－B

结果为：

ans=

900

574

0-46

（3）输入：

5*A

结果为：

ans=

501015

255045

451535

（4）输入：

A*B

结果为：

ans=

374743

8610374

727649

（5）输入

A’

结果为

ans=

1059

2103

397

例1.3已知矩阵















712

010

321

A，分别计算：（1）5A；（2）1A

解：输入：

A=[1,2,3;0,1,0;2,1,7]

结果为：

123

010

217

（1）输入：

A^5

结果为：

ans=

010

781

（2）输入：

inv(A)

或输入A^-1

结果都为：

ans=

7-11-3

010

-231

例1.4已知矩阵















192

250

596

A，















182

401

266

B，且满足BPA，BAQ，计算矩阵

P和Q。

解：方法一：利用求逆矩阵的方法，输入：

A=[6,9,5;0,5,2;2,9,1]

B=[6,6,2;1,0,4;2,8,1]

P=B*inv(A)

Q=inv(A)*B

方法二：利用MATLAB软件特有的矩阵“左除”和“右除”运算，输入：

A=[6,9,5;0,5,2;2,9,1]

B=[6,6,2;1,0,4;2,8,1]

P=B/A%矩阵右除

Q=AB%矩阵左除

两种方法的运算结果都为：

695

052

291

662

104

281

0.8043-1.30430.5870

0.57611.1739-1.2283

0.0435-0.04350.8696

0.60871.4565-1.2065

0.04350.78260.2174

0.3913-1.95651.4565

例1.5已知矩阵















107

026

305

A，

















254

603

312

B，分别按以下要求进行矩阵元素的群

运算：

（1）把矩阵A和矩阵B所有对应元素相乘，得到9个乘积，计算由这9个数所构成的同形

矩阵C。

（2）对矩阵A中的所有元素进行平方运算，得到矩阵D，求该矩阵。

解：Matlab软件提供了矩阵元素群运算的功能，输入：

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

B=[2,1,3;3,0,6;4,5,-2]

结果为：

503

620

701

213

306

45-2

（1）输入：

C=A.*B%在运算符号前加“.”，其含义即为矩阵元素的群运算

结果为：

1009

1800

280-2

（2）输入：

D=A.^2%在运算符号前加“.”，其含义即为矩阵元素的群运算

结果为：

2509

3640

4901

1.4实验习题

1．利用函数rand和函数round构造一个5×5的随机正整数矩阵A和B，验证以下等式是

否成立：

（1）BAAB；（2）22BABABA；（3）TT

TABAB

2．已知向量4321，0107，请计算它们的内积。要求：

（1）用矩阵相乘命令计算；（2）用矩阵元素群运算的方法计算。

3．已知

































0152

095

038

125

231

135

X，用求逆矩阵和矩阵右除两种方法求矩阵X。

4、已知ABBA，其中

















200

012

021

B，用求逆矩阵和矩阵左除两种方法求矩阵A。

5、已知















43215

32154

21543

15432

54321

A，求5A，1A。

实验2行列式与方程组的求解

2.1实验目的

1.掌握Matlab软件求行列式的命令；

2.掌握Matlab软件求矩阵秩的命令；

3.掌握Matlab软件对矩阵进行初等行变换的命令；

4.掌握Matlab软件求解满秩线性方程组的各种方法；

5.掌握Matlab软件符号变量的应用；

6.通过Matlab软件验证与行列式相关的各种公式和定理，从而加深对相关概念的理解。

2.2实验指导

本实验利用MATLAB软件来计算与行列式相关的各种运算问题，其中包括：行列式的

求解、矩阵秩的求解、矩阵逆的求解、利用克莱姆法则解方程组、验证行列式按行（列）展

开定理及符号变量在行列式中的应用等。

MATLAB软件不仅拥有简单明了的命令窗口，而且也提供了程序编辑器。在实验1图

1.1所示的MATLAB操作界面中，点击左上方的小按钮，即可打开MATLAB的M文件编

辑器窗口，如图2.1所示，在这个窗口中可以编写扩展名为M的文件。

图2.1MATLAB的M文件编辑器

表2.1给出了与本实验相关的MATLAB命令或函数。

表2.1与本实验相关的MATLAB命令或函数

命令功能说明位置

A=[….];

在赋值语句后，若有一个分号“;”，它的含义是不在窗口中显示矩阵A。例2.1

U=rref(A)

对矩阵A进行初等行变换，矩阵U为矩阵A的最简行阶梯矩阵。例2.1

clear

清除工作空间中的各种变量，往往写在一个程序最前。例2.1

n=input(\'…\')

数据输入函数，单引号内的字符串起说明作用。例2.1

if…elseif…end

条件语句，控制程序流程，和C语言功能类似。例2.1

[m,n]=size(A)

计算结果为一个二维行向量，m,n分别存放矩阵A的行数和列数。例2.1

关系运算符号：等于例2.1

关系运算符号：不等于例2.1

逻辑运算符号：逻辑或例2.1

disp(‘…….’)

显示单引号中的字符串例2.1

det(A)

计算矩阵A的行列式例2.1

B(:,i)=b

把向量b赋给矩阵B的第i列，要求矩阵B的列向量和向量b同形。例2.1

[A,eye(5)]

创建一个5×10的矩阵，前5列为矩阵A，后5列为单位矩阵I例2.2

B(:,1:5)

取矩阵B的第1列至第5列例2.2

rank(A)

计算矩阵A的秩例2.2

for…end

for循环语句，控制程序流程，和C语言功能类似。例2.2

T(1,:)=[]

把一个空行[]赋给矩阵T的第1行，即：删除矩阵T的第一行。例2.2

A(i,j)

引用矩阵A中第i行第j列的元素例2.2

formatshort

以短格式显示例2.2

formatlong

以长格式显示例2.2

symsx

定义x为符号变量例2.3

factor(D)

对符号变量多项式D进行因式分解例2.3

solve(D)

求符号变量多项式方程D＝0的解例2.3

randn(m,n)

创建m×n阶均值为0，方差为1的标准正态分布的随机矩阵例2.4

2.3实验内容

例2.1已知非齐次线性方程组：





















85103537

2227772

9021161153

591310732

8054326

54321

xxxxx

，要求用下列方法

求解该方程组。

（1）求逆矩阵法；

（2）矩阵左除法；

（3）初等行变换；

（4）克莱姆法则。

解：（1）把齐次线性方程组写为矩阵形式：bAx，则bAx1，直接在MATLAB的命

令窗口输入：

A=[6,2,3,4,5;2,-3,7,10,13;3,5,11,-16,21;2,-7,7,7,2;7,3,-5,3,10];

%赋值语句最后的分号“;”，表示不在窗口中显示矩阵A

b=[80;59;90;22;85];

x=inv(A)*b

%或：x=A^-1*b

计算结果为：

9.0000

3.0000

2.0000

1.0000

2.0000

（2）矩阵的乘法不遵守乘法交换律，Matlab软件定义了矩阵左除和矩阵右除运算，针对方

程组的矩阵形式bAx，等式两端同时左除A，得到：“bAx”。在MATLAB命令窗

口中输入：

A=[6,2,3,4,5;2,-3,7,10,13;3,5,11,-16,21;2,-7,7,7,2;7,3,-5,3,10];

b=[80;59;90;22;85];

x=Ab%符号“”即为左除运算，注意它的方向。

结果为：

9.0000

3.0000

2.0000

1.0000

2.0000

（3）用初等行变换，把方程组的增广矩阵变换为最简行阶梯形式，从而得到方程组的解。

在MATLAB命令窗口中输入：

A=[6,2,3,4,5;2,-3,7,10,13;3,5,11,-16,21;2,-7,7,7,2;7,3,-5,3,10];

b=[80;59;90;22;85];

U=rref([A,b])%对增广矩阵[A,b]进行行初等变换

%矩阵U为矩阵A的最简行阶梯形矩阵

运算结果为：

100009

010003

001002

000101

000012

矩阵U即为方程组增广矩阵的最简行阶梯形矩阵，矩阵的最后一列即为方程组的解。

（4）根据克莱姆法则，有：

，其中D是方程组的系数行列式，

D是用常数列向

量b代替系数行列式的第i列所得到的行列式。用Matlab的M文件编辑器，编写la01.m文

件如下：

%用克莱姆法则求解方程组

clear%清除变量

n=input(\'方程个数n＝\')%请用户输入方程个数

A=input(\'系数矩阵A=\')%请用户输入方程组的系数矩阵

b=input(\'常数列向量b=\')%请用户输入常数列向量

if(size(A)~=[n,n])|(size(b)~=[n,1])

%判断矩阵A和向量b输入格式是否正确

disp(\'输入不正确,要求A是n阶方阵，b是n维列向量\')

%disp：显示字符串

elseifdet(A)==0%判断系数行列式是否为零

disp(\'系数行列式为零，不能用克莱姆法则解此方程。\')

else

fori=1:n%计算x1,x2,...xn

B=A;%构造与A相等的矩阵B

B(:,i)=b;%用列向量b替代矩阵B中的第i列

x(i)=det(B)/det(A);%根据克莱姆法则计算x1,x2,...xn

end

x=x\'%以列向量形式显示方程组的解

end

在MATLAB命令窗口中输入：

la01

得到以下人机对话结果：

方程个数n＝5

系数矩阵A=

[6,2,3,4,5;2,-3,7,10,13;3,5,11,-16,21;2,-7,7,7,2;7,3,-5,3,10]

62345

2-371013

3511-1621

2-7772

73-5310

常数列向量b=[80;59;90;22;85]

显然，当方程组的系数矩阵不是方阵，或系数行列式等于零时，逆矩阵法和克莱姆法都

不能实现方程组的求解，而初等行变换的方法适合各种线性方程组的求解，在下一个实验中

将继续讨论rref命令的详细应用。关于矩阵的左除运算，有着多种的含义，在实验3和实验

5中将逐步讨论它的其它功能。

例2.2求矩阵

























41118307

39203

259113

113631

64627

A的逆，要求用以下方法：

（1）矩阵左除和右除运算；

（2）初等行变换；

（3）利用伴随矩阵*A求逆的公式

1。

解：在实验1中，已经介绍了两种求矩阵逆的方法：一种是命令法：inv（A），一种是幂运

算法：A^-1。下面分别给出求逆矩阵的其它方法。

（1）Matlab软件定义了矩阵的左除运算和右除运算，给矩阵运算带来了很大的方便。从公

式IAA1出发，等式两边同时右除矩阵A，得“AIA/1”，从而求得矩阵A的逆阵；

若从公式IAA1出发，等式两边同时左除矩阵A，得“IAA1”，同样可以求得矩

阵A的逆阵。需要注意的是，左除和右除不仅是矩阵的左右位置不同，其运算符号也是不

同的。

（2）用笔来计算一个具体n阶方阵A的逆阵时，往往采用初等行变换法，即：是把n阶矩

阵A和n阶单位矩阵I并排放在一起，然后对n×2n矩阵IA|做初等行变换，当其变为

最简行阶梯矩阵时，如果该矩阵的前n列为单位矩阵I，则该矩阵的后n列就是矩阵A的

逆矩阵1A；如果该矩阵的前n列不是单位矩阵I，则矩阵A不可逆。

（3）根据IAAAAA**，可以得到求逆矩阵公式：

1，其中*A为矩阵A的

伴随矩阵。

在MATLAB的M文件编辑器中，编写程序la02.m：

%逆矩阵各种求法：

clear

A=[-7,-2,-6,4,6;1,3,-6,3,11;3,-11,9,5,-2;-3,0,-2,9,-3;7,30,-18,11

,4];

%1.命令法：

An1=inv(A)

%2.幂运算法：

An2=A^-1

%3.右除法：

An3=eye(5)/A%eye(5)为5阶单位矩阵

%4.左除法：

An4=Aeye(5)

%5.初等行变换法：

B=rref([A,eye(5)]);%对矩阵[A,I]进行初等行变换

%B为矩阵A的最简行阶梯矩阵

if(rank(B(:,1:5))==5)%判断最简行阶梯矩阵B的前5列是否为单位阵

An5=B(:,6:10)%取出矩阵的后5列，并显示

else

disp(\'A不可逆\');

end

%6.伴随矩阵求逆法：

fori=1:5%构造伴随矩阵的5×5个元素

forj=1:5

T=A;%把矩阵A赋给矩阵T

T(i,:)=[];%删去矩阵T的第i行

T(:,j)=[];%删去矩阵T的第j列

%此时，|T|为矩阵A元素aij的余子式

AA(j,i)=(-1)^(i+j)*det(T);

%算出aij的代数余子式

%并放入矩阵AA的第j行、第i列

%当循环结束，矩阵AA即为A的伴随矩阵

end

ifdet(A)~=0

An6=AA/det(A)

else

disp(\'A不可逆\');

end

运算程序la02，前四个方法计算结果相同：

1.0e+004*

-1.58951.3448-1.06461.6206-0.6308

1.6298-1.37891.0916-1.66170.6468

2.5392-2.14831.7007-2.58891.0077

0.3631-0.30720.2432-0.37020.1441

0.9860-0.83420.6604-1.00530.3913

后两个方法计算结果相同：

-1589513448-1064616206-6308

16298-1378918

25392-214831777

3631-30722432-37021441

9860-83426604-100533913

从计算结果可以发现，前四个方法得到的是实数矩阵，而后两个方法得到的是整数矩阵。

如果在Matlab环境下，键入：formatlong（以上结果是在formatshort格式下得到的），然后

再重新运行该程序，会发现前四个方法的运算结果存在误差，这是计算机做数值运算时，存

在舍入误差的原因。为了进一步观察计算机做数值运算所产生的误差，现在用上述六种方法

来计算矩阵















131211

101010

321

A的逆，前四种方法得到以下类似结果：

Warning:Matrixisclosetosingularorbadlyscaled.

=2.135044e-018.

ans=

1.0e+015*

-4.5036-4.50364.5036

9.00729.0072-9.0072

-4.5036-4.50364.5036

显然此结果是不正确的，因为A不可逆。

例2.3解方程：0

2317

2315

1223

1123





。

解：Matlab软件定义了“符号变量”的概念。在MATLAB的M文件编辑器中，应用“符

号变量”编写程序la03.m：

%求解符号行列式方程

clearall%清除各种变量

symsx%定义x为符号变量

A=[3,2,1,1;3,2,2-x^2,1;5,1,3,2;7-x^2,1,3,2]

%给矩阵A赋值

D=det(A)%计算含符号变量矩阵A的行列式D

f=factor(D)%对行列式D进行因式分解

%从因式分解的结果，可以看出方程的解

X=solve(D)%求方程“D＝0”的解

在MATLAB的命令窗口输入：

la03

运行结果为：

[3,2,1,1]

[3,2,2-x^2,1]

[5,1,3,2]

[7-x^2,1,3,2]

-6+9*x^2-3*x^4

-3*(x-1)*(x+1)*(x^2-2)

[1]

[-1]

[2^(1/2)]

[-2^(1/2)]

向量x即为方程的解，MATLAB针对符号变量可以得出解析解。

例2.4请用Matlab软件验证行列式按行（列）展开公式：











jkikji

jiA

当

解：用Matlab程序构造了一个5阶随机数方阵A。首先，按第一行展开：

AaAaAas，验证s是否与A的行列式相等；其次，计算A的第一行

元素与第三行元素对应的代数余子式乘积之和：

351532123111

AaAaAas，验算s是

否为0。在MATLAB的M文件编辑器中，编写程序la04.m：

%验证行列式按行（列）展开公式

clear

A=round(10*randn(5));%构造5阶随机数方阵

D=det(A);%计算矩阵A的行列式

%矩阵A按第一行元素展开：s=a11*A11+a12*A12+…+a15*A15

s=0;

fori=1:5

T=A;

T(1,:)=[];%删去阵矩第1行

T(:,i)=[];%删去矩阵第i列

%此时，|T|为矩阵A元素a1i的余子式

s=s+A(1,i)*(-1)^(1+i)*det(T);

end

e=D-s%验算D与s是否相等

在MATLAB的命令窗口中输入：

la04

计算结果为：

在MATLAB的M文件编辑器中，编写程序la05.m：

%计算5阶方阵A的第一行元素与第三行元素对应的代数余子式乘积之和：

%s=a11*A31+a12*A32+…+a15*A35

clear

A=round(10*randn(5));%构造5阶随机数方阵

s=0;

fori=1:5

T=A;

T(3,:)=[];%删去矩阵第3行

T(:,i)=[];%删去矩阵第i列

%此时，|T|为矩阵A元素a3i的余子式

s=s+A(1,i)*(-1)^(3+i)*det(T);

end

s%验算s是否为0

在MATLAB命令窗口中输入：

la05

计算结果为：

2.4实验习题

1．利用函数randn和函数round构造一个6×6的随机整数矩阵A，验证矩阵A的行列式满

足下列性质：

（1）TAB，验证AB；

（2）把矩阵A的第二行和第五行进行对调后的矩阵为B，验证AB

（3）把矩阵A的第三列加到第一列后的矩阵为B，验证AB

2．构造5阶方阵A，验证公式：IAAAAA**，其中*A为矩阵A的伴随矩阵，I为

同阶单位阵。

3.求解非齐次线性方程组



















1201099

560074

21292954

11414622

54321

xxxxx

。

4．求下列矩阵的逆：

（1）























931154

16142022

11713511

26112176

55201610

（2）















101101012

00071

112306

1535182

107135

5．求下列含符号变量的行列式，并要求把结果因式分解。

（1）





11000

1100

0110

0011

0001

；（2）

abbb

abaa

aaba

bbba

；（3）

4444

2222

1111

dcba

实验3向量组的相关性及方程组的通解

3.1实验目的

1.掌握Matlab软件分析向量组线性相关性的方法；

2.掌握Matlab软件求解线性方程组通解的各种方法；

3.通过Matlab软件进一步理解和认识齐次线性方程组解空间的概念。

3.2实验指导

本实验利用MATLAB软件来分析向量组的线性相关性；向量的线性表示；齐次线性方

程组的通解及非齐次线性方程组的通解。

表3.1给出了与本实验相关的MATLAB命令或函数。

表3.1与本实验相关的MATLAB命令或函数

命令功能说明位置

[R,s]=rref(A)

把矩阵A的最简行阶梯矩阵赋给R；s是一个行向量，它的元素由R的基

准元素所在的列号构成。

例3.1

length(s)

计算向量s的长度，即向量s的维数。例3.1

end

矩阵的最大下标，即：最后一行或最后一列。例3.1

null(A,\'r\')

计算齐次线性方程组Ax＝0的基础解系例3.1

x0=Ab

求非齐次线性方程组Ax＝b的一个特解x0例3.1

fprintf

按指定格式写文件，和C语言功能类似。例3.2

find(s)

计算向量s中非零元素的下标例3.2

subs(A,k,n)

将A中的所有符号变量k用数值n来替代例3.3

3.3实验内容

例3.1求非齐次线性方程组



















4319252

23196463

7236263

216442

54321

xxxxx

的通解。

解：在MATLAB命令窗口，输入以下命令：

A=[2,4,-1,4,16;-3,-6,2,-6,-23;3,6,-4,6,19;1,2,5,2,19];

%输入系数矩阵A

b=[-2;7;-23;43];%输入常数列向量b

[R,s]=rref([A,b])%把增广矩阵的最简行阶梯矩阵赋给R

%而R的所有基准元素在矩阵中的列号构成了行向量s

计算结果为：

120293

001028

000000

在得出该方程组增广矩阵的最简行阶梯矩阵R后，根据线性代数知识可以得到该齐次

线性方程组的通解。下面，在MATLAB的M文件编辑器中，编写程序la06.m，可以给出该

方程组的一个具体特解和对应齐次方程组的通解。

%求齐次线性方程组的通解

clear

A=[2,4,-1,4,16;-3,-6,2,-6,-23;3,6,-4,6,19;1,2,5,2,19];

%输入系数矩阵A

b=[-2;7;-23;43];%输入常数列向量b

[R,s]=rref([A,b]);%把增广矩阵的最简行阶梯矩阵赋给R

%而R的所有基准元素在矩阵中的列号构成了行向量s

[m,n]=size(A);%矩阵A的行数、列数赋给了变量m、n

x0=zeros(n,1);%将特解x0初始化为N维零列向量

r=length(s);%矩阵A的秩赋给变量r

x0(s,:)=R(1:r,end);%将矩阵R的最后一列按基准元素的位置给特解x0赋值

disp(\'非齐次线性方程组的特解为：\')

x0%显示特解x0

disp(\'对应齐次线性方程组的基础解系为：\')

x=null(A,\'r\')%得到齐次线性方程组Ax＝0的基础解系x

在MATLAB命令窗口中输入：

la06

运算结果为：

非齐次线性方程组的特解为：

x0=

对应齐次线性方程组的基础解系为：

-2-2-9

100

00-2

010

001

则方程组的通解为：





























































321

kkk

此计算方法和与传统的笔算方法一致，所以其结果也是一致的。齐次线性方程组的特解

还可以用Matlab的矩阵左除运算来求得，直接在MATLAB命令窗口输入以下命令：

A=[2,4,-1,4,16;-3,-6,2,-6,-23;3,6,-4,6,19;1,2,5,2,19];

b=[-2;7;-23;43];

x0=Ab%用矩阵左除运算求得方程组特解x0

x=null(A,\'r\')%得到齐次线性方程组Ax＝0的基础解系x

运算结果为：

Warning:Rankdeficient,rank=2tol=4.3099e-014.

x0=

7.3333

0.3333

-2-2-9

100

00-2

010

001

其中特解x0与前一方法的特解不同。（注：如果欠定方程组有解，则它有无穷个特解，通解

中只需要任何一个特解即可）

方程组的通解为：





























































322

321

kkk

例3.2已知向量组















，















，















，















，

















，求出它的

最大无关组，并用该最大无关组来线性表示其它向量。

解：用笔计算的过程为：



















22132

91682

22000

23341

69231

,,,,

54321

A，对矩阵A进行初等行变换，最后变为最简

行阶梯矩阵：















00000

11000

20110

30101

，该矩阵基准元素所在的列号为1，2，4，则原向量组的一个最大无

关组为：

421

,,。根据该矩阵的第3列，可以得到：

4213

011；同理，

可以得到：

4215

123。

根据以上思路，编写Matlab程序la07.m：

%找向量组的最大无关组，并用它线性表示其它向量

clear

a1=[1;1;0;2;2];%输入5个列向量

a2=[3;4;0;8;3];

a3=[2;3;0;6;1];

a4=[9;3;2;1;2];

a5=[6;-2;2;-9;2];

A=[a1,a2,a3,a4,a5];%由5个列向量构造矩阵A

[R,s]=rref(A);%把矩阵A的最简行阶梯矩阵赋给了R

%而R的所有基准元素在矩阵中的列号构成了行向量s

%向量s中的元素即为最大无关组向量的下标

r=length(s);%最大无关组所含向量个数赋给r

fprintf(\'最大线性无关组为：\')%输出字符串

fori=1:r

fprintf(\'a%d\',s(i))%分别输出最大无关组的向量a1，…

end

fori=1:r%从矩阵A中取出最大无关组赋给A0

A0(:,i)=A(:,s(i));

end

A0%显示最大无关组矩阵A0

s0=[1,2,3,4,5];%构造行向量s0

fori=1:r

s0(s(i))=0;%s(i)是最大无关组的列号

end%若s0的某元素不为0，表示该元素为矩阵A中

%除最大无关组以外其它列向量的列号

s0=find(s0);%删除s0中的零元素

%此时s0中元素为其它向量的列号

fori=1:5-r%用最大无关组来线性表示其它向量

fprintf(\'a%d=\',s0(i))

forj=1:r

fprintf(\'%3d*a%d+\',R(j,s0(i)),s(j));

end

fprintf(\'bbn\');%去掉最后一个”+”

end

在MATLAB命令窗口中输入：

la07

运行结果为：

最大线性无关组为：a1a2a4

A0=

139

143

002

281

232

a3=-1*a1+1*a2+0*a4

a5=3*a1+-2*a2+1*a4

例3.3已知齐次线性方程组：























011333

03233

03323

033321

4321

xkxxx

xxkxx

xxxkx

xxxxk

，当k取何值时方程组有非

零解？在有非零解的情况下，求出其基础解系。

解：在MATLAB的M文件编辑器中，编写程序la08.m

%计算带符号变量的齐次线性方程组的解

clear

symsk%定义符号变量k

A=[1-2*k,3,3,3;3,2-k,3,3;3,3,2-k,3;3,3,3,11-k];%给系数矩阵赋值

D=det(A);%算出系数矩阵的行列式D

kk=solve(D);%解方程“D＝0”，得到解kk，即k值

fori=1:4

AA=subs(A,k,kk(i));%分别把k值代入系数矩阵A中

fprintf(\'当k=\');

disp(kk(i));%显示k的取值

fprintf(\'基础解系为：n\');

disp(null(AA))%计算齐次线性方程组“Ax=0”的基础解系

end

在MATLAB命令窗口中输入：

la08

运算结果为：

当k=7/2

基础解系为：

[1]

[2]

[-2]

当k=14

基础解系为：

[1]

[2]

[5]

当k=-1

基础解系为：

[-1,-1]

[1,0]

[0,1]

[0,0]

当k=-1

基础解系为：

[-1,-1]

[1,0]

[0,1]

[0,0]

3.4实验习题

1.求下列向量组的一个最大无关组，并把其余向量用此最大无关组线性表示。















，















，

















，

















，















，

















。

2.求非齐次线性方程组













12372

02324

43224

543236

54321

xxxxx

的通解。

3．已知齐次线性方程组：























02870

04523

0032

04422

4321

xkxxx

xxkxx

xxxkx

xxxxk

，当k取何值时方程组有非

零解？在有非零解的情况下，求出其基础解系。

实验4特征向量与二次型

4.1实验目的

1.掌握Matlab软件对向量组正交化的方法；

2.掌握Matlab软件求方阵特征值和特征向量的方法；

3.掌握Matlab软件分析方阵是否可对角化的方法；

4.掌握Matlab软件化二次型为标准型的方法；

5.掌握Matlab软件分析对称阵是否正定的方法；

6.了解MATLAB软件关于矩阵分解的命令。

4.2实验指导

本实验利用MATLAB软件来求解矩阵的特征值和特征向量；分析矩阵的对角化问题；

给出二次型标准化的方法；分析二次型的正定性；了解矩阵的分解命令。表4.1给出了与本

实验相关的MATLAB命令或函数。

表4.1与本实验相关的MATLAB命令或函数

命令功能说明位置

orth(A)

求出矩阵A的列向量组构成空间的一个正交规范基例4.1

P=poly(A)

计算矩阵A的特征多项式，P是一个行向量，其元素是多项式系数。例4.2

roots(P)

求该多项式P的零点例4.2

r=eig(A)

r为一列向量，其元素为矩阵A的特征值。例4.2

[V,D]=eig(A)

矩阵D为矩阵A特征值所构成的对角阵，矩阵V的列为矩阵A的单位特

征向量，它与D中的特征值一一对应。

例4.2

eval(lamda)

把符号形式转换为数值形式例4.2

[V,D]=schur(A)

同上例4.4

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

U、V都是正交矩阵，S是矩阵A的奇异值构成的对角矩阵，满足：A＝USVT例4.6

[L,U]=lu(A)

L为准下三角矩阵，U为上三角矩阵，满足：A＝LU例4.6

[Q,R]=qr(A)

Q为正交矩阵，R为上三角矩阵，满足：A＝QR例4.6

L=chol(A)

L为上三角矩阵，满足：A＝LTL（要求矩阵A为对称正定阵）例4.6

4.3实验内容

例4.1设向量组：















，

















，















，求由这三个向量生成的子空间V的

一个标准正交基。

解：在MATLAB命令窗口输入：

a1=[1;2;3];a2=[-1;1;2];a3=[5;1;0];

A=[a1,a2,a3]

P=orth(A)%将矩阵A的列向量组正交规范化，

%P的列构成了空间V的一个标准正交基

%P的列数反应了空间V的维数

运算结果为：

-0.92660.3359

-0.3116-0.4343

-0.2105-0.8358

从矩阵P的列数可以看出，原向量组是线性相关的，它生成的空间是二维的。

例4.2求矩阵





















7480

1648

119162

192022

A的特征值。

解：在MATLAB的M文件编辑器中，编写la09.m文件，它给出了三种求矩阵特征值的方

法：

%矩阵特征值的求解方法

clear

A=[2,-2,-20,-19;-2,16,-9,11;-8,4,-6,1;0,-8,-4,-7];

%方法一：

symsk%定义符号变量k

B=A-k*eye(length(A));%构造矩阵B=(A-kI）

D=det(B);%计算行列式：|A-kI|

lamda1=solve(D)%求|A-kI|=0的符号形式的解

%方法二：

P=poly(A);%计算矩阵A的特征多项式

%向量P的元素为该多项式的系数

lamda2=roots(P)%求该多项式的零点，即特征值

%方法三：

lamda3=eig(A)%直接求出矩阵A的特征值

在MATLAB命令窗口中输入：

la09

运算结果为：

lamda1=

[12]

[-1/3*(19585+120*22674^(1/2))^(1/3)-385/3/(19585+120*22674^(1/2))^(1/3)-7/3]

[1/6*(19585+120*22674^(1/2))^(1/3)+385/6/(19585+120*22674^(1/2))^(1/3)-7/3+1/2*i*3^(1/2)*

(-1/3*(19585+120*22674^(1/2))^(1/3)+385/3/(19585+120*22674^(1/2))^(1/3))]

[1/6*(19585+120*22674^(1/2))^(1/3)+385/6/(19585+120*22674^(1/2))^(1/3)-7/3-1/2*i*3^(1/2)*

(-1/3*(19585+120*22674^(1/2))^(1/3)+385/3/(19585+120*22674^(1/2))^(1/3))]

lamda2=

-17.3347

12.0000

5.1673+6.3598i

5.1673-6.3598i

lamda3=

-17.3347

5.1673+6.3598i

5.1673-6.3598i

12.0000

其中，方法一是根据笔算矩阵特征值的算法编写而成，MATLAB给出了一个符号形式

的解，可以进一步把符号解转化为数值解，输入以下命令：

lamda1=eval(lamda1)

结果为：

lamda1=

12.0000

-17.3347

5.1673-6.3598i

5.1673+6.3598i

例4.3求下列矩阵A的特征值和特征向量，并判断矩阵是否可以对角化，若能对角化，请

找出可逆矩阵V，使1VAVD。

(1)















211

312

321

A；（2）



















012

212

214

A；（3）

















101

034

011

解：（1）在MATLAB命令窗口输入：

A=[1,2,3;2,1,3;1,1,2];

[V,D]=eig(A)%矩阵D为矩阵A的特征值构成的对角阵

%矩阵V的列向量为矩阵A与特征值D对应的特征向量

运行结果为：

-0.6396-0.7071-0.5774

-0.63960.7071-0.5774

-0.42640.00000.5774

5.000000

0-1.00000

000.0000

从矩阵D可以看出，矩阵A有三个不同的特征值5，-1，0，则矩阵V的3个列向量线

性无关，所以矩阵A可以对角化。可以在MATLAB窗口中，验算结果1VAVD的正确

性。

（2）在MATLAB命令窗口输入：

A=[4,-1,-2;2,1,-2;2,-1,0];

[V,D]=eig(A)

运行结果为：

0.7276-0.57740.7437

0.4851-0.57740.2373

0.4851-0.57740.6250

2.000000

01.00000

002.0000

从矩阵D可以看出，2

是矩阵A的二重特征根，属于这个特征根的特征向量

是否存在两个线性无关的特征向量，需要分析齐次线性方程组02xIA解空间的维数

是否等于2，即分析3－rank2AI是否等于2。在MATLAB命令窗口继续输入：

if3-rank(A-2*eye(3))==2

%判断齐次线性方程组(A2－2I)x=0解空间的维数是否为2

disp(\'能对角化\');

%如果解空间维数为2，则存在2个线性无关特征向量

else

disp(\'不能对角化\');

end

运行结果为：

能对角化

从对角矩阵D和特征向量构成的矩阵V中也可以看出，矩阵V的第一列和第三列是属

于特征值2的特征向量，它们是线性无关的，则矩阵V的3个列向量线性无关，故矩阵A

可以对角化。在MATLAB命令窗口中，可以验证结果1VAVD的正确性。

（3）在MATLAB命令窗口输入：

A=[1,-1,0;4,-3,0;1,0,1];

[V,D]=eig(A)

运行结果为：

00.4364-0.4364

00.8729-0.8729

1.0000-0.21820.2182

1.000000

0-1.00000

00-1.0000

从矩阵D可以看出，1

是矩阵A的二重特征根，在MATLAB命令窗口继

续输入：

if3-rank(A+eye(3))==2

%判断齐次线性方程组(A2－2I)x=0解空间的维数是否为2

disp(\'能对角化\');

%如果解空间维数为2，则存在2个线性无关特征向量

else

disp(\'不能对角化\');

end

运行结果为：

不能对角化

其实，从对角阵D和特征向量构成的矩阵V中也可以看出，矩阵V的第二列和第三列

是属于特征值-1的特征向量，它们是线性相关的，即属于特征值-1的线性无关的特征向量

只有一个。故矩阵A不能对角化。

例4.4用正交变换法将二次型

1321

422,,xxxxxxxxf化为标准形。

解：在MATLAB的M文件编辑器中，编写la10.m文件：

%用正交变换法将二次型化为标准型

clear

A=[1,0,0;0,2,2;0,2,2];%输入二次型的矩阵A

[V,D]=eig(A);%其中矩阵V即为所求正交矩阵

%矩阵D为矩阵A的特征值构成的对角阵

%或：[V,D]=schur(A)%结果和eig()函数相同

disp(\'正交矩阵为：\');

disp(\'对角矩阵为：\');

disp(\'标准化的二次型为：\');

symsy1y2y3

f=[y1,y2,y3]*D*[y1;y2;y3]

在MATLAB命令窗口中输入：

la10

运行结果为：

正交矩阵为：

01.00000

-0.707100.7071

0.707100.7071

对角矩阵为为：

000

010

004

标准化的二次型为：

y2^2+4*y3^2

矩阵V为所求正交矩阵，把Vyx代入二次型

AxxfT，得：

2

4yyyyyAVVyAxxfTTTT。

例4.5判断下列矩阵的正定性：

（1）



















511

121

111

A；（2）

















513

122

321

B；

（3）



















541

452

121

C；（4）

















301

032

123

解：在MATLAB命令窗口输入：

A=[1,1,-1;1,2,-1;-1,-1,5];

B=[1,2,3;2,2,-1;3,-1,5];

C=[-3,2,1;2,-3,0;1,0,-3];

D=[1,2,-1;2,5,-4;-1,-4,5];

lamda_A=eig(A)

lamda_B=eig(B)

lamda_C=eig(C)

lamda_D=eig(D)

运行结果为：

lamda_A=

0.3542

2.0000

5.6458

lamda_B=

-1.8900

3.2835

6.6065

lamda_C=

-5.2361

-3.0000

-0.7639

lamda_D=

-0.0000

1.4689

9.5311

从矩阵特征值的正负，可以看出矩阵A正定，矩阵B不定，矩阵C负定，而矩阵D的

第一个特征值显示为“-0.0000”，继续查看为“-1.7351e-017”，那么它是等于零？还是小于

零呢？现在，利用例4.2中求特征值的方法一，再次求解矩阵D的特征值，在MATLAB命

令窗口输入：

symsk

D=[1,2,-1;2,5,-4;-1,-4,5];

d=det(D-k*eye(3));%计算行列式：|D-kI|

lamda_D=solve(d)%求特征方程|D-kI|=0的解

显示结果为：

lamda_D=

[0]

[11/2+1/2*65^(1/2)]

[11/2-1/2*65^(1/2)]

从以上符号解中可以看出：矩阵D的确有一个特征值为零，则它是半正定的。而计算

机在执行Matlab函数eig的运算中，产生了舍入误差。

以下是用求特征值符号形式的方法来判断对称阵正定性的通用程序la11.m：

%判断对称阵的正定性

clear

A=input(\'输入对称阵A:\');

ifA\'-A~=0%若矩阵不是对称阵，则退出

disp(\'输入错误\');

return;%退出该程序

end

n=length(A);%取矩阵A的阶数

symsk%定义k为符号变量

d=det(A-k*eye(n));%d为矩阵A的特征多项式

lamda=solve(d);%lamda为矩阵A的符号形式的特征根

lamda=eval(lamda);%把符号形式变为数值形式

iflamda>0%判断特征值是否全部大于零

disp(\'矩阵A为正定矩阵\');

elseiflamda>=0%判断特征值是否全部大于等于零

disp(\'矩阵A为半正定矩阵\');

elseiflamda<0%判断特征值是否全部小于零

disp(\'矩阵A为负定矩阵\');

elseiflamda<=0%判断特征值是否全部小于等于零

disp(\'矩阵A为半负定矩阵\');

else

disp(\'矩阵A为不定矩阵\');

end

在MATLAB命令窗口中输入：

la11

人机对话结果为：

输入对称阵A:[1,2,-1;2,5,-4;-1,-4,5]

矩阵A为半正定矩阵

例4.6已知对称矩阵















923

232

322

A，请按下列要求对矩阵A进行分解。

（1）特征值分解：找出正交矩阵V，使得：TVDVA，其中D为矩阵A的特征值构成

的对角阵；

（2）SVD分解：找出正交矩阵VU,，使得：TUSVA，其中S为矩阵A的奇异值构成

的对角阵；

（3）LU分解：找出一个准下三角矩阵L和一个上三角矩阵U，使得：LUA；

（4）QR分解：找出一个正交矩阵Q和一个上三角矩阵R，使得：QRA；

（5）Cholesky分解：找出一个上三角矩阵L，使得：LLAT；

解：在MATLAB命令窗口输入：

A=[2,2,3;2,3,2;3,2,9];

（1）特征值分解：输入：

[V,D]=eig(A)

%或[V,D]=schur(A)

结果为：

0.85510.36950.3637

-0.48560.81650.3121

-0.1817-0.44350.8777

0.226700

02.81870

0010.9546

（2）SVD分解：输入：

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

结果为：

-0.3637-0.3695-0.8551

-0.3121-0.81650.4856

-0.87770.44350.1817

10.954600

02.81870

000.2267

-0.3637-0.3695-0.8551

-0.3121-0.81650.4856

-0.87770.44350.1817

（3）LU分解：输入：

[L,U]=lu(A)

结果为：

0.66670.40001.0000

0.66671.00000

1.000000

3.00002.00009.0000

01.6667-4.0000

00-1.4000

把矩阵L的第一行与第三行交换，就变为了一个下三角矩阵。

（4）QR分解：输入：

[Q,R]=qr(A)

结果为：

-0.4851-0.0844-0.8704

-0.4851-0.80220.3482

-0.72760.59110.3482

-4.1231-3.8806-8.9738

0-1.39333.4620

001.2185

（5）Cholesky分解：输入：

L=chol(A)

结果为：

1.41421.41422.1213

01.0000-1.0000

001.8708

4.4实验习题

1．已知矩阵





















1102

121

2101

A，求下列矩阵的特征值和特征向量。

（1）A；（2）IAA32523；（3）132AI

2．把下面向量组正交化。

















，

















，



















3．判断下列矩阵，是否可以对角化，若能对角化，请找可逆矩阵V，使1VAVD。

（1）

















402

121

301

A；（2）

















2399

1611

0011

A；

4．用正交变换法将下列二次型化为标准形。

（1）

323121

1321

48444,,xxxxxxxxxxxxf

（2）

434232413121

14321

222222,,,xxxxxxxxxxxxxxxxxxxxf

5．判断下列矩阵的正定性。

（1）

















333

382

328

A；（2）

















212

122

221

A；（3）



















431

341

112

6．已知对称矩阵















631

321

111

A，请对矩阵A进行特征值分解；SVD分解；LU分解；QR

分解及Cholesky分解。并验证分解结果的正确性。

实验5线性代数的几何概念与Matlab作图

5.1实验目的

1.掌握Matlab软件制作二维和三维图形的基本方法；

2.利用Matlab软件的绘图功能，理解方程组解的几何意义；

3.利用Matlab软件的绘图功能，理解向量组线性相关的几何意义；

4.利用Matlab软件的绘图功能，理解二维向量线性变换的意义；

5.利用Matlab软件的绘图功能，理解二阶方阵特征值的几何意义；

6.利用Matlab软件的绘图功能，理解二次型正定性的几何意义。

5.2实验指导

本实验介绍了MATLAB的基本作图命令，通过二维和三维图形进一步来理解线性代数

中若干概念的几何意义，其中包括：线性方程组的解、向量的线性表示、线性变换、特征值

和特征向量及二次型的正定性等。表5.1给出了与本实验相关的MATLAB命令或函数。

表5.1与本实验相关的MATLAB命令或函数

命令功能说明位置

closeall

关闭所有图形例5.1

subplot(2,2,1)

准备画2×2个图形中的第一个图形例5.1

ezplot(\'x1+2*x2=5\')

绘制符号变量构成的直线方程x1+2*x2=5例5.1

holdon

准备在同一图中画多条直线，保留当前图形例5.1

title(\'\')

把单引号中字符串内容作为标题在图上方显示例5.1

gridon

在图中显示网格例5.1

x=Ab

求超定方程组Ax=b的最小二乘解例5.1

x=pinv(A)*b

求超定方程组Ax=b的最小二乘解，pinv(A)为对矩阵A进行

伪逆（广义逆）运算。

例5.1

plot(x,y,’*’)

在（x，y）位置画“*”，若x，y是向量，则画出一根曲线图。例5.1

ezmesh(\'x1+5*x2+1\')

绘制符号变量构成的平面方程：x3=x1+5*x2+1（即：

x1+5*x2-x3=-1）

例5.2

functiondrawvec(u)

定义函数drawvec，其参数为向量u，与C语言函数概念类

似。

例5.3

acos()

反余弦函数例5.3

norm(u)

计算矢量的范数，即矢量u的长度。例5.3

MATLAB定义的常数：圆周率例5.3

fill(x，y)填充二维多边形，多边形的顶点坐标分别放在向量x和y中。例5.3

axis([x1,x2,y1,y2])

定义坐标轴的刻度范围：x轴从x1到x2；y轴从y1到y2。例5.3

sin(),cos()

正弦，余弦函数例5.4

axisequal

定义x轴和y的刻度相等例5.4

eigshow(A1)

特征值和特征向量的动画程序例5.5

5.3实验内容

例5.1求解下列非齐次线性方程组，并画出二维图形来表示解的情况。

（1）











432

（2）









693

（3）











662

（4）















解：在MATLAB的M文件编辑器中，编写la12.m文件：

%绘制二元非齐次线性方程组解的图形

clear

closeall

symsx1x2%定义x1、x2为符号变量

U1=rref([1,2,5;2,-3,-4])%把方程组（1）的增广矩阵通过初等行变换

%变为最简阶梯矩阵，从而得到方程组的解

subplot(2,2,1)%准备画2×2个图形中的第一个

ezplot(\'x1+2*x2=5\')%绘制直线x1+2*x2=5

holdon%保留原来图形

ezplot(\'2*x1-3*x2=-4\')%再绘制直线2*x1-3*x2=-4

title(\'x1+2*x2=52*x1-3*x2=-4\')

%在图上标注x1+2*x2=52*x1-3*x2=-4

gridon%显示网格

U2=rref([1,3,2;3,9,6])

subplot(2,2,2)

ezplot(\'x1+3*x2=2\')

holdon

ezplot(\'3*x1+9*x2=6\')

title(\'x1+3*x2=23*x1+9*x2=6\')

gridon

U3=rref([1,-3,5;2,-6,-6])

subplot(2,2,3)

ezplot(\'x1-3*x2=5\')

holdon

ezplot(\'2*x1-6*x2=-6\')

title(\'x1-3*x2=52*x1-6*x2=-6\')

gridon

U4=rref([1,-2,3;2,1,2;1,3,5])

subplot(2,2,4)

ezplot(\'x1-2*x2=3\')

holdon

ezplot(\'2*x1+x2=2\')

ezplot(\'x1+3*x2=5\')

title(\'x1-2*x2=32*x1+x2=2x1+3*x2=5\')

gridon

holdoff

在MATLAB命令窗口中输入：

la12

运行结果为：

U1=

101

012

U2=

132

000

U3=

1-30

001

U4=

100

010

001

绘制图形如图5.1所示：

图5.1二元非齐次线性方程组解的几何意义

从运行结果可以看出，方程组（1）的解为









；方程组（2）的通解为：





























k；

方程组（3）和方程组（4）的最简行阶梯矩阵的最后一行为矛盾方程，则这两个方程组无解。

从图5.1中可以形象地看出，方程组（1）的两条直线有一个交点，故有唯一解；方程

组（2）的两条直线重合，则有无穷组解；方程组（3）的两条直线相平行，永远没有交点，

即无解；方程组（4）的三条直线不共点，则也无解。

Matlab软件针采用最小二乘法原理，针对方程组（4）这类超定方程组，给出了一种求

近似解的方法，以下给出了求解方程组（4）的程序la13.m：

%求超定方程组的最小二乘解，并用图形描述解的情况

clear

closeall

A=[1,-2;2,1;1,3];%给方程组系数矩阵A赋值

b=[3;2;5];%给常数列向量b赋值

x=Ab%求出矛盾方程组的最小二乘解

%或：x=pinv(A)*b%求出矛盾方程组的最小二乘解

symsx1x2%定义x1、x2为符号变量

ezplot(\'x1-2*x2=3\')%绘制直线x1-2*x2=3

holdon%保留原来图形

ezplot(\'2*x1+x2=2\')

ezplot(\'x1+3*x2=5\')

title(\'x1-2*x2=32*x1+x2=2x1+3*x2=5\')

%在该图上标注三个方程

plot(x(1,1),x(2,1),\'*\');%用“*”标出方程组的近似解x

gridon%显示网格

holdoff

在MATLAB命令窗口中输入：

la13

运算结果为：

1.8000

0.4000

绘制图形如图5.2所示：

图5.2二元超定方程组的近似解

针对不同种类的齐次线性方程组Ax=b，Matlab命令“x=Ab”有着不同的结果：

（1）当方程组为适定方程组时，x为方程组的唯一解；

（2）当方程组为欠定方程组时，x为方程组的一个特解；

（3）当方程组为超定方程组时，x为按照最小二乘法原理，得到的一组近似解。

例5.2求解下列线性方程组，并画出三维图形来表示解的情况。

（1）

















35.02

233

321

xxx

；（2）















321

xxx

；

（3）

















254

124

575

321

xxx

；（4）

















107

解：在MATLAB的M文件编辑器中，编写la14.m文件：

clear

closeall

A1=[1,5,-1;3,-3,-1;-2,-0.5,-1];

b1=[-1;2;-3];

U1=rref([A1,b1])

subplot(2,2,1);

ezmesh(\'x1+5*x2+1\')%绘制平面：x3=x1+5*x2+1（即x1+5*x2-x3=-1）

holdon

ezmesh(\'3*x1-3*x2-2\')

ezmesh(\'-2*x1-x2/2+3\')

title(\'方程组1\');%显示该图标题

A2=[8,1,-1;2,1,-1;-3,1,-1];

b2=[0;0;0];

U2=rref([A2,b2])

subplot(2,2,2);

ezmesh(\'8*x1+x2\');

holdon

ezmesh(\'2*x1+x2\');

ezmesh(\'-3*x1+x2\');

title(\'方程组2\');

A3=[5,-7,-1;1,4,-1;1,4,-1];

b3=[5;-12;25];

U3=rref([A3,b3])

subplot(2,2,3);

ezmesh(\'5*x1-7*x2-5\');

holdon

ezmesh(\'x1+4*x2+12\');

ezmesh(\'x1+4*x2-25\');

title(\'方程组3\');

A4=[0,5,-1;0,-7,-1;0,0,1];

b4=[8;10;15];

U4=rref([A4,b4])

subplot(2,2,4);

ezmesh(\'0*x1+5*x2-8\');

holdon

ezmesh(\'0*x1-7*x2-10\');

ezmesh(\'12\');

title(\'方程组4\');

holdoff;

在MATLAB命令窗口中输入：

la14

运行结果为：

U1=

1.0000000.9286

01.00000-0.1429

001.00001.2143

U2=

1000

01-10

0000

U3=

1.00000-0.40740

01.0000-0.14810

0001.0000

U4=

0100

0010

0001

绘制图形如图5.3所示：

图5.3三元非齐次线性方程组解的几何意义

从图5.3中可以看出，方程组（1）的解为三个平面的交点，U1的最后一列给出了该方

程组的唯一解；方程组（2）的三个平面刚好相交于同一条直线，这个齐次线性方程组的解

空间是一维的，从U2中可以得到该方程组的通解为：













k；分析矩阵U3和U4，可以得

出两个方程组都是矛盾方程，即无解，从图5.3中也可以看出方程组（3）和方程组（4）中

的三个平面没有共同的交点。

例5.3已知向量组：















u，

















v，















w，请用向量u和v来线性表示向量w，

并绘制向量图。

解：把向量

和

代入向量等式：wvxux

中，则有齐次线性方程组：







































1

。

在MATLAB命令窗口中输入：

u=[1;2];v=[2;-1];w=[8;1];

rref([u,v,w])

运行结果为：

ans=

102

013

最后一列给出了方程组的解：wvu32。

程序drawvec.m是一个绘制二维向量的函数，具体程序如下：

%函数drawvec（u）：绘制二维向量u

functiondrawvec(u)

plot([0;u(1)],[0;u(2)]);%画向量线段

holdon

theta=acos(u(1)/norm(u));%计算向量u与x轴夹角，0

if(u(2)<0)

theta=2*pi-theta;%当u向量在第三和第四象限时，theta>pi

end

fill([u(1)-0.5*cos(theta+pi/12),u(1),u(1)-0.5*cos(theta-pi/12)],[

u(2)-0.5*sin(theta+pi/12),u(2),u(2)-0.5*sin(theta-pi/12)],\'black\'

);

%给向量箭头填充

holdoff

程序la15.m绘制了向量

和向量

线性表示向量

的几何图形，其中调用的函数drawvec()，

具体程序如下：

%绘制向量的线性表示

clear

closeall

u=[1;2];v=[2;-1];w=[8;1];

u1=u*2;v1=v*3;

drawvec(u);holdon;%调用函数drawvec()，画向量u

drawvec(v);holdon;

drawvec(u1);holdon;

drawvec(v1);holdon;

drawvec(w);holdoff

axis([-210-56]),gridon;%定义坐标轴范围：-2

在MATLAB命令窗口中输入：

la15

绘制图形如图5.4所示：

图5.4向量的线性表示

从图中可以看出：向量

即为向量1u和向量1v的矢量和，满足平行四边形法则。

例5.4已知向量













x，矩阵

















A，

















A，















05.0

A，



















cossin

sincos

A，















。请分析经过线性变换xAy

后，向量

y与原向量

的几何关系4,3,2,1i，并绘制向量图。

解：在MATLAB的M文件编辑器中，编写la16.m文件

%线性变换的几何意义

x=[2;1];A1=[-1,0;0,1];A2=[1,0;0,-1];A3=[0.5,0;0,2];A4=[cos(pi/3),

sin(pi/3);-sin(pi/3),cos(pi/3)];

y1=A1*x;y2=A2*x;y3=A3*x;y4=A4*x;

subplot(2,2,1);

drawvec(x);holdon;drawvec(y1);axisequal;axis([-3,3,-1.5,2]);grid

on;

%axisequal要求横坐标与纵坐标的刻度相等

subplot(2,2,2);

drawvec(x);holdon;drawvec(y2);axisequal;axis([-3,3,-1.5,2]);grid

on;

subplot(2,2,3);

drawvec(x);holdon;drawvec(y3);axisequal;axis([-3,3,-1.5,2]);grid

on;

subplot(2,2,4);

drawvec(x);holdon;drawvec(y4);axisequal;axis([-3,3,-1.5,2]);grid

on;

在MATLAB命令窗口中输入：

la16

绘制图形如图5.5所示：

图5.5线性变换的几何意义

从图5.5中可以看出：矩阵

A对

进行线性变换的结果

y与原向量

关于y轴对称；

矩阵

A对

进行线性变换的结果

y与原向量

关于x轴对称；矩阵

A对

进行线性变换

的结果

y为把向量

的横坐标乘以0.5，把

的纵坐标乘以2得到的向量；矩阵

A对

进

行线性变换的结果

y为把向量

按顺时针方向旋转



所得到的向量。

例5.5已知矩阵：















A，

















A，















A，

















A，求

它们的特征值和特征向量，并绘制特征向量图，分析其几何意义。

解：针对矩阵

A，在MATLAB命令窗口中输入：

A1=[-1,3;2,5];

[V1,D1]=eig(A1)

eigshow(A1)%显示矩阵A1的特征值和特征向量

运行结果为：

V1=

-0.9602-0.4000

0.2794-0.9165

D1=

-1.87300

05.8730

绘制图形如图5.6所示：

图5.6运行eigshow函数的初始图

当用鼠标拖动向量

顺时针旋转时，Ax也开始旋转。向量

的轨迹为一个圆，而向量

Ax的轨迹一般情况为一个椭圆。同理，可以对其它三个矩阵进行同样的操作，绘制图形如

图5.7所示。

图5.7矩阵特征值的演示图形

函数eigshow(A)描述了向量Ax随向量

的变换关系，向量

为所有的二维单位向量，

Ax为用矩阵A对向量x进行线性变换的结果。当向量x在旋转的过程中，如果向量Ax与

向量

共线（包括同向和反向），则此时有等式xxA

成立，为一实数乘子，为正表

示两个向量同向，为负表示两个向量反向。人们把向量

与向量Ax共线的位置称为特征

位置，其中实数就称为矩阵A的特征值，而此时的

即为矩阵A的属于的特征向量。

针对矩阵

A，当向量x顺时针旋转时，向量xA

逆时针旋转，则矩阵

A必然存在实特

征值；针对矩阵

A，当向量x匀角速度顺时针旋转时，向量xA

也顺时针旋转，其角速度

一会变大，一会变小，存在四个特征位置；针对矩阵

A，当向量

匀角速度顺时针旋转时，

向量

沿一条过圆心的直线运动，此时矩阵

A有一个特征值为零；针对矩阵

A，当向量

顺时针旋转时，向量xA

也顺时针旋转，但它永远也追不上向量

，它们之间总保持着一

定的角度，则矩阵

A没有实特征值。

例5.6分析下列二次型的正定性及其对应的二次曲面

,ffxx；分析下列二次型标

准化前后所对应的二次曲线

,fxxc及

,fyyc。

（1）2

221

121

526,xxxxxxf

（2）2

221

121

243,xxxxxxf

（3）2

221

121

32,xxxxxxf

（4）2

121

3,xxxf

解：分析二次型的正定型并绘制其对应的二次曲面，在MATLAB的M文件编辑器中，编写

la17.m文件：

%绘制二次型对应的二次曲面

clear

closeall

A1=[6,-1;-1,5];A2=[3,2;2,-2];A3=[-2,0.5;0.5,-3];A4=[3,0;0,0];

lamda1=eig(A1),lamda2=eig(A2),lamda3=eig(A3),lamda4=eig(A4)

subplot(2,2,1);

ezmesh(\'6*x1^2-2*x1*x2+5*x2^2\');

%绘制二次曲面：f=6*x1^2-2*x1*x2+5*x2^2

subplot(2,2,2);

ezmesh(\'3*x1^2+4*x1*x2-2*x2^2\');

subplot(2,2,3);

ezmesh(\'-2*x1^2+x1*x2-3*x2^2\');

subplot(2,2,4);

ezmesh(\'3*x1^2\');

在MATLAB命令窗口中输入：

la17

运行结果为：

lamda1=

4.3820

6.6180

lamda2=

-2.7016

3.7016

lamda3=

-3.2071

-1.7929

lamda4=

绘制图形如图5.8所示：

图5.8二次型对应的二次曲面

为绘制二次型标准化前后所对应的二次曲线，在MATLAB的M文件编辑器中，编写

la18.m文件：

%绘制二次型标准化前后所对应的二次曲线

closeall

subplot(2,2,1);

ezplot(\'6*x1^2-2*x1*x2+5*x2^2-50\');

%绘制二次曲线：6*x1^2-2*x1*x2+5*x2^2=50

gridon;

subplot(2,2,2);

ezplot(\'4.3820*y1^2+6.6180*y2^2-50\');

%绘制对应标准型的二次曲线

gridon;

subplot(2,2,3);

ezplot(\'3*x1^2+4*x1*x2-2*x2^2-5\');

gridon;

subplot(2,2,4);

ezplot(\'3.7016*y1^2-2.7016*y2^2-5\');

gridon;

在MATLAB命令窗口中输入：

la18

绘制图形如5.9所示：

图5.9二次型标准化前后所对应的二次曲线

从图5.9中可以看出，用正交变换把二次型标准化，即是寻找一个新的直角坐标系

,yy，使得二次曲线在新坐标系中关于

y轴和

y轴对称。

5.4实验习题

1．利用最小二乘法求超定方程组



















的近似解，并绘制图形表示解的情况。

2．已知矩阵















193

632

321

A，请完成以下操作：

（1）用函数eig求A的特征值；

（2）用函数poly求A的特征多项式函数p；

（3）用ezplot函数绘制多项式函数p，观察特征多项式的零点（即特征值），从而定性地验

证eig的计算的结果。

3．已知向量组：















u，

















v，

















w，请用向量

和

来线性表示向量

，

并绘制向量图。（可以调用例5.3中的函数drawvec()）

4．下列数据表示12个点的坐标，要求：

x04610853.56.16.53.220

y.54.500

（1）把第一行的12个数放入行向量x中，把第二行的12个数放入行向量y中，请用函数

plot(x,y)绘制这12个点连成的折线；（要求：axisequal）

（2）把数据表理解为一个2×12的矩阵















X，请计算XAY

，其中

A为：















25.01

A；

















5.01

A、















A，



















cossin

sincos



















（3）请用函数plot()分别绘制矩阵

Y（4,3,2,1i）的12个点连成的折线；（要求：axisequal）

（4）用填充函数fill(x,y,’black’)替代plot函数，重新绘制图形。

5．构造一个随机二阶方阵，用函数eigshow分析其特征值和特征向量。

6．构造一个随机三阶对称矩阵，判断其正定性，绘制其对应的二次曲面；绘制其对应的的

二次曲线及标准化后的二次曲线。

应用篇

本篇包含11个线性代数的应用实验，分别介绍线性代数在数学、物理、化学、工程、

经济等各个领域中的应用，并给出了利用MATLAB来解决这些具体问题的方法。

实验6平板稳态温度的计算

6.1实验指导

本实验在简单的平板热传导模型下，利用线性方程组计算平板稳态温度的分布。当选

取节点较多时，MATLAB软件的应用就更显重要。最后利用MATLAB的绘图命令绘制平板

的温度分布，形象地表示出计算结果。表6.1给出了与本实验相关的MATLAB命令或函数。

表6.1与本实验相关的MATLAB命令或函数

命令功能说明位置

mod(i,N)

计算i除以N的余例6.2

ifmod(i,N)==0

判断i是否为N的整数倍例6.2

mesh(T)

T为一个矩阵，以T的行标为x轴值，以T的列标为y轴值，以T的

元素值为z轴值，绘制三维图形。

例6.2

6.2实验内容

例6.1在钢板热传导的研究中，常常用节点温度来描述钢板温度的分布。假设图6.1中钢板

已经达到稳态温度分布，上下、左右四个边界的温度值如图所示，而

4,321

,,TTTT

表示钢板

内部四个节点的温度。若忽略垂直于该截面方向的热交换，那么内部某节点的温度值可以近

似地等于与它相邻四个节点温度的算术平均值，如4/2010

321

TTT。请计算该钢

板的温度分布。

2020

图6.1钢板的节点分布（4个节点）

解：根据已知条件可以得到以下线性方程组：

















4/5030

4/5010

4/3020

4/2010

324

413

412

321

TTT

化简为标准的矩阵形式如下：









































4110

1401

1041

0114

在MATLAB命令窗口输入：

A=[4,-1,-1,0;-1,4,0,-1;-1,0,4,-1;0,-1,-1,4];

b=[30;50;60;80];

U=rref([A,b])

结果为：

1.000000021.2500

01.00000026.2500

001.0000028.7500

0001.000033.7500

得到方程组的解为：25.21

T℃，25.26

T℃，75.28

T℃，75.33

T℃。

例6.2在例6.1中，把钢板内部分成了2×2个节点，本例把钢板内部分为5×5个节点，如

图6.2所示。求钢板的稳态温度分布，并绘制温度分布图形。

30C

2345

78910

2122232425

...

............

30C

10C

50C50C50C50C50C

20C20C20C20C20C

图6.2钢板的节点分布（25个节点）

解：根据例6.1中的讨论，知：5×5个节点就构成了一个具有25个方程的线性方程组。在

MATLAB的M文件编辑器中编写la19.m文件：

%计算钢板的稳态温度分布

clear

closeall

N=input(\'N=\');%输入节点数，共有N×N个节点

t_u=input(\'temperatureup:\');%输入四个边界的温度值

t_d=input(\'temperaturedown:\');

t_l=input(\'temperatureleft:\');

t_r=input(\'temperatureright:\');

A=zeros(N*N);b=zeros(N*N,1);%构造N^2×N^2零矩阵A；构造N^2维零向量b

fori=1:N*N%矩阵A的主对角线元素都是4

A(i,i)=4;

end

fori=1:N*N;%给矩阵A和向量b赋值

ifi<=N%给向量b中和上边界节点对应的分量赋值

b(i)=t_u;

end

ifmod(i,N)==0%给向量b中和右边界节点对应的分量赋值

b(i)=b(i)+t_r;

end

ifmod(i,N)==1%给向量b中和左边界节点对应的分量赋值

b(i)=b(i)+t_l;

end

ifi>N*(N-1)%给向量b中和下边界节点对应的分量赋值

b(i)=b(i)+t_d;

end

ifi>N%给矩阵A中和上边界无关的节点所对应的元素赋值

A(i,i-N)=-1;

end

ifmod(i,N)~=1%给矩阵A中和左边界无关的节点所对应的元素赋值

A(i,i-1)=-1;

end

ifmod(i,N)~=0%给矩阵A中和右边界无关的节点所对应的元素赋值

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

end

ifi<=N*(N-1)%给矩阵A中和下边界无关的节点所对应的元素赋值

A(i,i+N)=-1;

end

U=rref([A,b]);%对增广矩阵进行行初等变换，化为最简行阶梯矩阵

fori=1:N%把矩阵U的最后一列按矩阵形式赋给B

forj=1:N

B(i,j)=U(N*(i-1)+j,N*N+1);

end

T(2:N+1,2:N+1)=B;%把钢板内部温度和和四周温度值放入矩阵T中

T(1,2:N+1)=t_u;

T(N+2,2:N+1)=t_d;

T(2:N+1,1)=t_l;

T(2:N+1,N+2)=t_r;

T([1,N+2],[1,N+2])=NaN;%矩阵四个角没有温度值，故把非数NaN放入

T%显示计算结果

mesh(T)%对钢板温度值绘制成曲面图形

在MATLAB的命令窗口中输入：

la19

人机对话及运算结果为：

N=5

temperatureup:20

temperaturedown:50

temperatureleft:10

temperatureright:30

NaN20.000020.000020.000020.000020.0000NaN

10.000016.565219.866721.840023.341525.312530.0000

10.000016.395821.060624.153826.212127.911130.0000

10.000017.958323.826127.500029.442630.119430.0000

10.000021.608728.787932.576933.939433.123330.0000

10.000029.687537.139540.083340.614038.434830.0000

NaN50.000050.000050.000050.000050.0000NaN

钢板的温度分布如图6.3所示。其中x、y坐标分别表示钢板横、纵方向的节点数，高

度表示节点的温度值，该三维图形形象地反映了钢板的温度分布。

图6.3钢板的温度分布

6.3实验习题

1、如图6.4所示，假设钢板已经达到稳态温度分布。分别用rTlTdTuT_,_,_,_来表示

钢板的上、下、左、右四个边界的温度值。请编写一个通用Matlab程序，求出钢板内部六

个节点的温度值，要求四个边界的温度值由用户输入。

（1）30_,20_,40_,10_rTlTdTuT，单位℃；

（2）60_,40_,80_,20_rTlTdTuT，单位℃；

（3）90_,60_,120_,30_rTlTdTuT单位℃。

123

456

图6.4钢板的节点分布（6个节点）

实验7化学方程式的配平

7.1实验指导

本实验给出了一个配平化学方程式的矩阵方法。针对复杂化学方程式，尤其是多套配

平化学方程式，MATLAB软件的应用极大地提高了计算效率。表7.1给出了与本实验相关

的MATLAB命令或函数。

表7.1与本实验相关的MATLAB命令或函数

命令功能说明位置

formatrat

定义输出格式为有理数逼近，即输出一个分数。例7.1

formatshort

定义输出格式为短格式，显示小数点后4位。缺省为formatshort

formatlong

定义输出格式为长格式，显示小数点后14位或15位。

7.2实验内容

例7.1化学反应方程式通用的配平方法有：“观察法”、“选配法”、“氧化数法”及“矩阵法”

等。请用“矩阵法”配平下列化学方程式：

（1）氧化铁和稀盐酸的化学反应：

OHFeClHClOFe

2332



（2）氯酸钾和浓盐酸的化学反应：

2223

ClOClOHKClHClKClO

解：矩阵法配平化学方程式就是用矩阵变换来求解齐次线性方程组的非零解，从而得到化学

方程式的配平系数。和其它配平方法相比，矩阵法可用于配平各种复杂的化学方程式。

（1）可以看出化学方程式

OHFeClHClOFe

2332

，共有四种元素：Fe、O、H、

Cl；四种物质：

OFe、HCl、

FeCl、OH

。设这四个分子式前的系数分别为：

x、

x和

x，则有：OHxFeClxHClxOFex

24332321

，按照方程的书写习惯，化学方

程式变为：

24332321

OHxFeClxHClxOFex

针对元素Fe，应满足方程：

00102

4321

xxxx；

针对元素O，应满足方程：01003

4321

xxxx；

针对元素H，应满足方程：02010

4321

xxxx；

针对元素Cl，应满足方程：00310

4321

xxxx；

从而构成了化学方程式的原子矩阵：















0310

2010

1003

0102F

2332

OHFeClHClOFe

元素

物质

令















x，则齐次线性方程组

















00310

02010

01003

00102

4321

xxxx

可以写成0Ax，系数矩阵A

即为原子矩阵。

在MATLAB命令窗口输入：

formatrat%用分数形式表述计算结果

A=[2,0,-1,0;3,0,0,-1;0,1,0,-2;0,1,-3,0]

U=rref(A)

20-10

300-1

010-2

01-30

100-1/3

010-2

001-2/3

0000

该矩阵对应的方程组为：













显然方程组的解不是唯一的，但化学方程式配平系数都是整数，且要求这些系数的最大公约

数为1。令3

x，则得到方程组的最小整数解为：













，配平的化学方程式为：

OHFeClHClOFe

2332

326

（2）化学方程式

2223

ClOClOHKClHClKClO，共有四种元素：K、Cl、O、

H；六种物质：

KClO、HCl、KCl、OH

、

Cl、

ClO，其原子矩阵为：























002010

201003

120111

000101

A。则在MATLAB命令窗口输入：

formatrat

A=[1,0,-1,0,0,0;1,1,-1,0,-2,-1;3,0,0,-1,0,-2;0,1,0,-2,0,0]

U=rref(A)

10-1000

11-10-2-1

300-10-2

010-200

1000-1/3-5/6

0100-2-1

0010-1/3-5/6

0001-1-1/2

此矩阵对应的方程组为：













654

653

652

651

xxx

令

0,3

xx，则有：





























，它对应的化学反应式为：

223

336ClOHKClHClKClO

令6,0

xx，则有：





























，它对应的化学方程式为：

223

63565ClOOHKClHClKClO

显然齐次线性方程组的解空间是2维的，它的通解是

kk。

当

kk，配平的化学方程式为：

2223

22242ClOClOHKClHClKClO

当

kk，配平的化学方程式为：

2223

857107ClOClOHKClHClKClO

当

kk，配平的化学方程式为：

2223

148121612ClOClOHKClHClKClO

……

这个化学方程式属于具有多套配平系数的氧化还原反应。

7.3实验习题

1、用“矩阵法”配平下列化学方程式：

（1）

OHCOOOHHC

22252



（2）

2242442224

OOHSOKMnSOSOHOHKMnO

实验8交通流量

8.1实验指导

本实验针对一个简单的城市交通模型，通过求解线性方程的方法求出车量的流量。当十

字路口增多时，计算量将增大，最后利用MATLAB软件给出了问题的答案。

8.2实验内容

例8.1某城市有如图8.1所示的交通图，每一条道路都是单行道，图中数字表示某一个时段

该路段的车流量。若针对每一个十字路口（节点），进入和离开的车辆数相等。请计算每两

个相邻十字路口间路段上的交通流量1,2,,4

xiL。

360

260

220

292320

251

260

357

图8.1单行道4节点交通流图

解：根据已知条件可以得到，四个节点的流通方程为：

节点A：x

360x

260

节点B：x

220x

292

节点C：x

320x

357

节点D：x

260x

251

将以上方程组进行整理，得：









12

=100

=72

=37







在MATLAB命令窗口输入：

A=[1,-1,0,0;0,1,-1,0;0,0,1,-1;-1,0,0,1];

b=[-100;72;37;-9];

U=rref([A,b])

计算结果为：

100-19

010-1109

001-137

00000

由于U的最后一行为全零，也就是说，四个方程中实际上只有三个独立方程。方程个

数比未知数个数少，所以该方程组为欠定方程，它存在无穷多组解。若以x

为自由变量，

方程组的解可以表示为：

109

















例8.2某城市有如图8.2所示的6节点交通图，每一条道路都是单行道，图中数字表示某一

个时段该路段的车流量。若针对每一个十字路口（节点），进入和离开的车辆数相等。请计

算每两个相邻十字路口间路段上的交通流量1,2,,7

xiL。

4088

x4x5

图8.2单行道6节点交通流图

解：本例要求计算出七个路段的交通流量，而交通图中只有六个节点（十字路口），即有六

个约束条件。显然没有唯一的解答。根据图中A、B、C、D、E及F六个节点四周的流量，

列出如下方程组：

Axb，其中

1000010

1100001

0110000

0011000

0001101

0000110



























，











，



















在MATLAB命令窗口中输入：

A=[-1,0,0,0,0,-1,0;1,-1,0,0,0,0,1;0,1,-1,0,0,0,0;0,0,1,-1,0,0,0;0

,0,0,1,-1,0,-1;0,0,0,0,1,1,0];

b=[-70;40;-58;69;-25;44];

U=rref([A,b])

计算结果为：

100001070

010001-130

001001-188

000101-119

000011044

00000000

由于U的最后一行为全零，即，六个方程中实际上只有五个独立方程。而方程组未知数个

数为七，则它存两个自由变化的量。若设

15;20xx，则方程组的解为：

267

367

467

7055

3035

8893

1924

4429

xxx

































8.3实验习题

1．某城市有如图8.3所示的9节点交通图，每一条道路都是单行道，图中数字表示某一个

时段该路段的车流量。若针对每一个十字路口，进入和离开的车辆数相等。请计算每两个相

邻十字路口间路段上的交通流量1,2,,12

xiL。

若已知AB段和FO段在修路，即：

0x，

0x，又已知

300x，

660x。

求此时各个路段的交流流量。

120

350

300100

266

500

x5x6

x9x12

CDE

x10

x11

100

23440050

100

图8.3单行道9节点交通流图

实验9情报检索问题

9.1实验指导

现代情报检索技术是在矩阵理论的基础上发展起来的。情报中心的数据库中存放着大量

文件，读者希望从中搜索到与自己特定关键词相匹配的文件。文件的类型可以是书籍、期刊

杂志、论文、研究报告、因特网上的网页、胶片等等。

假如数据库中包括了n个文件，而搜索所用的关键词有m个。关键词按字母顺序排列，

我们就可以把数据库用一个mn的矩阵A来表示。矩阵A的第j列表示第j个文件，矩阵A

的第i行表示第i个关键词，矩阵A的元素

表示第i个关键词出现在第j个文件中的相对

频率。用于检索的关键词清单用m维列向量

表示。假如检索策略中含有第i、j、k三个关

键词，则让搜索向量

的第i、j、k个元素为1，其它元素为0。当确定了数据库矩阵A和

搜索向量

后，进行线性变换运算：TyAx，其中向量y即为检索结果。

9.2实验内容

例9.1假如某数据库包含以下8种图书：B1：线性代数，B2：初等线性代数，B3：初等线性

代数及其应用，B4：线性代数及其应用，B5：矩阵代数及其应用，B6：矩阵理论，B7：线性

代数及MATLAB入门，B8：基于MATLAB的线性代数及其应用。而检索的7个关键词按拼音字

母次序排列为：“初等，代数，矩阵，理论，MATLAB，线性，应用”。读者1的检索策略为：

“代数，MATLAB”；读者2的检索策略是：“代数，MATLAB，应用”。请用矩阵运算来为这两

位读者检索图书。

解：因为这些关键词在书名中出现最多只有1次，所以其相对频率数不是0就是1。当第i

个关键词出现在第j本书名上时，矩阵A的元素

a就等于1，否则就等于0。这样数据库矩

阵就可以用表9.1来表示：

表9.1情报检索矩阵图

关键词

书

B1B2B3B4B5B6B7B8

初等01100000

代数11111011

矩阵00001100

理论00000100

MATLAB00000011

线性11110011

应用00111001

读者1输入的关键词是：“代数，MATLAB”；读者2输入的关键词是：“代数，MATLAB，

应用”，则数据库矩阵A、搜索向量

x及

x分别为：

01100000

11111011

00001100

00000100

00000011

11110011

00111001











，











，











搜索结果可以通过以下线性变换求得：

TyAx

，

TyAx，在MATLAB命令窗口输

入：

A=[0,1,1,0,0,0,0,0;1,1,1,1,1,0,1,1;0,0,0,0,1,1,0,0;0,0,0,0,0,1,0,

0;0,0,0,0,0,0,1,1;1,1,1,1,0,0,1,1;0,0,1,1,1,0,0,1];

x1=[0;1;0;0;1;0;0];

x2=[0;1;0;0;1;0;1];

y1=A\'*x1

y2=A\'*x2

计算结果为：

y1=

y2=

向量

,yy的各个分量分别表示八种图书与搜索向量匹配的程度。y1的第7个和第8个

分量都为2，说明后两本书B7和B8必然包含读者1的检索策略中的两个关键词，这两本书

就被认为具有最高的匹配度；同理，y2的第8个分量为3，说明图书B8包含了读者2的检

索策略中的三个关键词，这本书就被认为具有最高的匹配度。

9.3实验习题

1、假如一个数据库包含以下10种图书：B1：高等代数，B2:线性代数，B3：工程线性代数，

B4：初等线性代数，B5：线性代数及其应用，B6：MATLAB在数值线性代数中应用，B7：矩

阵代数及其应用，B8：矩阵理论，B9：线性代数及MATLAB入门，B10：基于MATLAB的线性

代数及其应用。而检索的6个关键词按拼音字母次序排列为：“代数，工程，矩阵，MATLAB，

数值，应用”。读者1的检索策略为：“代数，MATLAB”；读者2的检索策略是：“代数，应用”。

请用矩阵运算来为这两位读者检索图书。

实验10图与矩阵

10.1实验指导

图论是应用数学的一个重要分支，它广泛地应用于各个领域。本实验介绍了一个图论的

基本应用，并用MATLAB给出了具体解答。

图10.1描述了四个城市之间的航空航线图，其中1、2、3、4表示四个城市；带箭头线

段表示两个城市之间的航线。从图中可以看出：城市1到城市2有航班，而城市2到城市1

没有航班；城市1和城市4之间来回都有航班；城市3和城市4之间没有航班。为了描述这

四个城市航线的邻接关系，定义邻接矩阵A为：

图10.1航空航线图（四城市）

0111

0011

0000

1100











其中，第i行描述从城市i出发，可以到达各个城市的情况，若能到达第j个城市，则

A(i,j)=1，否则A(i,j)=0，规定A(i,i)=0,(其中i=1,2,3,4)。如第2行表示：从城市2出发可以到

达城市3和城市4，不能到达城市1和城市2。从数学上可以证明：矩阵2A表示一个人连续

坐两次航班可以到达的城市，矩阵3A表示连续坐三次航班可以到达的城市，如：

1222

0122

0000

2211











分析矩阵3A的第二行，可以得出：某人从城市2出发，连续坐三次航班可以到达城市2、

城市3和城市4，不能到达城市1，而到达城市3和城市4的方法各有两种。

表10.1给出了与本实验相关的MATLAB命令或函数。

表10.1与本实验相关的MATLAB命令或函数

命令功能说明位置

max(B)

计算矩阵B每一列的最大值，计算结果是一个行向量。例10.1

max(max(B))

计算矩阵B的列最大值构成向量的最大值，结果即为矩阵B的最大值。例10.1

[i,j]=find(B==m)

寻找矩阵B中值为m的元素的位置，i和j为同维列向量，它们分别

存放行值和列值。

例10.1

10.2实验内容

例10.1图10.2描述了六个城市之间的航空航线图，其中1、2、……、6表示六个城市，带

箭头线段表示两个城市之间的航线。用MATLAB软件完成以下操作：

（1）构造该图的邻接矩阵A；

（2）若某人连续乘坐五次航班，那么他从哪一个城市出发到达哪一个城市的方法最多？

（3）若某人可以乘坐一次、二次、三次或四次航班，那么他从哪一个城市出发总是不能达

到哪一个城市？

图10.2航空航线图（六城市）

解：（1）根据前面的讨论及图10.2，构造邻接矩阵A；

（2）计算矩阵5A，找出该矩阵的最大元素，并确定它所在的位置；

（3）计算矩阵234AAAA，找出该矩阵中零元素的位置。

在MATLAB软件的M编辑器中编写la20.m文件：

%图与矩阵

clear

A=[0,1,0,0,0,1;0,0,1,1,0,0;0,0,0,1,1,0;0,1,0,0,0,0;1,0,1,0,0,0;0,

1,0,0,1,0];%构造邻接矩阵

B=A^5;

C=A+A^2+A^3+A^4;

disp(\'邻接矩阵A为：\');

disp(A);

disp(\'矩阵A^5为：\');

disp(B);

m=max(max(B));%计算矩阵B的最大值

[m_i,m_j]=find(B==m);%寻找矩阵B中元素等于m的位置

fprintf(\'矩阵A^5最大值%d的位置在：n\',m);

disp([m_i,m_j]);

disp(\'矩阵A＋A^2＋A^3＋A^4为：\');

disp(C);

[z_i,z_j]=find(C==0);%寻找矩阵C中零元素的位置

disp(\'矩阵A＋A^2＋A^3＋A^4零元素的位置在：\');

disp([z_i,z_j]);

在MATLAB命令窗口中输入：

la20

计算结果为：

邻接矩阵A为：

010001

001100

000110

010000

101000

010010

矩阵A^5为：

255531

244320

235521

021321

264541

144742

矩阵A^5最大值7的位置在：

矩阵A＋A^2＋A^3＋A^4为：

265642

144631

254541

133310

356642

366541

矩阵A＋A^2＋A^3＋A^4零元素的位置在：

从计算结果中可以看出，矩阵A^5最大值出现在矩阵的第六行第四列，说明：这个人

如果从城市6出发连续乘坐五次航班后到达城市4，他可以选择的乘机路线最多，共有7种

不同的方法。

矩阵A＋A^2＋A^3＋A^4的零元素出现在第四行第六列，说明：这个人如果从城市4

出发他乘坐一次、二次、三次或四次航班，都无法到达城市6。

10.3实验习题

1、5个小朋友玩传球游戏。游戏规则：任意两个人之间都可以相互传球，但自己不能给自

己传。请用MATLAB完成以下操作：

（1）把五个小朋友看成五个节点，构造这五个节点的邻接矩阵A；

（2）假设从第一个小朋友开始传球，经过四次传球后，球又传回到第一个小朋友手里。问

共有多少种不同的传法。

（3）假设从第一个小朋友开始传球，经过一次，或者二次，或者三次传球，球传给了第二

个小朋友。问共有多少种传法。

实验11利用行列式计算面积和体积

11.1实验指导

本实验利用二阶行列式和三阶行列式的几何意义分别计算平行四边形的面积及平行六

面体的体积。并用MATLAB软件来实现具体运算。

二阶行列式和三阶行列式的几何意义如下：由向量

,uab和

,vab所构成的

平行四边形的面积为行列式11

的绝对值，如图11.1（1）所示；由向量

111

,,uabc，



222

,,vabc和

333

,,wabc所构成的平行六面体的体积为行列式

111

222

333

abc

的绝对

值，如图11.1（2）所示；

（1）向量构成的平行四边形（2）向量构成的平行六面体

图11.1行列式的几何意义

表11.1给出了与本实验相关的MATLAB命令或函数。

表11.1与本实验相关的MATLAB命令或函数

命令功能说明位置

function

s=cal_are3(a,b,c)

cal_are3为函数名，参数为三个2维向量，s为函数返回值例11.1

&逻辑运算符号：逻辑与，与C语言&&功能类似

例11.1

逻辑运算符号：逻辑或，与C语言||功能类似例11.1

return

结束当前程序，返回调用处。例11.1

A=[ab;ac]

把行向量ab和行向量ac分别放到矩阵A的第一行和第二行例11.1

A=[ab;ac]

把列向量ab和列向量ac分别放到矩阵A的第一列和第二列例11.1

abs

针对复数取幅值函数；针对实数取绝对值。例11.1

11.2实验内容

例11.1根据二阶行列式的几何意义，利用MATLAB软件编写通用函数，计算三角形ABC

的面积。其中三角形三个顶点坐标以参数形式输入。

（1）已知三角形ABC三个顶点的坐标为：（1，5）、（4，3）、（2，-1），利用通用函数计算

该三角形的面积；

（2）已知凸五边形ABCDE五个顶点的坐标分别为：（-3，4）、（1，5）、（5，1），（-1，-1）、

（-6，1），利用通用函数计算该五边形的面积。

（3）在平面坐标系中画出以上三角形和五边形。

解：设三角形ABC顶点坐标分别为：

,ab、

,ab和

,ab，则向量

ABOBOA

uuuruuuruuur



2121

,aabb，向量ACOCOA

uuuruuuruuur



3131

,aabb，而三角形ABC

的面积就等于向量AB

uuur

和向量AC

uuur

所构成平行四边形面积的一半，如图11.2所示。根据行

列式的几何意义，可以得到计算三角形ABC面积的公式如下：

2121

3131

0.5

aabb

Sabs

aabb











，abs为取绝对值。

图11.2向量构成的三角形

在MATLAB软件的M编辑器中编写计算三角形面积的通用函数cal_area3.m：

%函数cal_area3()利用行列式计算三角形面积

functions=cal_are3(a,b,c)

ifsize(a)~=[1,2]&size(a)~=[2,1]|size(a)~=size(b)|size(a)~=size(c)

%确定参数格式是否正确

%a,b,c应为同形的2维行向量或列向量

disp(\'参数格式错误\');

return;%结束当前程序，返回调用处

end

ab=b-a;%计算向量AB

ac=c-a;%计算向量AC

ifsize(ab)==[1,2]%判读向量AB是否为行向量

A=[ab;ac];%构造矩阵A

else

A=[ab,ac];

end

s=abs(det(A))/2;%根据公式计算三角形面积

（1）在MATLAB的M文件编辑器中编写程序la21.m：

%计算三角形的面积

a=input(\'请输入A点坐标：\');%输入三角形顶点坐标

b=input(\'请输入B点坐标：\');

c=input(\'请输入C点坐标：\');

S=cal_area3(a,b,c);%调用函数cal_area3

fprintf(\'三角形的面积为：%d\',S);%输出计算结果计算结果为：

在MATLAB的命令窗口输入：

la21

人机对话结果为：

请输入A点坐标：[1,5]

请输入B点坐标：[4,3]

请输入C点坐标：[2,-1]

三角形的面积为：8

（2）在MATLAB软件的M文件编辑器中，编写计算凸多边形面积的通用函数cal_arean.m：

%函数cal_arean()利用行列式计算凸多边形面积

functions=cal_arean(A)

[m,n]=size(A);

if(m~=2&n~=2)|(m<3&n<3)%确定参数格式是否正确

disp(\'参数格式错误\');%A应为2×n或m×2矩阵

return;%且n和m不能同时小于3

end

ifm==2%如果输入矩阵为2×n

A=A\';%把矩阵转置

[m,n]=A;%重新取得矩阵行与列

end

s=0;%给结果变量赋初值零

%计算凸m边形的面积

fori=1:m-2%累加m-2个三角形面积

s=s+cal_area3(A(1,:),A(1+i,:),A(2+i,:));

%调用计算三角形面积函数cal_area3()

end

在MATLAB的M文件编辑器中编写程序la22.m：

%计算凸多边形的面积

A=input(\'请输入凸多边形顶点的坐标：\');%输入凸多边形顶点坐标

S=cal_arean(A);%调用函数cal_arean

fprintf(\'凸多边形的面积为：%d\',S);%输出计算结果计算结果为：

在MATLAB的命令窗口输入：

la22

人机对话结果为：

请输入凸多边形顶点的坐标：[-3,4;1,5;5,1;-1,-1;-6,1]

凸多边形的面积为：3.750000e+001

矩阵A的五行分别表示凸五边形ABCDE五个顶点的坐标，该多边形的面积为37.5。

（3）在MATLAB的M文件编辑器中编写绘制三角形和五边形的程序la23.m：

%画三角形

closeall

A=[1,5;4,3;2,-1;1,5];%矩阵A的最后一行和第一行相同，目的是画出闭合图形

subplot(1,2,1);

plot(A(:,1),A(:,2));%以矩阵A的第一列为横坐标，以矩阵A的第二列为纵坐标

axisequal;

axissquare;

gridon;

%画五边形

A=[-3,4;1,5;5,1;-1,-1;-6,1;-3,4];

subplot(1,2,2);

plot(A(:,1),A(:,2));

axisequal;

axissquare;

gridon;

在MATLAB命令窗口中输入：

la23

绘制图形如图11.1所示。

图11.1凸多边形的绘制图

11.3实验习题

1、如图11.3所示，求向量1,2,3u、3,1,0v、0,5,1w所构成的四面体体积。

图11.3三个3维向量所构成的四面体

2、计算凸十边形的面积。其中十个顶点的坐标为：（0，8）、（3，7）、（4，5）、（4，0）、（3，

-4）、（1，-5）、（-3，-5）、（-6，0）、（-5，6）、（-2，8）。

实验12中成药的配制问题

12.1实验指导

本实验给出了一个中成药的配制问题，其中涉及到了向量组的线性相关性、向量组的最

大无关组、向量的线性表示及向量空间等线性代数知识。并利用MATLAB软件完成具体计

算。

12.2实验内容

例12.1某中药厂用9种中草药（A、B、……、I），根据不同的比例制成了7种特效药。表

12.1给出了每种特效药每包所需各种成分的质量（单位：克）。

表12.1特效药的成分含量（单位：克）

1号成药2号成药3号成药4号成药5号成药6号成药7号成药

A10

B12

C531105140

D79255154735

E012255336

F255355355550

G94172523925

H651610103510

I821200620

（1）某医院要购买这7种特效药，但药厂的第3号和第6号特效药已经卖完，请问能否用

其它特效药配制出这两种脱销的药品。

（2）现在该医院想用这9中种草药配制三种新的特效药，表12.2给出新药所需的成分质量

（单位：克）。请问如何配制。

表12.2三种新药的成分含量（单位：克）

1号新药2号新药3号新药

A4016288

B6214167

C14278

D4410251

E53607

F5015580

G7111838

H416821

I145230

解：（1）把每一种特效药都看成一个9维列向量，分析向量组

127

,,,uuuL的线性相关性。

若向量组线性无关，则无法配制脱销的特效药；若向量组相关，并且能找到不含

u和

u的

一个最大无关组，则可以利用现有的特效药来配制第3号和第6号药品。

在MATLAB的命令窗口中输入：

u1=[10;12;5;7;0;25;9;6;8];

u2=[2;0;3;9;1;5;4;5;2];

u3=[14;12;11;25;2;35;17;16;12];

u4=[12;25;0;5;25;5;25;10;0];

u5=[20;35;5;15;5;35;2;10;0];

u6=[38;60;14;47;33;55;39;35;6];

u7=[100;55;0;35;6;50;25;10;20];

U=[u1,u2,u3,u4,u5,u6,u7]

[U0,r]=rref(U)

计算结果为：

531105140

79255154735

012255336

255355355550

94172523925

651610103510

821200620

U0=

1010000

0120030

0001010

0000110

0000001

0000000

12457

从最简行阶梯矩阵U0中可以看出：特效药向量组是线性相关的，它的秩是5，它的一

个最大无关组是：

12457

,,,,uuuuu，则可以用现有特效药来配制出3号和6号药品：

312

2uuu，

6245

3uuuu。

（2）三种新药分别用向量

,vv和

v表示，把特效药向量组和新药向量组放入同一个矩阵中：



1234567123

,,,,,,,,,Uuuuuuuuvvv，经过初等行变化，矩阵U变为最简行阶梯矩阵U0，

从U0中的后3列就可以获得答案。在MATLAB命令窗口中输入：

u1=[10;12;5;7;0;25;9;6;8];

u2=[2;0;3;9;1;5;4;5;2];

u3=[14;12;11;25;2;35;17;16;12];

u4=[12;25;0;5;25;5;25;10;0];

u5=[20;35;5;15;5;35;2;10;0];

u6=[38;60;14;47;33;55;39;35;6];

u7=[100;55;0;35;6;50;25;10;20];

v1=[40;62;14;44;53;50;71;41;14];

v2=[162;141;27;102;60;155;118;68;52];

v3=[88;67;8;51;7;80;38;21;30];

U=[u1,u2,u3,u4,u5,u6,u7,v1,v2,v3];

[U0,r]=rref(U)

计算结果为：

U0=

1010000130

1245710

从最简行阶梯矩阵U0的后3列可以看出：

1124

32vuuu，

21247

342vuuuu,

而

v不能由前七种特效药线性表示，即无法配制出第三种新药。

12.3实验习题

1、一个混凝土生产企业可以生产出三种不同型号的混凝土，它们的具体配方比例如表11.3

所示。

表11.3混凝土的配方

型号1混凝土型号2混凝土型号3混凝土

水101010

水泥222618

砂323129

石子536450

灰058

（1）分析这三种混凝土是否可以用其中的两种来配出第三种？

（2）现在有甲、乙两个用户要求混凝土中含水、水泥、砂、石子及灰的比例分别为：24，

52，73，133，12和36，75，100，185，20。那么，能否用这三种型号混凝土配出满足甲和

乙要求的混凝土？如果需要这两种混凝土各500吨，问三种混凝土各需要多少？

实验13人口迁徙问题

13.1实验指导

本实验从人口迁移问题出发，讨论了方程组、矩阵乘法、特制值特征向量、及矩阵对

角化问题。并利用MATLAB软件对线性代数理论进行了验证。表13.1给出了与本实验相关

的MATLAB命令或函数。

表13.1与本实验相关的MATLAB命令或函数

命令功能说明位置

symsn

定义符号变量n例13.1

lamda.^n

群运算，对矩阵lamda中所有元素进行幂运算例13.1

13.2实验内容

例13.1假设一个城市的总人口数是固定不变，但人口的分布情况变化如下：每年都有5％的

市区居民搬到郊区；而有15％的郊区居民搬到市区。若开始有700000人口居住在市区，

300000人口居住在郊区。请利用MATLAB软件分析：

（1）10年后市区和郊区的人口各是多少？

（2）30年后、50年后市区和郊区的人口各是多少？

（3）分析（2）中数据相似的原因。

解：设第n年市区人数和郊区人数分别为：

x、

y，则第n＋1年的市区和郊区人数为：

0.950.15

0.050.85

nnn

xxy

yxy













。用矩阵表示为：1

0.950.15

0.050.85

















，进一步写为：

AXX

1

，矩阵

0.950.15

0.050.85









，向量

X描述第n年该城市的市区和郊区人数。

（1）根据

AXX

1

，及开始市区和郊区的人口数

700000

300000









，可以得到10年后市

区和郊区的人口分布：210

10980

XAXAXAXL，在MATLAB命令窗口中输入：

A=[0.95,0.15;0.05,0.85];

X0=[700000;300000];

X10=A^10*X0

计算结果为：

X10=

1.0e+005*

7.4463

2.5537

可以知道，10年后市区和郊区人口数约为：744630和255370。

（2）同理可以得到30年后和50年后的人口分布：30

300

XAX和50

500

XAX。在MATLAB

命令窗口中输入：

A=[0.95,0.15;0.05,0.85];

X0=[700000;300000];

X30=A^30*X0

X50=A^50*X0

计算结果为：

X30=

1.0e+005*

7.4994

2.5006

X50=

1.0e+005*

7.5000

2.5000

从计算结果可以看出30年后和50年后的市区和郊区人口分布非常接近。

（3）若矩阵A可以对角化，则可以利用矩阵对角化的方法来计算矩阵nA。

设矩阵A的特征值构成的对角矩阵为，则总存在可逆矩阵P，使得：1APP，

则1nnAPP，其中矩阵P为A的特征值所对应的特征列向量所构成的矩阵。那么有：

XAXPPX。可以利用MATLAB的符号运算来计算n年后该城市的人口分布

情况

X。在MATLAB的M文件编辑器中编写程序la24.m：

%分析n年后城市人口分布

A=[0.95,0.15;0.05,0.85];

X0=[700000;300000];

[P,lamda]=eig(A);

symsn%定义符号变量n

Xn=P*lamda.^n*inv(P)*X0%.^n群运算，对矩阵lamda中所有元素进行幂运算

在MATLAB命令窗口输入：

la24

计算结果为：

Xn=

[750000-50000*(4/5)^n]

[250000+50000*(4/5)^n]

从计算结果可以看出：当n越大，(4/5)^n就越近于零，则

X就趋近于

750000

250000







，这就是

为什么

X与

X接近的原因。

13.3实验习题

1、李博士培养了一罐细菌，在这个罐子里存放着A、B、C三类不同种类的细菌，最开始A、

B、C三种细菌分别有810、2×810、3×810个。但这些细菌每天都要发生类型转化，转化

情况如下：A类细菌一天后有5％的变为B类细菌、15％的变为C细菌；B类细菌一天后有

30％的变为A类细菌、10％的变为C类细菌；C类细菌一天后有30％的变为A类细菌、20％

的变为B类细菌。请利用MATLAB软件分析：

（1）一周后李博士的A、B、C类细菌各有多少个？

（2）两周后和三周后李博士的A、B、C类细菌各有多少个？

（3）分析在若干周后，李博士的各种细菌的个数几乎不发生变化的原因。

实验14多项式插值与曲线拟合

14.1实验指导

本实验介绍了多项式插值和曲线拟合的概念。利用求解适定线性方程组实现函数插值；

利用求解超定线性方程组来实现曲线拟合。并用MATLAB软件绘制出插值和拟合的图形。

表14.1给出了与本实验相关的MATLAB命令或函数。

表14.1与本实验相关的MATLAB命令或函数

命令功能说明位置

x.^4

群运算，对向量x中所有元素进行幂运算。例14.1

xi=linspace(-1,9,100)

构造具有100个元素的一维数组，其值从-1到9均匀分布。例14.1

plot(x,y,’o’)

在图中用符号“O”来标注位置x,y。例14.1

p=polyfit(x,y,2)

对数据x,y进行2次多项式拟合，p是拟合结果，按从高次

幂到低次幂排列的系数向量。

例14.1

14.2实验内容

例14.1表14.2给出了平面坐标系中五个点的坐标。

表14.2五点数据表

x01234

y-270210-75

（1）请过这五个点作一个四次多项式函数234

401234

()pxaaxaxaxax，并求当

5x时的函数值

5p。用MATLAB绘制多项式函数

()px曲线、已知点及插值点（5，



5p）。

（2）请根据这五个点，拟合一个二次多项式函数2

2012

()pxaaxax，并用MATLAB

绘制多项式函数

()px曲线及已知的五个点。

解：（1）根据已知条件，把五个点的坐标值分别代入四次多项式函数，可以得到如下线性方

程组：

234

01234

234

01234

234

01234

234

01234

234

01234

000027

11110

222221

33330

444475

aaaaa

























对应矩阵等式为：Aay，其中

234

10000

11111

12222

13333

14444











，











，















系数矩阵A的行列式为范德蒙行列式，且五个坐标点的横坐标各不相同，则该行列式不等

于零，所以方程组有唯一解。

（2）根据已知条件，把五个点的坐标值分别代入二次多项式函数，可以得到如下线性方程

组：

012

0027

110

2221

330

4475

aaa

























对应矩阵等式为：Aay，其中

100

111

122

133

144











，











，















该方程组有三个未知数，但有五个方程，进一步分析可以得到该方程组无解，即不存在一个

二次多项式曲线刚好能过已知的五个点。MATLAB软件提供了一个利用最小二乘法解决超

定方程组近似解的方法。即可以找到一条二次曲线来近似地描述已知5点的变化情况。

在MATLAB软件M文件编辑器中编写程序la25.m：

%多项式插值和函数逼近

clear

closeall

x=[0;1;2;3;4];%输入已知点坐标

y=[-27;0;21;0;-75];

A=[x.^0,x.^1,x.^2,x.^3,x.^4];%构造范德蒙矩阵

a=Ay;%得到适定方程组的唯一解a，即确定了多项式函数

%或p=polyfit(x,y,4)%p是按从高次幂到低次幂排列的系数向量；

disp(\'四次多项式系数为：\')

disp(a);

xi=linspace(-1,9.5,100);%构造数组xi，从-1到9.5均匀取100个值

yi=a(1)+a(2)*xi+a(3)*xi.^2+a(4)*xi.^3+a(5)*xi.^4;

%计算对应xi的多项式函数值yi

x0=5;

y0=a(1)+a(2)*x0+a(3)*x0^2+a(4)*x0^3+a(5)*x0^4;

%计算插值点函数值

disp(\'四次多项式函数插值点p(5)=\');

disp(y0);

subplot(1,2,1);

plot(xi,yi,x,y,\'o\',x0,y0,\'*\');%绘制四次多项式函数、已知五点及插值点

axissquare;%使坐标轴为正方形

axis([-19-400100])%确定x轴和y轴范围

gridon;%显示网格

A=[x.^0,x.^1,x.^2];

a=Ay;%根据最小二乘法得到超定方程组的近似解a

%或p=polyfit(x,y,2)%p是按从高次幂到低次幂排列的系数向量；

disp(\'二次多项式系数为：\')

disp(a);

xi=linspace(-1,5,100);%构造数组xi，从-1到5均匀取100个值

yi=a(1)+a(2)*xi+a(3)*xi.^2;%计算对应xi的多项式函数值yi

subplot(1,2,2);

plot(xi,yi,x,y,\'o\');%绘制二次多项式函数及已知五点

axissquare;

axis([-15-15050])

gridon;

在MATLAB命令窗口输入：

la25

计算结果为：

四次多项式系数为：

-27

-12

四次多项式函数插值点p(5)=

-192

二次多项式系数为：

-32.1429

60.6857

-17.5714

图14.1给出了MATLAB绘制的图形。从图中可以形象地看出插值和拟合的区别。

图14.1插值和拟合的示意图

14.3实验习题

1、表14.3给出了平面坐标系中六个点的坐标。

表14.3六点数据表

x012345

y-7540

（1）请过这六个点作一个五次多项式函数2345

5012345

()pxaaxaxaxaxax

，并

求当1x时的函数值

1p。用MATLAB绘制多项式函数

()px曲线、已知点及插值

点（－1，

1p）。

（2）请根据这六个点，拟合一个三次多项式函数23

30123

()pxaaxaxax。并用

MATLAB绘制多项式函数

()px曲线及已知点。

实验15刚体的平面运动

15.1实验指导

本实验介绍利用矩阵乘法来实现刚体的平面运动。并用MATLAB对一个具体实例进行

分析，最后绘制刚体运动前后的图形。

用平面坐标系中的一个闭合图形来描述刚体，用一个矩阵X来表示它。X的一列表示

刚体一个顶点的坐标。为了使图形闭合，X的最后一列和第一列相同；为了实现刚体的平移

运算，给矩阵X添加元素值都为1的一行，使矩阵X的形状为3×n。

若有矩阵：

001











，

cossin0

sincos0

001

Rtt













，且：

YMX，

YRX

可以证明，矩阵

Y是刚体X沿x轴正方向平移

c，沿y轴正方向平移

c后的结果；矩阵

是刚体X以坐标原点为中心逆时针转动t弧度的结果。

15.2实验内容

例15.1用下列数据表示大写字母A。对图形A进行以下平面运动，并绘制移动前后的图形。

x04610853.56.16.53.220

y.54.500

（1）向上移动15，向左移动30；

（2）逆时针转动

；

（3）先逆时针转动

pi，然后向上移动30，向右移动20。

解：（1）构造刚体矩阵

04610853.56.16.53.220

.54.500











及平移

矩阵

1030

0120

001













。

（2）构造转动矩阵

cossin0

sincos0

001















。

（3）构造转动矩阵

cossin0

sincos0

001















及平移矩阵

1020

0130

001











。

在MATLAB的M文件编辑器中编写程序la26.m

%刚体的平面运动

closeall

X=[0,4,6,10,8,5,3.5,6.1,6.5,3.2,2,0;0,14,14,0,0,11,6,6,4.5,4.5,0,0;ones(1,12)];

%构造刚体矩阵

M=[1,0,-30;0,1,15;0,0,1];%构造平移矩阵

Y1=M*X;%计算平移结果

plot(X(1,:),X(2,:));%绘制原来刚体

holdon

axisequal

fill(Y1(1,:),Y1(2,:),\'red\');%绘制平移后刚体

R=[cos(pi/3),-sin(pi/3),0;sin(pi/3),cos(pi/3),0;0,0,1];

%构造转动矩阵

Y2=R*X;

fill(Y2(1,:),Y2(2,:),\'blue\');%绘制转动后刚体

M=[1,0,20;0,1,30;0,0,1];

R=[cos(3*pi/4),-sin(3*pi/4),0;sin(3*pi/4),cos(3*pi/4),0;0,0,1];

Y3=M*R*X;

fill(Y3(1,:),Y3(2,:),\'black\');%绘制转动及平移后刚体

gridon

holdoff

在MATLAB命令窗口中输入：

la26

绘制图形如图15.1所示。

图15.1刚体的平面运动

15.3实验习题

1、用MATLAB软件实现以下操作：

（1）构造一个直角三角形刚体矩阵X；

（2）先对刚体逆时针转动



，然后再向下移动20，向右移动20；

（3）先对刚体向下移动20，向右移动20，然后对刚体逆时针转动



。

实验16投入产出问题

16.1实验指导

本实验根据经济活动各部门投入与产出平衡关系的思想，利用矩阵方程计算该平衡关

系下的各个要素的分布。当生产消耗关系矩阵较大时，MATLAB软件的应用就更显重要。

最后利用MATLAB得到计算结果。表16.1给出了与本实验相关的MATLAB命令或函数。

表16.1与本实验相关的MATLAB命令或函数

命令功能说明位置

sum(A)

对矩阵A的列求和例16.1

16.2实验内容

一般地，价值型投入产出如表1.2，其中

11121

21222

nnnn

aaa











MMLM

为中间产品矩阵，

,,,

DdddL为固定资产折旧向量，

,,,

ZzzzL为新创造价

值向量，

,,,

YyyyL为最终产品向量．

表1.2价值型投入产出表

产出

投入

中间产品最终产品总

产

品

量

消

费

积

累

出

口

小

计

资

料

补

偿

价

值

n

固定资产折旧

新

创

造

价

值

劳动报酬

纯收入

小计

总投入

注：

a表示第j部门在生产过程中消耗第i部门的产品数量

例16.1在问题中若已知某地区在某一生产周期内各部门之间的生产消耗关系矩阵A，固定

资产折旧向量D及新创造价值向量Z分别为

18376543

26358473

69242135

47825137











，15201017D，25Z

求出总投入向量

1234

Xxxxx、最终产品向量

1234

Yyyyy及直接消耗系数

矩阵

Bb（其中ij

）注：

b表示第j部门消耗第i部门的产品价值在第j部门

的总产品价值中所占比率．

解：根据各部门的总投入1,2,3,4

xj等于生产资料补偿价值与新创造价值之和，得到

总投入向量sum(A)+D+ZX；最终产品向量(sum(A))YX



；直接消耗系数矩阵

ij

Bb（其中ij

）；

在MATLAB的M文件编辑器中编写程序la27.m：

A=[18376543;26358473;69242135;47825137];

D=[15201017];

Z=[25];

X=sum(A)+D+Z;

Y=(X-sum(A\'))\';

B=zeros(size(A));

fori=1:size(A)

forj=1:size(A)

B(i,j)=A(i,j)/X(j);

end

在Matlab命令窗口输入：

la27

计算结果为：

425578893589

262

360

744

372

0.04240.06400.07280.0730

0.06120.06060.09410.1239

0.16240.04150.02350.0594

0.11060.14190.05710.0628

例16.2.若已知四个部门的直接消耗系数矩阵B与总产值X及固定资产折旧D分别为

00.150.550

0.250.050.10.25

0.1500.050.35

0.10.150.150.1











，

360

240

180

300











，5151020D

求出各部门新造价值Z，各部门最终产品Y，各部门中间产品

a．

解：由题意，生产消耗关系矩阵为

()

Aa

（其中

ijijj

abx

）；由总投入等于总产值，新创

造价值向量sum(A)DZX



；最终产品向量(sum(A))YX



；

在MATLAB的M文件编辑器中编写程序la28.m：

B=[0,0.15,0.55,0;0.25,0.05,0.1,0.25;0.15,0,0.05,0.35;0.1,0.15,0.1

5,0.1];

X=[360;240;180;300];

D=[5,15,10,20];

A=zeros(size(B));

fori=1:size(B)

forj=1:size(B)

A(i,j)=B(i,j)*X(j);

end

Z=X\'-sum(A)-D;

Y=(X\'-sum(A\'))\';

在Matlab命令窗口中输入：

la28

计算结果为：

036.000099.00000

90.000012.000018.000075.0000

54.000009.0000105.0000

36.000036.000027.000030.0000

1751411770

225

171

16.3实验习题

某地区经济系统最近周期的中间产品价值形成情况如表1.3所示，假设其计划周期的最

终产品向量为1000,1500,1000Y，那么计划周期的价值型投入产出表应如何编制？

表1.3价值型投入产出表（单位：亿元）

产出

投入

中间产品

最终产品Y总产

出X

1农业2工业3其他消费积累出口小计

生产资料

补偿价值

1农业6

2工业22426617403054

3其他474822375544

固定资

产折旧

115320107

新创造价值

劳动

报酬

520510414

纯收入330770470

小计8501280884

总投入X4

参考文献：

谢云荪、张志让．数学实验．北京：科学出版社，1998年

实验17Hill密码的加密、解密与破译

17.1实验指导

本实验研究了Hill密码的加密、解密与破译的数学模型及求解方法，当码值较长时，

借助数学软件Matlab可以高效地进行求解。

Hill密码是一种传统的密码体系，它的加密过程可以用以下流程来描述：











解码器明文

明文加密器密文普通信道

密码分析（被敌方截获）

在这个过程中，用到的数学手段是矩阵运算，加密过程的具体步骤如下：

a.根据明文字母的表值，将明文信息用数字表示，明文信息可以由26个字母A~Z表

示，也可以不止26个字母，可以含数字或符号．通信双方事先约定这26个字母的表值．

b.选择一个二阶可逆整数矩阵A，称为

Hill密码的加密矩阵，它是这个加密体系的

“密钥”．

c.将明文字母依次逐对分组，

Hill密码的加密矩阵为二阶矩阵，则明文字母2个一组

（

Hill

是

个字母分成一组），若最后一组字母数不足，则补充一些没有意义的哑字母，这

样使每一组都有2个字母组成，查出每个明文字母的表值，构成一个向量



．

d.A乘以，得一新的2维列向量A，由的两个分量反查字母表值得到的两

个字母即为密文字母．

以上四步即为

Hill密码的加密过程，解密过程为上述过程的逆过程．

也就是说如何在模意义下求解方程组

A（1）

在线性代数中，式（1）有唯一解的充要条件是矩阵A可逆，但模意义下矩阵逆的含义

却有所不同．记0,1,2,,1

GmL，

为一整数．下面给出模

下矩阵可逆的定义．

定义1对于一个元素属于集合

G的

阶方阵A，若存在一个元素属于

G的方阵B，

使得

(mod)ABBAEm（2）

则称A为模

可逆，B为A的模

的逆矩阵，记为1(mod)BAm．

定义2设

G的一个整数x，存在

G的一个整数y，使得1(mod)xym，则称y为

x的乘法逆（或者称为倒数），记1(mod)yxm．

可以证明，如果

与

无公共素数因子，则

有唯一的模

倒数．据此不加证明地给

出如下命题．

命题元素属于

G的方阵A模m可逆的充要条件是：m和detA没有公因子．

易见，所选加密矩阵必须符合模A可逆的条件.

前面的加密与解密过程类似与二维向量空间进行线性变换及其逆变换．每个明文向量是一个

G上的二维向量，乘以加密矩阵后，仍为

G上的二维向量．由于加密矩阵A为可逆矩阵，

所以知道了两个线性无关的二维明文向量与其对应的密文向量，就可以求出它的加密矩阵A

及1A．

17.2实验内容

例17.1（1）甲方收到与之有秘密通信往来的乙方的一个密文信息，密文内容如下：

AXSTZOSAOPBSTKSANKOPSAHAUUNSUUAKGAUZCKOPDO

按照甲方和乙方的约定，他们之间的密文采用

Hill密码，密钥为二阶矩阵







且26个英文字母与0~25的整数建立一一对应关系，称为字母表值，具体的表值见表17.1，

问这段密文的原意是什么？

表17.1字母表值

ABCDEFGHIJKLM

111213

NOPQRSTUVWXYZ

表17.2给出了模26的倒数表．

表17.2模26倒数表

19212325

1

31151725

（2）如果丙方截获了一段密文：

DWNINQNZAUUJBEIYILNHILONIYGKIYILN

经分析这段密文是用

Hill密码编译的，且这段密文中的字母ONIY对应字母OQIN，

能否破译这段密文的内容呢？

解：问题（1）所选的明文字母共有26个，将密文相邻2个字母分成一组，如下

AXSTZOSAOPBSTKSANKOPSAHAUUNSUUAKGA

UZCKOPDO根据表4.1列出的字母表值构造21组二维列向量：

1191519

,,,,,,,,,,

2421161







,…

3154

111615







。

在Matlab命令窗口输入：

A=[1,1;0,3];

a=det(A)%检验A是否可逆

alpha=mod(inv(A).*27,26)%对矩阵A的逆进行模26运算

Miwen=[1,24;19,20;0,15;19,1;15,16;2,19;20,11;19,1;14,11;15,16;19,

1;8,1;21,21;14,19;21,21;1,11;7,1;21,0;3,7;15,16;4,15]\';

%21组二维列向量组成的矩阵

Mingwen=zeros(2,21);

fori=1:21

Mingwen(:,i)=mod(alpha*Miwen(:,i),26);%求向量[1;24]对应的明文

end

结果为:

alpha=

117

Mingwen=

125

824599975

综上，得到了如下明文与密文对照表17.3：

表17.3明文与密文对照表

序

号

密文密文

表值

向量

明文序

号

密文密文

表值

向量

明文

1AX124198SH12HA81259YI

2ST19202124UX13UU2121147NG

3ZO015215UE14NS14192525YO

4SA191109JI15UU2121147NG

5OP1516114AN16AK111621FU

6BS2191315MO17GA71249XI

7TK20112521YU18UZ210210UZ

8SA191109JI19CK311821HU

9NK14111921SU20OP1516114AN

10OP1516114AN21DO415255YE

11SA191109JI

将译出的明文依据汉语拼音写出，经组合得到

SHUXUEJIANMOYUJISUANJIYINGYONGFUXIUZHUANYE

问题（2）属于破译问题．密文中只出现一些字母，当然它可以是汉语拼音，英文字母

或其他语言字母，所以可猜测秘密信息是由26个字母组成．设26m，通常由破译部门通

过大量的统计分析与语言分析确定表值．例如，所确定的表值就是表17.1，则有

OOII

NQYN









由表17.1得

111

1515

1417











222

2514











在模26的意义下，

159

,(mod26)5

1714

，它有模26的倒数，所以在模26

的意义下，

,是线性无关的．类似地，可以验证

,在模26下也是线性无关的．

借助线性代数的一些运算可以求得密钥







,其后的运算完全等价于第（1）问。

得这段密文的明文为DONGNANDAXUEBAINIANXIAOQINGJINIAN

17.3实验习题

1.（1）甲方收到与之有秘密通信往来的乙方的一个密文信息，密文内容如下：

WOWUYSBACPGZSAVCOVKPEWCPADKPPABUJCQLYXQEZAACPP

已知密钥为







能否知道这段密文的意思？

（2）甲方截获了一段密文：

OJWPISWAZUXAUUISEABAUCRSIPLBHAAMMLPJJOTENH

经分析这段密文是用

Hill密码编译的，且这段密文中的字母UCRS依次代表字母

TACO，问能否破译这段密文的内容呢？

参考文献：

乐经良．数学实验．北京：高等教育出版社，1999年．

👁️ 阅读量：0

🔖 本文标签：

🔥 最新发布文章

下三角矩阵的逆简单算法 矩阵单独一行乘一个数变吗