Matlab在线性代数中的应用

王涛 纪维强

【摘要】针对线性代数中比较典型的问题，文章借助于Matlab工具，结合具体的例子，展示Matlab在线性代数中的实际应用.

【关键词】Matlab；行列式；矩阵；特征值

线性代数中许多问题可以借助于Matlab软件来求解；本文结合线性代数中相关典型问题，给出了Matlab求解这些问题的相关用法，以供大家参考.

一、计算行列式

文中计算行列式的Matlab命令：det（A），其中A为方阵.

例1计算行列式axxxxaxxxxaxxxxa.

Matlab中输入：

clear

symsax；

A=[axxx；xaxx；xxax；xxxa]；

det（A）

得到：ans=

a^4-6*a^2*x^2+8*a*x^3-3*x^4

注：如果再输入：factor（det（A）），可得到因式分解形式下的行列式结果：

ans=

（a+3*x）*（a-x）^3

例2求方程1111123x149x21827x3=0的全部根.

Matlab中输入：

clear

symsx；

A=[1111；123x；149x^2；1827x^3]；

det（A）

得到：ans=

2*x^3-12*x^2+22*x-12

再输入：solve（‘2*x^3-12*x^2+22*x-12=0），可求得原方程的根为：

ans=

1

2

3

二、将矩阵化为行最简形矩阵

文中Matlab命令为：B=rref（A），rref（A）表示求A的行最简形矩阵B.

说明：根据文中思想，求出矩阵的行最简形后，就比较容易求出矩阵的秩、一个最高阶的非零子式、列（行）向量组的一个最大无关组及用最大无关组表示其余向量等等相关问题.

例3求矩阵A=1-130-21-21-1-152的行最簡形矩阵.

Matlab中输入：

clear

A=[1-130；-21-21；-1-152]；

B=rref（A）

得到：B=

10-1-101-4-10000

三、求逆矩阵

例4求A=1000120021301214的逆矩阵.

解法1根据文中相关知识，用行变换（A，E）：（E，A-1），将A化为行最简形矩阵，右端自然就出现A-1.

Matlab中输入：

clear

formatrat%设置显示格式为有理数

A=[1000；1200；2130；1214]；

B=[Aeye（4）]；%eye（m，n）为m×n的单位矩阵

C=rref（B）；

Ainv=C（：，5：8）

得到：Ainv=

1000-1/21/200-1/2-1/61/301/8-5/24-1/121/4

从而

A-1=1000-1/21/200-1/2-1/61/301/8-5/24-1/121/4.

解法2直接应用文中Matlab求逆矩阵命令：inv（A）也可得上述结果.

例5求A=abcd的逆矩阵（abcd≠0）.

Matlab中输入：

clear

symsabcd；

A=diag（[abcd]）；%构造对角形矩阵

B=inv（A）

得到：B=

[1/a，0，0，0]

[0，1/b，0，0]

[0，0，1/c，0]

[0，0，0，1/d]

说明对于一般形式的矩阵（若它是可逆的），都可以按照例4的方法来求解它的逆矩阵.

四、求解齐次线性方程组

在Matlab中，函数null用来求解零空间，即满足AX=0的解空间，实际上是求出解空间的一组基（基础解系）.基本格式：

z=null（A）%z的列向量为方程组AX=0的规范正交基，满足zTz=E；

z=null（A，′r′）%z的列向量为方程AX=0的有理基.

例6求齐次方程组x1+2x2-2x3+2x4-x5=0，x1+2x2-x3+3x4-2x5=02x1+4x2-7x3+x4+x5=0的通解.

Matlab中输入：

clear

formatrat

A=[12-22-1；12-13-2；24-711]；

B=null（A，r）

得到：B=

-2-431000-11010001

再输入：symsk1k2k3；

X=sym（B）*[k1；k2；k3]%或者输入X=k1*B（：，1）+k2*B（：，2）+k3*B（：，3）求出通解形式

得到通解形式：

X=

-2*k1-4*k2+3*k3

k1

-k2+k3

k2

k3

五、求解非齐次线性方程组

根据文的思路，非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解.因此，步骤为：

第一步：判断AX=b是否有解，若有解则进行第二步；

第二步：求AX=b的一个特解；

第三步：求对应齐次方程组AX=0的通解；

第四步：根据非齐次方程组通解结构（即AX=b的一个特解+对应齐次方程组AX=0的通解），求得通解形式.

例7求解非齐次方程组x1+2x2-x3+3x4=2，2x1+4x2-2x3+5x4=1，-x1-2x2+x3-x4=4.

Matlab中输入：

A=[12-13；24-25；-1-21-1]；

b=[214]；

B=[Ab]；

n=4；

R_A=rank（A）

R_B=rank（B）

formatrat

得到：R_A=

2R_B=

2

%根据方程组解的判定定理判定非齐次方程解的情形

再输入：symsk1k2；%齐次方程组的基础解系含有2个向量，选定2个自由常数

X=null（sym（A））*[k1；k2]+sym（Ab），E=A*X-b%此处X为通解形式，E为A*X与b的差向量值，目的在于验证X的求解是否准确

得到：Warning：Rankdeficient，rank=2，tol=5.2545e-015.

X=

-2*k1+k2

k1-7/2

k2

3

E=

0

0

0

%此处E为零向量，说明X为原方程组AX=b的精确解

另外文中，Matlab也可求解矩阵方程组，有如下命令：

①若矩阵方程形式为AX=B，在方程组有解的条件下，可用Matlab命令：X=AB求解；

②若矩阵方程形式为XA=B，在方程组有解的条件下，可用Matlab命令：X=B/A求解.

六、求矩阵的特征值和特征向量

文Matlab中求矩阵Am×n的特征值和特征向量的命令为：eig（A）或[V，D]=eig（A）

例8求矩阵A=1-333-536-64的特征值和特征向量.

解法1Matlab中输入：

clear

A=[1-33；3-53；6-64]；

eig（sym（A））

运行结果为：ans=

4

-2

-2

如果运行[V，D]=eig（sym（A））命令，得到：

V=

[1，-1，1]

[1，0，1]

[2，1，0]

D=

[4，0，0]

[0，-2，0]

[0，0，-2]

说明eig（A）仅显示A的特征值，而[V，D]=eig（A）不仅显示对角型矩阵D（对角线元素即为A的特征值），还求解出相应的特征向量构成的矩阵V.

解法2用求方程组基础解系的方法来求对应特征值的特征向量.

Matlab中输入：

clear

A=[1-33；3-53；6-64]；

eig（sym（A））

P1=sym（null（A-4*eye（3）））；%求（A-4E）x=0的基礎解系，即属于特征值λ1=4的线性无关的特征向量，sym允许含根号的形式

P2=sym（null（A+2*eye（3）））；%求属于特征值λ2=-2的线性无关的特征向量

P=[P1P2]%也可用disp（[P1P2]）

得到：ans=

4

-2

-2

P=

[sqrt（1/6），-sqrt（2/3），0]

[sqrt（1/6），-sqrt（1/6），-sqrt（1/2）]

[sqrt（2/3），sqrt（1/6），-sqrt（1/2）]

说明一般而言，解法2中求基础解系的方法，确定出来的特征向量与实际较为符合，误差较小.

七、求使得对称矩阵对角化的正交矩阵

例9求把A=22-225-4-2-45对角化的正交矩阵P

Matlab中输入：

clear

formatrat

A=[22-2；25-4；-2-45]；

f=poly（A）；%得到A的特征多项式f的标量形式

f=poly2sym（f）%得到A的特征多项式f的变量形式

solve（f）%求得特征多项式的根

运行结果为：f=

x^3-12*x^2+21*x-10

ans=

10

1

1

再运行：p1=sym（null（A-10*eye（3）））；

p2=sym（null（A-eye（3）））；

P=[p1p2]

得到：P=

[1/3，sqrt（8/9），0]

[2/3，-sqrt（1/18），sqrt（1/2）]

[-2/3，sqrt（1/18），sqrt（1/2）]

说明根据文中思想，借助于Matlab求得对称矩阵A的特征值后，可以比较容易判定出矩阵A的正定性；如果再进一步求出使得A对角化的正交矩阵P，则二次型f=xTAx采用正交变换x=Py化为标准形的问题也得到解决.

【参考文献】

[1]何正风.Matlab在数学方面的应用[M].北京：清华大学出版社，2012（1）：105-127.

[2]艾冬梅，刘琳，等.Matlab与数学实验[M].北京：机械工业出版社，2010.

[3]同济大学数学系，工程数学.线性代数[M].北京：高等教育出版社，2007（5）：124-133.