矩阵方程求解内置函数

特征值问题的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'：表示实数分解。

海森伯格形式的分解

如果矩阵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)

1.0000 0 0

0 -0.9987 0.0502

0 0.0502 0.9987

-149.0000 42.2037 -156.3165

-537.6783 152.5511 -554.9272

0 0.0728 2.4489

H的第一子对角元素是H(3,1)=0。

线性方程的组的求解

我们将线性方程的求解分为两类：一类是方程组求唯一解或求特解，另一类是方程组求无穷解即通解。可以通过系数矩阵的秩来判断：

若系数矩阵的秩r=n（n为方程组中未知变量的个数），则有唯一解；

若系数矩阵的秩r<n，则可能有无穷解；

线性方程组的无穷解 = 对应齐次方程组的通解+非齐次方程组的一个特解；其特解的求法属于解的第一类问题，通解部分属第二类问题。

求线性方程组的唯一解或特解（第一类问题）

这类问题的求法分为两类：一类主要用于解低阶稠密矩阵 —— 直接法；另一类是解大型稀疏矩阵 —— 迭代法。

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 =

2.2662

-1.7218

1.0571

-0.5940

0.3188

这就是方程组的解。

或用函数rref求解：

>> C=[A,B] %由系数矩阵和常数列构成增广矩阵C

>> R=rref(C) %将C化成行最简行

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)

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)

显示结果如下：

0.3636 -0.5000 1.0000

0.2727 1.0000 0

1.0000 0 0

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

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)

说明这三种分解，在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。

求线性齐次方程组的通解

在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') %求解空间的有理基

运行后显示结果如下：

2 5/3

-2 -4/3

1 0

0 1

或通过行最简行得到基：

>> B=rref(A)

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) %让通解表达式更加精美

运行后结果如下：

[ 2*k1+5/3*k2]

[ -2*k1-4/3*k2]

[ k1]

[ k2]

% 下面是其简化形式

[2k1 + 5/3k2 ]

[ ]

[-2k1 - 4/3k2]

[ ]

[ k1 ]

[ ]

[ k2 ]

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

非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解。

因此，步骤为：

第一步：判断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];

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' %判断无解

运行后结果显示：

R_A =

R_B =

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];

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'

运行后结果显示为：

R_A =

R_B =

Warning: Rank deficient, rank = 2 tol = 8.8373e-015.

> In D:\Matlab\pujun\lx0723.m at line 11

-8/15

3/2 -3/4

3/2 7/4

1 0

0 1

所以原方程组的通解为X=k ₁ +k ₂ +

解法二：用rref求解

A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];

b=[1 4 0]';

B=[A b];

C=rref(B) %求增广矩阵的行最简形，可得最简同解方程组。

运行后结果显示为：

1 0 -3/2 3/4 5/4

0 1 -3/2 -7/4 -1/4

0 0 0 0 0

对应齐次方程组的基础解系为：，非齐次方程组的特解为：所以，原方程组的通解为：X=k ₁ +k ₂ + 。

线性方程组的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表示每次迭代共轭梯度残差的范数

双共轭梯度法解方程组

函数 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。

iter = %表明解法不当，使得初始估计值0向量比后来所有迭代值都好。

resvec = （略） %每次迭代的残差。

>> semilogy(0:20,resvec/norm(b),'-o') %作每次迭代的相对残差图形，结果如下图。

>> xlabel('iteration number') %x轴为迭代次数。

>> ylabel('relative residual') %y轴为相对残差。

图1-1 双共轭梯度法相对误差图

稳定双共轭梯度方法解方程组

函数 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。

[x,flag,relres] = cgs(A,b,…)

[x,flag,relres,iter] = cgs(A,b,…)

[x,flag,relres,iter,resvec] = cgs(A,b,…)

调用方式和返回的结果形式与命令bicg一样。

共轭梯度的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 %计算终止时的迭代次数

广义最小残差法

函数 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时收敛。

最小残差法解方程组

函数 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 =

relres =

4.6537e-014

iter =

resvec =（略）

resveccg =（略）

预处理共轭梯度方法

函数 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为对称正定矩阵。

准最小残差法解方程组

函数 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 = %计算终止时的迭代次数

resvec2 = %每次迭代的残差

1.0e+005 *

7.0557

7.1773

3.4032

1.7220

0.0000

图1-3 准最小残差法相对误差图

特征值与二次型

工程技术中的一些问题，如振动问题和稳定性问题，常归结为求一个方阵的特征值和特征向量。

特征值与特征向量的求法

设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)

结果显示：

-0.7071 -0.2425 0.3015

0 0 0.9045

-0.7071 -0.9701 0.3015

-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)

结果显示为

0 0.4082 -0.4082

0 0.8165 -0.8165

1.0000 -0.4082 0.4082

2 0 0

0 1 0

0 0 1

说明当特征值为1 (二重根)时，对应特征向量都是k (0.4082 0.8165 -0.4082) ^T ，k为任意常数。

提高特征值的计算精度

函数 balance

格式 [T,B] = balance(A) %求相似变换矩阵T和平衡矩阵B，满足。

B = balance(A) %求平衡矩阵B

复对角矩阵转化为实对角矩阵

函数 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)

1.0000 -0.0191 - 0.4002i -0.0191 + 0.4002i

0 0 - 0.6479i 0 + 0.6479i

0 0.6479 0.6479

1.0000 0 0

0 4.0000 + 5.0000i 0

0 0 4.0000 - 5.0000i

>> [V,D]=cdf2rdf(v,d)

1.0000 -0.0191 -0.4002

0 0 -0.6479

0 0.6479 0

1.0000 0 0

0 4.0000 5.0000

0 -5.0000 4.0000

命令 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

则显示结果为

1.0000 0 0

0 0.7071 -0.7071

0 0.7071 0.7071

1.0000 0 0

0 1.0000 0

0 0 1.0000

例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

运行后结果显示如下：

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

1 0 0 0

0 1 0 0

0 0 -3 0

0 0 0 1

[ .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]

y1^2+y2^2-3*y3^2+y4^2

即 f =

秩与线性相关性

矩阵和向量组的秩以及向量组的线性相关性

矩阵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)

由于秩为3 < 向量个数，因此向量组线性相关。

求行阶梯矩阵及向量组的基

行阶梯使用初等行变换，矩阵的初等行变换有三条：

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]

1 -2 -1 0 2

-2 4 2 6 -6

2 -1 0 2 3

3 3 3 3 4

>> [R,jb]=rref(A)

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

1 2 4

>> A(:,jb)

ans =

1 -2 0

-2 4 6

2 -1 2

3 3 3

即：a1 a2 a4为向量组的一个基。

稀疏矩阵技术

稀疏矩阵的创建

函数 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)对应位置元素为s _i ，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)

(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)

(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

将稀疏矩阵转化为满矩阵

函数 full

格式 A=full(S) %S为稀疏矩阵，A为满矩阵。

例1-97

>> S=sparse(1:5,1:5,4:8)

(1,1) 4

(2,2) 5

(3,3) 6

(4,4) 7

(5,5) 8

>> A=full(S)

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

稀疏矩阵非零元素的索引

函数 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]’

外部数据转化为稀疏矩阵

函数 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]

1 2 3

2 5 4

3 4 6

3 6 7

>> S=spconvert(D)

(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];

1 2 3 4

2 5 4 0

3 4 6 9

3 6 7 4

>> S=spconvert(D)

(1,2) 3.0000 + 4.0000i

(3,4) 6.0000 + 9.0000i

(2,5) 4.0000

(3,6) 7.0000 + 4.0000i

基本稀疏矩阵

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)

41 11 0

52 22 0

63 33 13

74 44 24

-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．稀疏矩阵非零元素的个数

函数 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阶稀疏矩阵对角矩阵

(1,1) 1

(2,2) 2

(3,3) 3

(4,4) 4

f= spfun('exp',S) %即指数e的非零元素方幂。

(1,1) 2.7183

(2,2) 7.3891

(3,3) 20.0855

(4,4) 54.5982

6．把稀疏矩阵的非零元素全换为1

函数 spones

格式 R = spones(S) %将稀疏矩阵S中的非零元素全换为1

画稀疏矩阵非零元素的分布图形

函数 spy

格式 spy(S) %画出稀疏矩阵S中非零元素的分布图形。S也可以是满矩阵。

spy(S,markersize) % markersize为整数，指定点阵大小。

spy(S,'LineSpec') %'LineSpec'指定绘图标记和颜色

spy(S,'LineSpec',markersize) %参数与上面相同

函数 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]

11 0 13 0

41 0 0 44

0 22 0 24

0 0 63 0

>> [p,q,r]=dmperm(A)

3 2 1 4

2 4 1 3

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)

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 稀疏对称最小度排列图

稀疏矩阵的近似欧几里得范数和条件数

命令矩阵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相同。

稀疏矩阵的分解

命令不完全的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]

11 0 13 0

41 0 0 44

0 22 0 24

0 0 63 0

>> S=sparse(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')

(1,1) 1.0000

(4,1) 0.2683

(2,2) 1.0000

(3,3) 1.0000

(4,3) 0.2063

(4,4) 1.0000

(1,1) 41

(2,2) 22

(3,3) 63

(1,4) 44

(2,4) 24

(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无穷分解

稀疏矩阵的特征值分解

函数 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,…)

看后请点赞