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节的方程组求解。
