• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号

MATLAB矩阵运算(3)

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

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   


鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
上一篇:
delphi批量更改文件名发布时间:2022-07-22
下一篇:
delphi视频播放器发布时间:2022-07-22
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap