ch4_2例题

1 4.1 矩阵和代数方程 4.1.1 矩阵运算和特征参数 一 矩阵运算 表4.2-1 矩阵运算含义及相应符号

术 语 数 学 含 义 MATLAB表达 矩阵加减 nmijijnmnmba][BA A+B,

A-B 标量与矩阵加减 nmijnmijnmnmbabaa

][][][B

a+B,

a-B

矩阵乘积 

nllmnmlkkjiknlijlmijnmijnmbabacBAC 1 ][ A*B

标量与矩阵相乘 nmijnmbaa][B a*B

矩阵的转置 HAB

IjkRjkIkjRkjaiabib)()()()(

A'

〖说明〗  在MATLAB中，“标量与矩阵加减运算”是一种“行业性”说法，在教科书中没有这种表述。本表格数学含义栏对“标量与矩阵加减运算”给出了完整的定义。  MATLAB的矩阵运算是在复数域定义的。在本表所列操作中，矩阵转置可全称为“矩阵的共轭转置（Hermition transpose）”，在数学表述上采用上标H表示。

【例 4.2-1】矩阵乘法算例。本例演示：矩阵乘法的含义；如何减少和避免循环；如何判断

两个双精度矩阵是否相等。 （1）实现矩阵乘的“标量程式”

clear rand('state',12) A=rand(2,4);B=rand(4,3); 2

%---------用三重循环体实现矩阵乘法--------- C1=zeros(size(A,1),size(B,2)); %为C预分配内存，可加快计算 <4> for ii=1:size(A,1) for jj=1:size(B,2) for k=1:size(A,2) C1(ii,jj)=C1(ii,jj)+A(ii,k)*B(k,jj); end end end % <13>

（2）实现矩阵乘的“saxpy程式”（saxpy是LAPACK术语） %---------用二重循环体实现矩阵乘法--------- C2=zeros(size(A,1),size(B,2)); for jj=1:size(B,2) for k=1:size(B,1) C2(:,jj)=C2(:,jj)+A(:,k)*B(k,jj); end end

（3）MATLAB中矩阵乘的最好实现方式 C3=A*B; %直接运用MATLAB的矩阵乘

（3）计算结果比较 C3_C3=norm(C3-C1,'fro') %采用计算矩阵F范数的方式比较两个矩阵之间的“差” C3_C2=norm(C3-C2,'fro')

〖说明〗  “标量程式”的效率最差，“saxpy”是MATLAB提高计算效率一系列措施中的一个，直接采用MATLAB中的矩阵乘算符“*”实现矩阵乘最好。  由于数值计算是有限精度计算，不同计算程式所得结果的比较一般采用矩阵范数进行。

【例 4.2-2】观察矩阵的转置操作和数组转置操作的差别。本例演示：共轭转置和非共轭转置；有理表示的显示格式

 p.4上的思考题 3

 本章习题都要做，下次课堂抽查。 二 矩阵的标量特征参数 表4.2-2 计算矩阵标量特征参数的指令 术 语 数 学 含 义 MATLAB指令

秩 Rank

可采用以下任何一种表述：  矩阵A中线性无关列（或行）向量组中最大向量数；  矩阵A中最高非零子行列式的阶数；  矩阵A中最高非奇异子矩阵的维数。 rank(A)

迹 Trace



),min(1nmiiia

，即矩阵主对角元素之

和。 trace(A)

行列式Determinant

njjjjnna1111)1(AA

，

式中j1A是元素ja1对应的子行列式 det(A)

【例4.2-3】矩阵标量特征参数计算示例。本例演示：rank, det, trace的使用；子行

列式的计算。

【例4.2-4】矩阵标量特征参数的性质。本例演示：两相乘矩阵位置交换迹不变；两

同阶矩阵相乘位置交换行列式不变。 format short g %采用默认显示格式

rand('state',0) %设置均匀分布随机发生器的初始状态 A=rand(3,3); %生成(3*3)随机阵 B=rand(3,3); %生成另一个(3*3)随机阵 C=rand(3,4); D=rand(4,3);

%任何符合矩阵乘法规则的两个矩阵的乘积的“迹”不变 tAB=trace(A*B) %两个同阶矩阵乘积 4

tBA=trace(B*A) tCD=trace(C*D) %两个“内维”相等矩阵的乘积 tDC=trace(D*C)

tAB =

3.6697 tBA = 3.6697 tCD = 2.1544 tDC = 2.1544

%同阶矩阵乘积行列式等于各矩阵行列式之乘积 d_A_B=det(A)*det(B) dAB=det(A*B) dBA=det(B*A)

d_A_B =

0.0846 dAB = 0.0846 dBA = 0.0846

%对于非同阶矩阵，两个矩阵的乘积行列式随相乘次序不同而不同，假如矩阵乘积存在的话。 %当(m*n)矩阵与(n*m)矩阵之积，或是(m*m)阵，或是(n*n)阵。 %若m>n，则积(m*m)阵一定行列式为零 dCD=det(C*D) dDC=det(D*C) %(m*n)

4.1.2 矩阵的变换和特征值分解 [R, ci] = rref(A) 借助初等变换把A变换成行阶梯矩阵R X = null(A) A矩阵零空间的全部正交基，满足0AX。 Z = orth(A) A矩阵值空间的全部正交基，满足)()(AZspanspan。 5

[V, D] = eig(A) A矩阵的特征值、特征向量分解，使VDAV。 〖说明〗  指令rref的输出量ci是个行数组。它的元素表示了矩阵A中线性独立“列”的序号。因此，length(ci)就是矩阵A的秩；A(:,ci)的所有列向量构成矩阵A的值空间。但length(ci)法计算矩阵秩远不如rank(A)准确、可靠；A(:,ci)决定的值空间也不如orth(A)可靠。 本指令主要用于教学目的，因为它的结果形式比较接近高校教科书提供的求秩方法。  指令eig只适用于特征值各异的矩阵。假如存在相同特征值，则应尝试采用实现Jordan分解的指令jordan。

【例4.2-7】简单实阵的特征值分解。本例演示：eig, cdf2rdf指令。 A=[1,-3;2,2/3];[V,D]=eig(A) V = 0.7746 0.7746 0.0430 - 0.6310i 0.0430 + 0.6310i D = 0.8333 + 2.4438i 0 0 0.8333 - 2.4438i

%把复数特征值对角阵D转换成实数块对角阵 [VR,DR]=cdf2rdf(V,D) %满足VR*DR/VR=A

【补例4.2-1】严谨的矩阵特征值分解指令condeig。本例演示：condeig的意

义和严谨性；特征条件数与eps1相当时，被分解矩阵可能存在“重特征值”，此时特征值分解不可信，而

应采用Jordan分解。 （1）特征值分解

A=gallery(5) %MATLAB提供的特殊矩阵，它具有五重特征值。 [V,D,c_eig]=condeig(A); %严谨型特征值分解指令，输出量c_eig是特征条件数 <2> diag(D)' %为紧凑地显示特征值而写 c_eig %若存在数量级接近或大于1e8的分量，说明特征值分解不可靠。

（2）Jordan分解 [VJ,DJ]=jordan(A); %求出准确的广义特征值，使A*VJ=VJ*D成立。 6

DJ

（3）分解结果验算 Ae=V*D/V %理论上应有1

VDVA

AJ=VJ*DJ/VJ %理论上应有1

JJJVDVA

【补例4.2-2】矩阵秩计算的几种方法。本例演示：指令的使用；比较几种求秩方法；

矩阵各组基所张空间相同。 （1）矩阵的分解

%format rat %有理格式显示 A=magic(4) %仅为产生一个试验矩阵 [R,ci]=rref(A) %行阶梯分解 Z=orth(A) %值空间 S=null(A) %零空间 [V,D]=eig(A) %特征值、特征向量分解

（2）确定矩阵秩的不同方法 r1=length(ci) %ci是A线性独立列的序号向量，独立列的数目即秩。 r2=size(Z,2) %A值空间所有基的数目就是秩 r3=rank(A) %采用更精确算法得到的秩。 r4=length(A)-size(S,2) %A的维数与零空间维数之差，就是秩 r5=sum(diag(D)~=0) %当A具有不同特征值时，其非零特征值数就是秩。 %本例由于计算截断误差，理论0特征值消失。 %由此判断出错误的矩阵秩。

（3）矩阵列向量组的基 % rref求出的ci给出了A中所有独立的列向量的序号 rank(A(:,ci))%验证ci所指示向量的独立性。这是A的一组基。

%Z和A(:,ci)所张空间相同 Z'*Z %理论上应是单位阵 rank([Z,A(:,ci)])