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=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)    %求增广矩阵的行最简形，可得最简同解方程组。

运行后结果显示为：

C =

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

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

    0        0        0        0         0 

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

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)对应位置元素为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)

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