1.3.9 特征值問題的QZ分解
函數 qz
格式 [AA,BB,Q,Z,V] = qz(A,B) %A、B為方陣，產生上三角陣AA和BB，正交矩陣Q、Z或其列變換形式，V為特征向量陣。且滿足：Q*A*Z= AA 和Q*B*Z = BB。
[AA,BB,Q,Z,V] = qz(A,B,flag) %產生由flag決定的分解結果，flag取值為：'complex'：表示復數分解（默認），取值為'real'：表示實數分解。
1.3.10 海森伯格形式的分解
如果矩陣H的第一子對角線下元素都是0，則H為海森伯格(Hessenberg)矩陣。如果矩陣是對稱矩陣，則它的海森伯格形式是對角三角陣。MATLAB可以通過相似變換將矩陣變換成這種形式。
函數 hess
格式 H = hess(A) %返回矩陣A的海森伯格形式
[P,H] = hess(A) %P為酉矩陣，滿足：A = PHP' 且P'P = eye(size(A))。
例1-75
>> A=[-149 -50 -154;537 180 546;-27 -9 -25];
>> [P,H]=hess(A)
P =
1.0000 0 0
0 -0.9987 0.0502
0 0.0502 0.9987
H =
-149.0000 42.2037 -156.3165
-537.6783 152.5511 -554.9272
0 0.0728 2.4489
H的第一子對角元素是H(3,1)=0。
1.4 線性方程的組的求解
我們將線性方程的求解分為兩類：一類是方程組求唯一解或求特解，另一類是方程組求無窮解即通解。可以通過系數矩陣的秩來判斷：
若系數矩陣的秩r=n（n為方程組中未知變量的個數），則有唯一解；
若系數矩陣的秩r<n，則可能有無窮解；
線性方程組的無窮解 = 對應齊次方程組的通解+非齊次方程組的一個特解；其特解的求法屬於解的第一類問題，通解部分屬第二類問題。
1.4.1 求線性方程組的唯一解或特解（第一類問題）
這類問題的求法分為兩類：一類主要用於解低階稠密矩陣 —— 直接法；另一類是解大型稀疏矩陣 —— 迭代法。
1．利用矩陣除法求線性方程組的特解（或一個解）
方程：AX=b
解法：X=A\b
例1-76 求方程組 的解。
解：
>>A=[5 6 0 0 0
1 5 6 0 0
0 1 5 6 0
0 0 1 5 6
0 0 0 1 5];
B=[1 0 0 0 1]';
R_A=rank(A) %求秩
X=A\B %求解
運行后結果如下
R_A =
5
X =
2.2662
-1.7218
1.0571
-0.5940
0.3188
這就是方程組的解。
或用函數rref求解：
>> C=[A,B] %由系數矩陣和常數列構成增廣矩陣C
>> R=rref(C) %將C化成行最簡行
R =
1.0000 0 0 0 0 2.2662
0 1.0000 0 0 0 -1.7218
0 0 1.0000 0 0 1.0571
0 0 0 1.0000 0 -0.5940
0 0 0 0 1.0000 0.3188
則R的最后一列元素就是所求之解。
例1-77 求方程組 的一個特解。
解：
>>A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
>>B=[1 4 0]';
>>X=A\B %由於系數矩陣不滿秩，該解法可能存在誤差。
X =[ 0 0 -0.5333 0.6000]’（一個特解近似值）。
若用rref求解，則比較精確：
>> A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
B=[1 4 0]';
>> C=[A,B]; %構成增廣矩陣
>> R=rref(C)
R =
1.0000 0 -1.5000 0.7500 1.2500
0 1.0000 -1.5000 -1.7500 -0.2500
0 0 0 0 0
由此得解向量X=[1.2500 – 0.2500 0 0]’（一個特解）。
2．利用矩陣的LU、QR和cholesky分解求方程組的解
（1）LU分解：
LU分解又稱Gauss消去分解，可把任意方陣分解為下三角矩陣的基本變換形式（行交換）和上三角矩陣的乘積。即A=LU，L為下三角陣，U為上三角陣。
則：A*X=b 變成L*U*X=b
所以X=U\(L\b) 這樣可以大大提高運算速度。
命令 [L，U]=lu (A)
例1-78 求方程組 的一個特解。
解：
>>A=[4 2 -1;3 -1 2;11 3 0];
>>B=[2 10 8]';
>>D=det(A)
>>[L,U]=lu(A)
>>X=U\(L\B)
顯示結果如下：
D =
0
L =
0.3636 -0.5000 1.0000
0.2727 1.0000 0
1.0000 0 0
U =
11.0000 3.0000 0
0 -1.8182 2.0000
0 0 0.0000
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.018587e-017.
> In D:\Matlab\pujun\lx0720.m at line 4
X =
1.0e+016 *
-0.4053
1.4862
1.3511
說明 結果中的警告是由於系數行列式為零產生的。可以通過A*X驗證其正確性。
（2）Cholesky分解
若A為對稱正定矩陣，則Cholesky分解可將矩陣A分解成上三角矩陣和其轉置的乘積，即： 其中R為上三角陣。
方程 A*X=b 變成
所以
（3）QR分解
對於任何長方矩陣A，都可以進行QR分解，其中Q為正交矩陣，R為上三角矩陣的初等變換形式，即：A=QR
方程 A*X=b 變形成 QRX=b
所以 X=R\(Q\b)
上例中 [Q, R]=qr(A)
X=R\(Q\B)
說明 這三種分解，在求解大型方程組時很有用。其優點是運算速度快、可以節省磁盤空間、節省內存。
1.4.2 求線性齊次方程組的通解
在Matlab中，函數null用來求解零空間，即滿足A·X=0的解空間，實際上是求出解空間的一組基（基礎解系）。
格式 z = null % z的列向量為方程組的正交規范基，滿足 。
% z的列向量是方程AX=0的有理基
例1-79 求解方程組的通解：
解：
>>A=[1 2 2 1;2 1 -2 -2;1 -1 -4 -3];
>>format rat %指定有理式格式輸出
>>B=null(A,'r') %求解空間的有理基
運行后顯示結果如下：
B =
2 5/3
-2 -4/3
1 0
0 1
或通過行最簡行得到基：
>> B=rref(A)
B =
1.0000 0 -2.0000 -1.6667
0 1.0000 2.0000 1.3333
0 0 0 0
即可寫出其基礎解系（與上面結果一致）。
寫出通解：
syms k1 k2
X=k1*B(:,1)+k2*B(:,2) %寫出方程組的通解
pretty(X) %讓通解表達式更加精美
運行后結果如下：
X =
[ 2*k1+5/3*k2]
[ -2*k1-4/3*k2]
[ k1]
[ k2]
% 下面是其簡化形式
[2k1 + 5/3k2 ]
[ ]
[-2k1 - 4/3k2]
[ ]
[ k1 ]
[ ]
[ k2 ]
1.4.3 求非齊次線性方程組的通解
非齊次線性方程組需要先判斷方程組是否有解，若有解，再去求通解。
因此，步驟為：
第一步：判斷AX=b是否有解，若有解則進行第二步
第二步：求AX=b的一個特解
第三步：求AX=0的通解
第四步：AX=b的通解= AX=0的通解+AX=b的一個特解。
例1-80 求解方程組
解：在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<n %判斷有無窮解
X=A\b %求特解
C=null(A,'r') %求AX=0的基礎解系
else X='equition no solve' %判斷無解
end
運行后結果顯示：
R_A =
2
R_B =
3
X =
equition no solve
說明 該方程組無解
例1-81 求解方程組的通解：
解法一：在Matlab編輯器中建立M文件如下：
A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
b=[1 4 0]';
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<n
X=A\b
C=null(A,'r')
else X='Equation has no solves'
end
運行后結果顯示為：
R_A =
2
R_B =
2
Warning: Rank deficient, rank = 2 tol = 8.8373e-015.
> In D:\Matlab\pujun\lx0723.m at line 11
X =
0
0
-8/15
3/5
C =
3/2 -3/4
3/2 7/4
1 0
0 1
所以原方程組的通解為X=k1 +k2 +
解法二：用rref求解
A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
b=[1 4 0]';
B=[A b];
C=rref(B) %求增廣矩陣的行最簡形，可得最簡同解方程組。
運行后結果顯示為：
C =
1 0 -3/2 3/4 5/4
0 1 -3/2 -7/4 -1/4
0 0 0 0 0
對應齊次方程組的基礎解系為： ， 非齊次方程組的特解為： 所以，原方程組的通解為：X=k1 +k2 + 。
1.4.4 線性方程組的LQ解法
函數 symmlq
格式 x = symmlq(A,b) %求線性方程組AX=b的解X。A必須為n階對稱方陣，b為n元列向量。A可以是由afun定義並返回A*X的函數。如果收斂，將顯示結果信息；如果收斂失敗，將給出警告信息並顯示相對殘差norm(b-A*x)/norm(b)和計算終止的迭代次數。
symmlq(A,b,tol) %指定誤差tol，默認值是1e-6。
symmlq(A,b,tol,maxit) %maxit指定最大迭代次數
symmlq(A,b,tol,maxit,M) %M為用於對稱正定矩陣的預處理因子
symmlq(A,b,tol,maxit,M1,M2) %M=M1×M2
symmlq(A,b,tol,maxit,M1,M2,x0) %x0為初始估計值，默認值為0。
[x,flag] = symmlq(A,b,…) %flag的取值為：0表示在指定迭代次數之內按要求精度收斂；1表示在指定迭代次數內不收斂；2表示M為壞條件的預處理因子；3表示兩次連續迭代完全相同；4表示標量參數太小或太大；5表示預處理因子不是對稱正定的。
[x,flag,relres] = symmlq(A,b,…) % relres表示相對誤差norm(b-A*x)/norm(b)
[x,flag,relres,iter] = symmlq(A,b,…) %iter表示計算x的迭代次數
[x,flag,relres,iter,resvec] = symmlq(A,b,…) %resvec表示每次迭代的殘差：norm(b-A*x0)
[x,flag,relres,iter,resvec,resveccg] = symmlq(A,b,…) %resveccg表示每次迭代共軛梯度殘差的范數
1.4.5 雙共軛梯度法解方程組
函數 bicg
格式 x = bicg(A,b) %求線性方程組AX=b的解X。A必須為n階方陣，b為n元列向量。A可以是由afun定義並返回A*X的函數。如果收斂，將顯示結果信息；如果收斂失敗，將給出警告信息並顯示相對殘差norm(b-A*x)/norm(b)和計算終止的迭代次數。
bicg(A,b,tol) %指定誤差tol，默認值是1e-6。
bicg(A,b,tol,maxit) %maxit指定最大迭代次數
bicg(A,b,tol,maxit,M) %M為用於對稱正定矩陣的預處理因子
bicg(A,b,tol,maxit,M1,M2) %M=M1×M2
bicg(A,b,tol,maxit,M1,M2,x0) %x0為初始估計值，默認值為0。
[x,flag] = bicg(A,b,…) %flag的取值為：0表示在指定迭代次數之內按要求精度收斂；1表示在指定迭代次數內不收斂；2表示M為壞條件的預處理因子；3表示兩次連續迭代完全相同；4表示標量參數太小或太大。
[x,flag,relres] = bicg(A,b,…) % relres表示相對誤差norm(b-A*x)/norm(b)
[x,flag,relres,iter] = bicg(A,b,…) %iter表示計算x的迭代次數
[x,flag,relres,iter,resvec] = bicg(A,b,…) %resvec表示每次迭代的殘差：norm(b-A*x0)
例1-83 調用MATLAB6.0數據文件west0479。
>> load west0479
>> A=west0479; %將數據取為系數矩陣A。
>> b=sum (A,2); %將A的各行求和，構成一列向量。
>> X=A\b; %用“\”求AX=b的解。
>> 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
resvec = （略） %每次迭代的殘差。
>> semilogy(0:20,resvec/norm(b),'-o') %作每次迭代的相對殘差圖形，結果如下圖。
>> xlabel('iteration number') %x軸為迭代次數。
>> ylabel('relative residual') %y軸為相對殘差。
圖1-1 雙共軛梯度法相對誤差圖
1.4.6 穩定雙共軛梯度方法解方程組
函數 bicgstab
格式 x =bicgstab(A,b)
bicgstab(A,b,tol)
bicgstab(A,b,tol,maxit)
bicgstab(A,b,tol,maxit,M)
bicgstab(A,b,tol,maxit,M1,M2)
bicgstab(A,b,tol,maxit,M1,M2,x0)
[x,flag] = bicgstab(A,b,…)
[x,flag,relres] = bicgstab(A,b,…)
[x,flag,relres,iter] = bicgstab(A,b,…)
[x,flag,relres,iter,resvec] = bicgstab(A,b,…)
穩定雙共軛梯度法解方程組，調用方式和返回的結果形式和命令bicg一樣。
例1-84
>>load west0479;
>>A=west0479;
>>b=sum(A,2);
>>[x,flag]=bicgstab(A,b)
顯示結果是x的值全為0，flag=1。表示在默認誤差和默認迭代次數（20次）下不收斂。
若改為：
>>[L1,U1] = luinc(A,1e-5);
>>[x1,flag1] = bicgstab(A,b,1e-6,20,L1,U1)
即指定誤差，並用A的不完全LU分解因子L和U作為預處理因子M=L*U，其結果是x1的值全為0，flag=2表示預處理因子為壞條件的預處理因子。
若改為
>>[L2,U2]=luinc(A,1e-6); %稀疏矩陣的不完全LU分解。
>>[x2,flag2,relres2,iter2,resvec2]=bicgstab(A,b,1e-15,10,L2,U2)
%指定最大迭代次數為10次，預處理因子M=L*U。
>>semilogy(0:0.5:iter2,resvec2/norm(b),'-o') %每次迭代的相對殘差圖形，見圖1-2。
|
圖1-2 穩定雙共軛梯度方法的相對誤差圖 |
>>xlabel('iteration number')
>>ylabel('relative residual')
結果為
x2=（其值全為1，略）
flag2 =
0 %表示收斂
relres2 =
2.8534e-016 %收斂時的相對誤差
iter2 =
6 %計算終止時的迭代次數
resvec2 = %每次迭代的殘差
1.4.7 復共軛梯度平方法解方程組
函數 cgs
格式 x = cgs(A,b)
cgs(A,b,tol)
cgs(A,b,tol,maxit)
cgs(A,b,tol,maxit,M)
cgs(A,b,tol,maxit,M1,M2)
cgs(A,b,tol,maxit,M1,M2,x0)
[x,flag] = cgs(A,b,…)
[x,flag,relres] = cgs(A,b,…)
[x,flag,relres,iter] = cgs(A,b,…)
[x,flag,relres,iter,resvec] = cgs(A,b,…)
調用方式和返回的結果形式與命令bicg一樣。
1.4.8 共軛梯度的LSQR方法
函數 lsqr
格式 x = lsqr(A,b)
lsqr(A,b,tol)
lsqr(A,b,tol,maxit)
lsqr(A,b,tol,maxit,M)
lsqr(A,b,tol,maxit,M1,M2)
lsqr(A,b,tol,maxit,M1,M2,x0)
[x,flag] = lsqr(A,b,…)
[x,flag,relres] = lsqr(A,b,…)
[x,flag,relres,iter] = lsqr(A,b,…)
[x,flag,relres,iter,resvec] = lsqr(A,b,…)
調用方式和返回的結果形式與命令bicg一樣。
例1-85
>>n = 100;
>>on = ones(n,1);
>>A = spdiags([-2*on 4*on -on],-1:1,n,n); %產生一個對角矩陣
>>b = sum(A,2);
>>tol = 1e-8; %指定精度
>>maxit = 15; %指定最大迭代次數
>>M1 = spdiags([on/(-2) on],-1:0,n,n);
>>M2 = spdiags([4*on -on],0:1,n,n); %M1*M2=M，即產生預處理因子
>>[x,flag,relres,iter,resvec] = lsqr(A,b,tol,maxit,M1,M2,[])
結果顯示
x的值全為1
flag =
0 %表示收斂
relres =
3.5241e-009 %表示相對殘差
iter =
12 %計算終止時的迭代次數
1.4.9 廣義最小殘差法
函數 qmres
格式 x = gmres(A,b)
gmres(A,b,restart)
gmres(A,b,restart,tol)
gmres(A,b,restart,tol,maxit)
gmres(A,b,restart,tol,maxit,M)
gmres(A,b,restart,tol,maxit,M1,M2)
gmres(A,b,restart,tol,maxit,M1,M2,x0)
[x,flag] = gmres(A,b,…)
[x,flag,relres] = gmres(A,b,…)
[x,flag,relres,iter] = gmres(A,b,…)
[x,flag,relres,iter,resvec] = gmres(A,b,…)
除參數restart（缺省時為系數方陣A的階數n）可以給出外，調用方式和返回的結果形式與命令bicg一樣。
例1-86
>>load west0479;
>>A = west0479;
>>b = sum(A,2);
>>[x,flag] = gmres(A,b,5)
結果顯示flag=1，表示在默認精度和默認迭代次數下不收斂。
若最后一行改為
>>[L1,U1] = luinc(A,1e-5);
>>[x1,flag1] = gmres(A,b,5,1e-6,5,L1,U1)
結果flag1=2，說明該方法失敗，原因是U1的對角線上有0元素，構成壞條件的預處理因子。
若改為：
[L2,U2] = luinc(A,1e-6);
tol = 1e-15;
[x4,flag4,relres4,iter4,resvec4] = gmres(A,b,4,tol,5,L2,U2)
[x6,flag6,relres6,iter6,resvec6] = gmres(A,b,6,tol,3,L2,U2)
[x8,flag8,relres8,iter8,resvec8] = gmres(A,b,8,tol,3,L2,U2)
結果flag4=0，flag6=0，flag8=0表示參數restart分別取為4，6，8時收斂。
1.4.10 最小殘差法解方程組
函數 minres
格式 x = minres(A,b)
minres(A,b,tol)
minres(A,b,tol,maxit)
minres(A,b,tol.maxit,M)
minres(A,b,tol,maxit,M1,M2)
minres(A,b,tol,maxit,M1,M2,x0)
[x,flag] = minres(A,b,…)
[x,flag,relres] = minres(A,b,…)
[x,flag,relres,iter] = minres(A,b,…)
[x,flag,relres,iter,resvec] = minres(A,b,…)
[x,flag,relres,iter,resvec,resveccg] = minres(A,b,…)
這里A為對稱矩陣，這種方法是尋找最小殘差來求x。調用方式和返回的結果形式與命令bicg一樣。
例1-87
>>n = 100; on = ones(n,1);
>>A = spdiags([-2*on 4*on -2*on],-1:1,n,n);
>>b = sum(A,2);
>>tol = 1e-10;
>>maxit = 50;
>>M1 = spdiags(4*on,0,n,n);
>> [x,flag,relres,iter,resvec,resveccg] = minres(A,b,tol,maxit,M1,[],[])
結果顯示：
flag =
0
relres =
4.6537e-014
iter =
49
resvec =（略）
resveccg =（略）
1.4.11 預處理共軛梯度方法
函數 pcg
格式 x = pcg(A,b)
pcg(A,b,tol)
pcg(A,b,tol,maxit)
pcg(A,b,tol,maxit,M)
pcg(A,b,tol,maxit,M1,M2)
pcg(A,b,tol,maxit,M1,M2,x0)
pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)
[x,flag] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)
[x,flag,relres] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)
[x,flag,relres,iter] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)
[x,flag,relres,iter,resvec] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)
調用方式和返回的結果形式與命令bicg一樣，這里A為對稱正定矩陣。
1.4.12 准最小殘差法解方程組
函數 qmr
格式 x = qmr(A,b)
qmr(A,b,tol)
qmr(A,b,tol,maxit)
qmr(A,b,tol,maxit,M)
qmr(A,b,tol,maxit,M1,M2)
qmr(A,b,tol,maxit,M1,M2,x0)
qmr(afun,b,tol,maxit,m1fun,m2fun,x0,p1,p2,…)
[x,flag] = qmr(A,b,…)
[x,flag,relres] = qmr(A,b,…)
[x,flag,relres,iter] = qmr(A,b,…)
[x,flag,relres,iter,resvec] = qmr(A,b,…)
調用方式和返回的結果形式與命令bicg一樣，這里A為方陣。
例1-88
>>load west0479;
>>A = west0479;
>>b = sum(A,2);
>>[x,flag] = qmr(A,b)
結果顯示flag=1，表示在缺省情況下不收斂
若改為
>>[L1,U1] = luinc(A,1e-5);
>>[x1,flag1] = qmr(A,b,1e-6,20,L1,U1)
結果顯示 flag=2，表示是壞條件的預處理因子，不收斂。
若改為
>>[L2,U2] = luinc(A,1e-6);
>>[x2,flag2,relres2,iter2,resvec2] = qmr(A,b,1e-15,10,L2,U2)
>>semilogy(0:iter2,resvec2/norm(b),'-o') %每次迭代的相對殘差圖形
>>xlabel('iteration number')
>>ylabel('relative residual')
結果為
x=（全為1，略）
flag2 =
0 %表示收斂
relres2 = %表示相對殘差
2.8715e-016
iter2 = %計算終止時的迭代次數
8
resvec2 = %每次迭代的殘差
1.0e+005 *
7.0557
7.1773
3.4032
1.7220
0.0000
0.0000
0.0000
0.0000
0.0000
圖1-3 准最小殘差法相對誤差圖
1.5 特征值與二次型
工程技術中的一些問題，如振動問題和穩定性問題，常歸結為求一個方陣的特征值和特征向量。
1.5.1 特征值與特征向量的求法
設A為n階方陣，如果數“ ”和n維列向量x使得關系式 成立，則稱 為方陣A的特征值，非零向量x稱為A對應於特征值“ ”的特征向量。
詳見1.3.5和1.3.6節：特征值分解問題。
例1-89 求矩陣 的特征值和特征向量
解：
>>A=[-2 1 1;0 2 0;-4 1 3];
>>[V,D]=eig(A)
結果顯示：
V =
-0.7071 -0.2425 0.3015
0 0 0.9045
-0.7071 -0.9701 0.3015
D =
-1 0 0
0 2 0
0 0 2
即：特征值-1對應特征向量(-0.7071 0 -0.7071)T
特征值2對應特征向量(-0.2425 0 -0.9701)T和(-0.3015 0.9045 -0.3015)T
例1-90 求矩陣 的特征值和特征向量。
解：
>>A=[-1 1 0;-4 3 0;1 0 2];
>>[V,D]=eig(A)
結果顯示為
V =
0 0.4082 -0.4082
0 0.8165 -0.8165
1.0000 -0.4082 0.4082
D =
2 0 0
0 1 0
0 0 1
說明 當特征值為1 (二重根)時，對應特征向量都是k (0.4082 0.8165 -0.4082)T，k為任意常數。
1.5.2 提高特征值的計算精度
函數 balance
格式 [T,B] = balance(A) %求相似變換矩陣T和平衡矩陣B，滿足 。
B = balance(A) %求平衡矩陣B
1.5.3 復對角矩陣轉化為實對角矩陣
函數 cdf2rdf
格式 [V,D] = cdf2rdf (v,d) %將復對角陣d變為實對角陣D，在對角線上，用2×2實數塊代替共軛復數對。
例1-91
>> A=[1 2 3;0 4 5;0 -5 4];
>> [v,d]=eig(A)
v =
1.0000 -0.0191 - 0.4002i -0.0191 + 0.4002i
0 0 - 0.6479i 0 + 0.6479i
0 0.6479 0.6479
d =
1.0000 0 0
0 4.0000 + 5.0000i 0
0 0 4.0000 - 5.0000i
>> [V,D]=cdf2rdf(v,d)
V =
1.0000 -0.0191 -0.4002
0 0 -0.6479
0 0.6479 0
D =
1.0000 0 0
0 4.0000 5.0000
0 -5.0000 4.0000
1.5.4 正交基
命令 orth
格式 B=orth(A) %將矩陣A正交規范化，B的列與A的列具有相同的空間，B的列向量是正交向量，且滿足：B'*B = eye(rank(A))。
例1-92 將矩陣 正交規范化。
解：
>>A=[4 0 0; 0 3 1; 0 1 3];
>>B=orth(A)
>>Q=B'*B
則顯示結果為
P =
1.0000 0 0
0 0.7071 -0.7071
0 0.7071 0.7071
Q =
1.0000 0 0
0 1.0000 0
0 0 1.0000
1.5.5 二次型
例1-93 求一個正交變換X=PY，把二次型
化成標准形。
解：先寫出二次型的實對稱矩陣
在Matlab編輯器中建立M文件如下：
A=[0 1 1 -1;1 0 -1 1;1 -1 0 1;-1 1 1 0];
[P,D]=schur(A)
syms y1 y2 y3 y4
y=[y1;y2;y3;y4]；
X=vpa(P,2)*y %vpa表示可變精度計算，這里取2位精度
f=[y1 y2 y3 y4]*D*y
運行后結果顯示如下：
P =
780/989 780/3691 1/2 -390/1351
780/3691 780/989 -1/2 390/1351
780/1351 -780/1351 -1/2 390/1351
0 0 1/2 1170/1351
D =
1 0 0 0
0 1 0 0
0 0 -3 0
0 0 0 1
X =
[ .79*y1+.21*y2+.50*y3-.29*y4]
[ .21*y1+.79*y2-.50*y3+.29*y4]
[ .56*y1-.56*y2-.50*y3+.29*y4]
[ .50*y3+.85*y4]
f =
y1^2+y2^2-3*y3^2+y4^2
即 f =
1.6 秩與線性相關性
1.6.1 矩陣和向量組的秩以及向量組的線性相關性
矩陣A的秩是矩陣A中最高階非零子式的階數；向量組的秩通常由該向量組構成的矩陣來計算。
函數 rank
格式 k = rank(A) %返回矩陣A的行（或列）向量中線性無關個數
k = rank(A,tol) %tol為給定誤差
例1-94 求向量組 ， ， ， ， 的秩，並判斷其線性相關性。
>>A=[1 -2 2 3;-2 4 -1 3;-1 2 0 3;0 6 2 3;2 -6 3 4];
>>k=rank(A)
結果為
k =
3
由於秩為3 < 向量個數，因此向量組線性相關。
1.6.2 求行階梯矩陣及向量組的基
行階梯使用初等行變換，矩陣的初等行變換有三條：
1．交換兩行 （第i、第j兩行交換）
2．第i行的K倍
3．第i行的K倍加到第j行上去
通過這三條變換可以將矩陣化成行最簡形，從而找出列向量組的一個最大無關組，Matlab將矩陣化成行最簡形的命令是rref或rrefmovie。
函數 rref或rrefmovie
格式 R = rref(A) %用高斯—約當消元法和行主元法求A的行最簡行矩陣R
[R,jb] = rref(A) %jb是一個向量，其含義為：r = length(jb)為A的秩；A(:, jb)為A的列向量基；jb中元素表示基向量所在的列。
[R,jb] = rref(A,tol) %tol為指定的精度
rrefmovie(A) %給出每一步化簡的過程
例1-95 求向量組a1=(1,-2,2,3)，a2=(-2,4,-1,3)，a3=(-1,2,0,3)，a4=(0,6,2,3)，a5=(2,-6,3,4)的一個最大無關組。
>> a1=[1 -2 2 3]';
>>a2=[-2 4 -1 3]';
>>a3=[-1 2 0 3]';
>>a4=[0 6 2 3]';
>>a5=[2 -6 3 4]';
A=[a1 a2 a3 a4 a5]
A =
1 -2 -1 0 2
-2 4 2 6 -6
2 -1 0 2 3
3 3 3 3 4
>> [R,jb]=rref(A)
R =
1.0000 0 0.3333 0 1.7778
0 1.0000 0.6667 0 -0.1111
0 0 0 1.0000 -0.3333
0 0 0 0 0
jb =
1 2 4
>> A(:,jb)
ans =
1 -2 0
-2 4 6
2 -1 2
3 3 3
即：a1 a2 a4為向量組的一個基。
1.7 稀疏矩陣技術
1.7.1 稀疏矩陣的創建
函數 sparse
格式 S = sparse(A) %將矩陣A轉化為稀疏矩陣形式，即由A的非零元素和下標構成稀疏矩陣S。若A本身為稀疏矩陣，則返回A本身。
S = sparse(m,n) %生成一個m×n的所有元素都是0的稀疏矩陣
S = sparse(i,j,s) %生成一個由長度相同的向量i，j和s定義的稀疏矩陣S，其中i，j是整數向量，定義稀疏矩陣的元素位置(i,j)，s是一個標量或與i，j長度相同的向量，表示在(i,j)位置上的元素。
S = sparse(i,j,s,m,n) %生成一個m×n的稀疏矩陣，(i,j)對應位置元素為si，m = max(i)且n =max(j)。
S = sparse(i,j,s,m,n,nzmax) %生成一個m×n的含有nzmax個非零元素的稀疏矩陣S，nzmax的值必須大於或者等於向量i和j的長度。
例1-96
>> S=sparse(1:10,1:10,1:10)
S =
(1,1) 1
(2,2) 2
(3,3) 3
(4,4) 4
(5,5) 5
(6,6) 6
(7,7) 7
(8,8) 8
(9,9) 9
(10,10) 10
>> S=sparse(1:10,1:10,5)
S =
(1,1) 5
(2,2) 5
(3,3) 5
(4,4) 5
(5,5) 5
(6,6) 5
(7,7) 5
(8,8) 5
(9,9) 5
(10,10) 5
1.7.2 將稀疏矩陣轉化為滿矩陣
函數 full
格式 A=full(S) %S為稀疏矩陣，A為滿矩陣。
例1-97
>> S=sparse(1:5,1:5,4:8)
S =
(1,1) 4
(2,2) 5
(3,3) 6
(4,4) 7
(5,5) 8
>> A=full(S)
A =
4 0 0 0 0
0 5 0 0 0
0 0 6 0 0
0 0 0 7 0
0 0 0 0 8
1.7.3 稀疏矩陣非零元素的索引
函數 find
格式 k = find(x) %按行檢索X中非零元素的點，若沒有非零元素，將返回空矩陣。
[i,j] = find(X) %檢索X中非零元素的行標i和列標j
[i,j,v] = find(X) %檢索X中非零元素的行標i和列標j以及對應的元素值v
例1-98 上例中
>> [i,j,v]=find(S)
則顯示為：
I =[ 1 2 3 4 5]’
j =[ 1 2 3 4 5]’
v =[4 5 6 7 8]’
1.7.4 外部數據轉化為稀疏矩陣
函數 spconvert
格式 S=spconvert(D) %D是只有3列或4列的矩陣
說明：先運用load函數把外部數據（.mat文件或.dat文件）裝載於MATLAB內存空間中的變量T；T數組的行維為nnz或nnz+1，列維為3（對實數而言）或列維為4（對復數而言）；T數組的每一行（以[i,j,Sre,Sim]形式）指定一個稀疏矩陣元素。
例1-99
>> D=[1 2 3;2 5 4;3 4 6;3 6 7]
D =
1 2 3
2 5 4
3 4 6
3 6 7
>> S=spconvert(D)
S =
(1,2) 3
(3,4) 6
(2,5) 4
(3,6) 7
>> D=[1 2 3 4; 2 5 4 0;3 4 6 9;3 6 7 4];
D =
1 2 3 4
2 5 4 0
3 4 6 9
3 6 7 4
>> S=spconvert(D)
S =
(1,2) 3.0000 + 4.0000i
(3,4) 6.0000 + 9.0000i
(2,5) 4.0000
(3,6) 7.0000 + 4.0000i
1.7.5 基本稀疏矩陣
1．帶狀（對角）稀疏矩陣
函數 spdiags
格式 [B,d] = spdiags(A) %從矩陣A中提取所有非零對角元素，這些元素保存在矩陣B中，向量d表示非零元素的對角線位置。
B = spdiags(A,d) %從A中提取由d指定的對角線元素，並存放在B中。
A = spdiags(B,d,A) %用B中的列替換A中由d指定的對角線元素，輸出稀疏矩陣。
A = spdiags(B,d,m,n) %產生一個m×n稀疏矩陣A，其元素是B中的列元素放
在由d指定的對角線位置上。
例1-100
>>A = [11 0 13 0
0 22 0 24
0 0 33 0
41 0 0 44
0 52 0 0
0 0 63 0
0 0 0 74];
>>[B,d] = spdiags(A)
B =
41 11 0
52 22 0
63 33 13
74 44 24
d =
-3 %表示B的第1列元素在A中主對角線下方第3條對角線上
0 %表示B的第2列在A的主對角線上
2 %表示B的第3列在A的主對角線上方第2條對角線上
例1-101
>> B=[1 2 3 4
5 6 7 8
9 10 11 12
13 14 15 16];
>> d=[-2 0 1 3];
>> A=spdiags(B,d,4,4);
>> full(A)
ans =
2 7 0 16
0 6 11 0
1 0 10 15
0 5 0 14
2．單位稀疏矩陣
函數 speye
格式 S = speye(m,n) %生成m×n的單位稀疏矩陣
S = speye(n) %生成n×n的單位稀疏矩陣
3．稀疏均勻分布隨機矩陣
函數 sprand
格式 R = sprand(S) %生成與S具有相同稀疏結構的均勻分布隨機矩陣
R = sprand(m,n,density) %生成一個m×n的服從均勻分布的隨機稀疏矩陣，非零元素的分布密度是density。
R = sprand(m,n,density,rc) %生成一個近似的條件數為1/rc、大小為m×n的均勻分布的隨機稀疏矩陣。
4．稀疏正態分布隨機矩陣
函數 sprandn
格式 R = sprandn(S) %生成與S具有相同稀疏結構的正態分布隨機矩陣。
R = sprandn(m,n,density) %生成一個m×n的服從正態分布的隨機稀疏矩陣，非零元素的分布密度是density。
R = sprandn(m,n,density,rc) %生成一個近似的條件數為1/rc、大小為m×n的均勻分布的隨機稀疏矩陣。
5．稀疏對稱隨機矩陣
函數 sprandsym
格式 R = sprandsym(S) %生成稀疏對稱隨機矩陣，其下三角和對角線與S具有相同的結構，其元素服從均值為0、方差為1的標准正態分布。
R = sprandsym(n,density) %生成n×n的稀疏對稱隨機矩陣，矩陣元素服從正態分布，分布密度為density。
R = sprandsym(n,density,rc) %生成近似條件數為1/rc的稀疏對稱隨機矩陣
R = sprandsym(n,density,rc,kind) %生成一個正定矩陣，參數kind取值為kind=1表示矩陣由一正定對角矩陣經隨機Jacobi旋轉得到，其條件數正好為1/rc；kind=2表示矩陣為外積的換位和，其條件數近似等於1/rc；kind=3表示生成一個與矩陣S結構相同的稀疏隨機矩陣，條件數近似為1/rc ，density被忽略。
1.7.6 稀疏矩陣的運算
1．稀疏矩陣非零元素的個數
函數 nnz
格式 n = nnz(X) %返回矩陣X中非零元素的個數
2．稀疏矩陣的非零元素
函數 nonzeros
格式 s = nonzeros(A) %返回矩陣A中非零元素按列順序構成的列向量
例1-102
>>A =[ 2 7 0 16
0 6 11 0
1 0 10 15
0 5 0 14];
>> s=nonzeros(A)
結果為：
s =[ 2 1 7 6 5 11 10 16 15 14]’
3．稀疏矩陣非零元素的內存分配
函數 nzmax
格式 n = nzmax(S) %返回非零元素分配的內存總數n
4．稀疏矩陣的存貯空間
函數 spalloc
格式 S = spalloc(m,n,nzmax) %產生一個m×n階只有nzmax個非零元素的稀疏矩陣，這樣可以有效減少存貯空間和提高運算速度。
5．稀疏矩陣的非零元素應用
函數 spfun
格式 f = spfun('function',S) %用S中非零元素對函數'function'求值，如果'function'不是對稀疏矩陣定義的，同樣可以求值。
例1-103 4階稀疏矩陣對角矩陣
S =
(1,1) 1
(2,2) 2
(3,3) 3
(4,4) 4
f= spfun('exp',S) %即指數e的非零元素方冪。
結果為
f =
(1,1) 2.7183
(2,2) 7.3891
(3,3) 20.0855
(4,4) 54.5982
6．把稀疏矩陣的非零元素全換為1
函數 spones
格式 R = spones(S) %將稀疏矩陣S中的非零元素全換為1
1.7.7 畫稀疏矩陣非零元素的分布圖形
函數 spy
格式 spy(S) %畫出稀疏矩陣S中非零元素的分布圖形。S也可以是滿矩陣。
spy(S,markersize) % markersize為整數，指定點陣大小。
spy(S,'LineSpec') %'LineSpec'指定繪圖標記和顏色
spy(S,'LineSpec',markersize) %參數與上面相同
|
圖1-4 稀疏矩陣非零元素的分布圖形 |
例1-104
>> load west0479
>> A=west0479;
>> spy(A,'ro',3)
結果如圖1-4所示。
1.7.8 矩陣變換
1．列近似最小度排序
函數 colamd
格式 p = colamd(S) %返回稀疏矩陣S的列的近似最小度排序向量p
2．列最小度排序
函數 colmmd
格式 p = colmmd(S) %返回稀疏矩陣S的列的最小度排序向量p，按p排列后的矩陣為S(: , p)。
例1-105 比較稀疏矩陣S與排序后的矩陣S(: , p)
>> load west0479;
>> S=west0479;
>> p=colmmd(S);
>> subplot(2,2,1),spy(S)
>> subplot(2,2,2),spy(S(:,p))
>> subplot(2,2,3),spy(lu(S))
>> subplot(2,2,4),spy(lu(S(:,p)))
圖1-5 稀疏矩陣的排序圖
3．非零元素的列變換
函數 colperm
格式 j = colperm(S) %返回一個稀疏矩陣S的列變換的向量。列按非0元素升序排列。有時這是LU分解前有用的變換：LU(S(:,j))。如果S是一個對稱矩陣，對行和列進行排序，有利於Cholesky分解：chol(S(j,j))
4．Dulmage-Mendelsohn分解
函數 dmperm
格式 p = dmperm (A) %返回A的行排列向量p，這樣，如果A滿列秩，就使得A(p,:)是具有非0對角線元素的方陣。
[p,q,r] = dmperm(A) %A為方陣，p為行排列向量，q為列排列向量，使得A(p,q)
是上三角塊形式，r為索引向量。
[p,q,r,s] = dmperm(A) %A不是方陣，p,q,r含義與前面相同，s也是索引向量。
例1-106
>> A=[11 0 13 0;41 0 0 44;0 22 0 24;0 0 63 0]
A =
11 0 13 0
41 0 0 44
0 22 0 24
0 0 63 0
>> [p,q,r]=dmperm(A)
p =
3 2 1 4
q =
2 4 1 3
r =
1 2 3 4 5
>> A(p,q)
ans =
22 24 0 0
0 44 41 0
0 0 11 13
0 0 0 63
5．整數的隨機排列
函數 randperm
格式 p = randperm (n) %對正整數1，2，3，…，n的隨機排列，可以用來創建隨機變換矩陣。
例1-107
>> p=randperm(6)
p =
3 4 6 5 1 2
6．稀疏對稱近似最小度排列
函數 symamd
格式 p = symamd(S) %S為對稱正定矩陣，返回排列向量p。
7．稀疏逆Cuthill-McKee排序
函數 symrcm
格式 r = symrcm (S) %返回S的對稱逆Cuthill-McKee排序r，使S的非0元素集中在主對角線附近。
8．稀疏對稱最小度排列
函數 symmmd
格式 p = symmmd(S) %返回S的對稱最小度排列向量p，S為對稱正定矩陣。
例1-108
>> B = bucky+4*speye(60);
>>r = symrcm(B);
>>p = symmmd(B);
>>R = B(r,r);
>>S = B(p,p);
>>subplot(2,2,1), spy(R), title('B(r,r)')
>>subplot(2,2,2), spy(S), title('B(s,s)')
>>subplot(2,2,3), spy(chol(R)), title('chol(B(r,r))')
>>subplot(2,2,4), spy(chol(S)), title('chol(B(s,s))')
圖1-6 稀疏對稱最小度排列圖
1.7.9 稀疏矩陣的近似歐幾里得范數和條件數
命令 矩陣A條件數的1-范數估計值
函數 condest
格式 c = condest(A) %計算方陣A的1-范數中條件數的下界值c
[c,v] = condest(A) %方陣A的1-范數中條件數的下界值c和向量v，使得 ，即：
norm(A*v,1) = norm(A,1)*norm(v,1)/c。
命令 2-范數估計值
函數 normest
格式 nrm = normest(S) %返回矩陣S的2-范數的估計值，相對誤差為10-6。
nrm = normest(S,tol) %tol為指定的相對誤差，而不用默認誤差10-6。
[nrm,count] = normest(…) %count為給出的計算范數迭代的次數
其條件數的計算與前面1.2.12相同。
1.7.10 稀疏矩陣的分解
命令 不完全的LU分解
函數 luinc
格式 [L,U] = luinc(X,'0') %X為稀疏方陣；L為下三角矩陣的置換形式；U為上三角矩陣；'0'是一種分解標准。
[L,U,P] = luinc(X,'0') %L為下三角陣，其主對角線元素為1；U為上三角陣，p為單位矩陣的置換形式。
[L,U] = luinc(X,options) %options取值為：droptol表示指定的舍入誤差；milu表示改變分解以便於從上三角分解因子中抽取被去掉的列元素。ugiag為1表示用droptol值代替上三角因子中的對角線上的零元素，默認值為0。thresh為中心臨界值。
[L,U] = luinc(X,droptol) %droptol表示指定不完全分解的舍入誤差
[L,U,P] = luinc(X,options)
[L,U,P] = luinc(X,droptol)
例1-109
>> S=[11 0 13 0;41 0 0 44;0 22 0 24;0 0 63 0]
S =
11 0 13 0
41 0 0 44
0 22 0 24
0 0 63 0
>> S=sparse(S)
S =
(1,1) 11
(2,1) 41
(3,2) 22
(1,3) 13
(4,3) 63
(2,4) 44
(3,4) 24
>> luinc(S,'0')
ans =
(1,1) 41.0000
(4,1) 0.2683
(2,2) 22.0000
(3,3) 63.0000
(4,3) 0.2063
(1,4) 44.0000
(2,4) 24.0000
>> [L,U,p]=luinc(S,'0')
L =
(1,1) 1.0000
(4,1) 0.2683
(2,2) 1.0000
(3,3) 1.0000
(4,3) 0.2063
(4,4) 1.0000
U =
(1,1) 41
(2,2) 22
(3,3) 63
(1,4) 44
(2,4) 24
p =
(4,1) 1
(1,2) 1
(2,3) 1
(3,4) 1
命令 稀疏矩陣的不完全Cholesky分解
函數 cholinc
格式 R = (X,droptol) %稀疏矩陣X的不完全Cholesky分解，droptol為指定誤差。
R = cholinc(X,options) %options取值為：droptol表示舍入誤差；michol表示如果michol=1，則從對角線上抽取出被去掉的元素。rdiag表示用droptol值代替上三角分解因子中的對角線上的零元素。
R = cholinc(X,'0') %'0'是一種分解標准
[R,p] = cholinc(X,'0') %不產生任何出錯信息，如果R存在，則p=0；如果R不存在，則p為正整數。
R = cholinc(X,'inf') %進行Cholesky無窮分解
1.7.11 稀疏矩陣的特征值分解
函數 eigs
格式 d = eigs(A) %求稀疏矩陣A的6個最大特征值d，d以向量形式存放。
d = eigs(A,B) %求稀疏矩陣的廣義特征值問題。滿足AV=BVD，其中D為特征值對角陣，V為特征向量矩陣，B必須是對稱正定陣或Hermitian正定陣。
d = eigs(A,k) %返回k個最大特征值
d = eigs(A,B,k) %返回k個最大特征值
d = eigs(A,k,sigma) %sigma取值：'lm' 表示最大數量的特征值；'sm' 最小數量特征值；對實對稱問題：'la'表示最大特征值；'sa'為最小特征值；對非對稱和復數問題：'lr' 表示最大實部；'sr' 表示最小實部；'li' 表示最大虛部；'si'表示最小虛部。
d = eigs(A,B,k,sigma) %同上
d = eigs(A,k,sigma,options) % options為指定參數：參見eigs幫助文件。
d = eigs(A,B,k,sigma,options) %同上。以下的參數k、sigma、options相同。
d = eigs(Afun,n) %用函數Afun代替A，n為A的階數，D為特征值。
d = eigs(Afun,n,B)
d = eigs(Afun,n,k)
d = eigs(Afun,n,B,k)
d = eigs(Afun,n,k,sigma)
d = eigs(Afun,n,B,k,sigma)
d = eigs(Afun,n,k,sigma,options)
d = eigs(Afun,n,B,k,sigma,options)
[V,D] = eigs(A,…) %D為6個最大特征值對角陣，V的列向量為對應特征向量。
[V,D] = eigs(Afun,n,…)
[V,D,flag] = eigs(A,…) %flag表示特征值的收斂性，若flag=0，則所有特征值都收斂，否則，不是所有都收斂。
[V,D,flag] = eigs(Afun,n,…)
1.7.12 稀疏矩陣的線性方程組
參見1.4節的方程組求解。