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.