《数字信号处理》教学参考资料（MATLAB 5手册）第7章 线性方程系统

China-bub.com 下载 第7章 线性方程系统 线性方程系统是最常见的计算问题，几乎在所有的应用中都把它作为子问题提出来。通常 MATLAB用运算符1从左分开来求解线性方程系统，超定系统的解决方法同样可用于欠定系统。 在线性系统理论中重要的概念是行列式、逆和矩阵的秩。首先定义这些MATLAB命令， 再在72节中开始介绍求解的方法。在范数和条件数定义之后，介绍一些因数分解。最后一节 涉及超定和欠定系统。 注意，本章中的所有命令只能用于二维矩阵。 7.1行列式、逆和秩 下列命令用来计算矩阵A的行列式、逆和矩阵的秩。 命令集67 矩阵函数 det(A) 求方阵A的行列式。 rank(A) 求A的秩，即A中线性无关的行数和列数。 inv(A) 求方阵A的逆矩阵。如果A是奇异矩阵或者近似奇异矩阵，则会给 出一个错误信息。 pinv(A) 求矩阵A的伪逆。如果A是m×n的矩阵，则伪逆的大小为n×m。对 于非奇矩阵A来说，有pinv(A)=inv(A)· trace(A) 求矩阵A的迹，也就是对角线元素之和。 ■例7.1 假设有下列矩阵： A1=(及42=(（3)A3-(3；) 将命令det、inv、rank和一些其他命令对上述矩阵进行操作。 (a)det1 det(A1),det2 det(A2),det3 det(A3) give: det1 -2 det2 0 ??Error using ==det Matrix must be square

下载 第7章 线性方程系统 线性方程系统是最常见的计算问题，几乎在所有的应用中都把它作为子问题提出来。通常 M AT L A B用运算符\从左分开来求解线性方程系统，超定系统的解决方法同样可用于欠定系统。 在线性系统理论中重要的概念是行列式、逆和矩阵的秩。首先定义这些 M AT L A B命令， 再在7 . 2节中开始介绍求解的方法。在范数和条件数定义之后，介绍一些因数分解。最后一节 涉及超定和欠定系统。 注意，本章中的所有命令只能用于二维矩阵。 7.1 行列式、逆和秩 下列命令用来计算矩阵A的行列式、逆和矩阵的秩。 命令集6 7 矩阵函数 d e t ( A ) 求方阵A的行列式。 r a n k ( A ) 求A的秩，即A中线性无关的行数和列数。 i n v ( A ) 求方阵A的逆矩阵。如果 A是奇异矩阵或者近似奇异矩阵，则会给 出一个错误信息。 p i n v ( A ) 求矩阵A的伪逆。如果A是m×n的矩阵，则伪逆的大小为 n×m。对 于非奇矩阵A来说，有p i n v ( A ) = i n v ( A )。 t r a c e ( A ) 求矩阵A的迹，也就是对角线元素之和。 ■ 例7 . 1 假设有下列矩阵： 将命令d e t、i n v、r a n k和一些其他命令对上述矩阵进行操作

96 MATLAB5手册 China-pub.com 下载 仅为方阵定义行列式。 (b)仅为方阵定义逆： Inv1 inv(A1),Inv2 =inv(A2),Inv3 inv(A3) 结果为： Invi -2.0000 1.5000 1.0000 -0.5000 Warning:Matrix is singular to working precision. Inv2 Inf Inf Inf Inf ??Error using ==inv Matrix must be square. 对所有矩阵定义伪逆： Pinv1 pinv(A1),Pinv2 pinv(A2),Pinv3 pinv(A3) 结果为： Pinv1 -2.0000 1.5000 1.0000 -0.5000 Pinv2 0.0200 0.0400 0.0600 0.1200 Pinv3 0.9048 -0.3333 1.0476 -0.3333 -1.5238 0.6667 注意，A1的逆矩阵和它的伪逆是一样的。 1 (C)如果A的逆矩阵存在，那么它的行列式det(A-)就等于： det(4) detinv1 det(inv(A1)) detinv1 -0.5000 (矩阵的秩和它的转置的秩相同： ranki rank(A1),rank2 rank(A2),rank3 rank(A3) 结果为： rank1 2 rank2 1

仅为方阵定义行列式。 (b) 仅为方阵定义逆： 结果为： 对所有矩阵定义伪逆： 结果为： 注意，A 1的逆矩阵和它的伪逆是一样的。 (c) 如果A的逆矩阵存在，那么它的行列式 d e t ( A－1 )就等于： (d) 矩阵的秩和它的转置的秩相同： 结果为： 1 det(A) 9 6 M ATLAB 5 手册 下载

China-pub.Com 第7章线性方程系统 下载 rankT3 2 和 rankT1 =rank(A1'),rankT2 rank(A2'),rankT3 rank(A3') 结果为： rankT1 2 rankT2 1 rankT3= 2 (©)实数矩阵的行列式和它的转置的行列式相同： detT1 det(A1'),detT2 det(A2,),detT3 =det(A3') detT1= -2 detT2 0 ??Error using ==>det Matrix must be square. 与线性系统相联系的两个子空间是值域和零空间。如果A为m×n的矩阵，它的秩为r,那 么A的向量空间就是由A的列划分的线性空间，这个空间的维数是r,也就是A的秩。如果=， 则A的列线性无关。MATLAB命令orth用来求A的空间的正交基。 A的零空间是由满足条件Ax=0的所有向量x组成的线性子空间。在MATLAB中可以用命令 null来求得零空间的正交基。 假设有一个向量集v,V,V。,可以通过定义矩阵B=(v,V,…v)来判断它们是否线性相 关。例如，如果B的秩是n一1，那么其中的一个向量v可以用其他向量线性表示。 用命令subspace.来求得两个向量或者两个子空间的夹角。 命令集68 值域、零空间和子空间的夹角 orth(A) 求A空间的正交基，它的列数等于A的秩。 null(A) 求A的零空间的正交基，它的列数等于零空间的维数。 subspace(x,y) 求列向量x和y的夹角，向量的长度必须一样。 subspace(A,B)求由矩阵A和B的列划分的子空间的夹角，列的长度必须一样。 ■例7.2 用命令orth、nulI和subspace求解例7.1中的矩阵的正交基、零空间和夹角。 (a)先求正交基： Range1 orth(A1),Range2 =orth(A2),Range3 orth(A3)

和 结果为： (e) 实数矩阵的行列式和它的转置的行列式相同： 与线性系统相联系的两个子空间是值域和零空间。如果 A为m×n的矩阵，它的秩为 r，那 么A的向量空间就是由A的列划分的线性空间，这个空间的维数是 r，也就是A的秩。如果r=n， 则A的列线性无关。M AT L A B命令o r t h用来求A的空间的正交基。 A的零空间是由满足条件 A x = 0的所有向量x组成的线性子空间。在M AT L A B中可以用命令 n u l l来求得零空间的正交基。 假设有一个向量集v1，v2，. . .vn，可以通过定义矩阵B= (v1 v2 . . . vn )来判断它们是否线性相 关。例如，如果B的秩是n－1，那么其中的一个向量vi可以用其他向量线性表示。 用命令s u b s p a c e来求得两个向量或者两个子空间的夹角。 命令集6 8 值域、零空间和子空间的夹角 o r t h ( A ) 求A空间的正交基，它的列数等于 A的秩。 n u l l ( A ) 求A的零空间的正交基，它的列数等于零空间的维数。 s u b s p a c e ( x , y ) 求列向量x和y的夹角，向量的长度必须一样。 s u b s p a c e ( A , B ) 求由矩阵A和B的列划分的子空间的夹角，列的长度必须一样。 ■ 例7 . 2 用命令o r t h、n u l l和s u b s p a c e求解例7 . 1中的矩阵的正交基、零空间和夹角。 (a) 先求正交基： 第7章 线性方程系统 9 7 下载 ■

98 China-pub.com MATLAB5手册 下载 结果为： Range1 0.5760 0.8174 0.8174 -0.5760 Range2 0.4472 0.8944 Range3 0.3667 -0.9303 0.9303 0.3667 (b)再求它们的秩： ranki rank(orth(A1)),... rank2 rank(orth(A2)),rank3 rank(orth(A3)) 结果为： rank1 2 rank2 1 rank3 2 当然，矩阵的向量空间的秩就等于矩阵本身的秩。 (c)然后求它们的零空间： nullSpace1 null(A1),... nullSpace2 null(A2),nullSpace3 null(A3) 结果为： nullSpace1 Empty matrix:2-by-0 nullSpace2 -0.9487 0.3162 nullSpace3 -0.8729 0.4364 -0.2182 这里的Empty Matrix表示零空间是空的。对AI求它的零空间时结果就得到零向量。 (d)求它们转置矩阵的零空间： nullSpaceT1 null(A1'),... nullSpaceT2 null(A2'),nullSpaceT3 null(A3') 结果为：

结果为： (b) 再求它们的秩： 结果为： 当然，矩阵的向量空间的秩就等于矩阵本身的秩。 (c) 然后求它们的零空间： 结果为： 这里的Empty Matrix表示零空间是空的。对A 1求它的零空间时结果就得到零向量。 (d) 求它们转置矩阵的零空间： 结果为： 9 8 M ATLAB 5 手册 下载

China-pub.com 第7章线性方程系统 下载 nullSpaceT1 Empty matrix:2-by-0 nullSpaceT2 -0.8944 0.4472 nullSpaceT3 Empty matrix:2-by-0 (e)用orth求正交基： RangeT1 orth(A1'),... RangeT2 orth(A2'),RangeT3 orth(A3') 结果为： RangeT1 0.4046 0.9145 0.9145 -0.4046 RangeT2 0.3162 0.9487 RangeT3= 0.2197 -0.4357 0.7508 -0.4958 0.6229 0.7513 (①最后求矩阵的夹角： angle subspace(null(A2),orth(A2')) 结果为： angle 1.5708 ■ 角度为π/2，说明这两个空间是正交的。命令subspace要求列的长度相同。 7.2线性系统的求解和LU因式分解 MATLAB中用运算符1求解线性系统，这个运算符的功能很强大，而且具有智能性。通常， 仔细研究计算过程是有价值的，在MATLAB中有几个专门这样的命令。 令A是nXm的矩阵，b和x是有n个元素的列向量，B和X是n行p列的矩阵。MATLAB用如 下命令求解系统Ax=b: x=A\b 求解更一般的系统AX=B,也用同样的方法，其中B=(b,b,…b): X=A\B 如果A是一个奇异矩阵，或者是近似奇异矩阵，则会给出一个错误信息。 ■例7.3 在MATLAB中求解下列方程组：

(e) 用o r t h求正交基： 结果为： (f) 最后求矩阵的夹角： 结果为： 角度为 / 2，说明这两个空间是正交的。命令 s u b s p a c e要求列的长度相同。 7.2 线性系统的求解和L U因式分解 M AT L A B中用运算符\求解线性系统，这个运算符的功能很强大，而且具有智能性。通常， 仔细研究计算过程是有价值的，在 M AT L A B中有几个专门这样的命令。 令A是n×m的矩阵，b和x是有n个元素的列向量， B和X是n行p列的矩阵。M AT L A B用如 下命令求解系统A x = b： x = A \ b 求解更一般的系统A X = B，也用同样的方法，其中B= (b1 b2 ... bp )： X = A \ B 如果A是一个奇异矩阵，或者是近似奇异矩阵，则会给出一个错误信息。 ■ 例7 . 3 在M AT L A B中求解下列方程组： 第7章 线性方程系统 9 9 下载 ■

100 MATLAB5手册 China-pub.com 下载 2x1+3x2=7 4x1+x2=9 先找出系数矩阵和右边的b: a=(任) b=(g) 解向量是×=(x,x),可以用x=Ab来求得： X 2 ■ MATLAB依据系数矩阵A的不同而相应地使用不同的方法求解线性系统。如果可能， MATLAB先分析矩阵的结构。例如，如果A是对称且正定的，则使用Cholesky分解。 如果没有找到可以替代的方法，则采用高斯消元法和部分主元法。主要是对矩阵进行LU 因式分解或LU分解。这种方法就是令A=LU,其中U是一个上三角矩阵，L是一个带有单位对 角线的下三角矩阵。 然而为了保证计算的稳定性可以使用部分主元法。也就是说，L通常是一个改变序列的下三 角矩阵，即有些行进行互换。这样，L就可能显得结构不完整，将这些排列定义为交换矩阵P。 交换矩阵P的大小为n×n,它实际上是一个单位矩阵，按照交换的顺序来交换列向量。交 换矩阵的逆等于它本身的转置。 PA =LIU LU因式分解可以用非交换下三角矩阵L表示出来： 即交换矩阵L由L=PL给出。 MATLAB中用命令可以求得U和交换或非交换下三角矩陶，后一种情况也可给出交换矩御。 命令集69 LU分解 [L,U]=1u(A) 求上三角矩阵U和交换下三角矩阵L。L是一个带有单位对角线 的下三角矩阵和交换矩阵，即P的逆矩阵的乘积，见下个命令。 [L,U,P]=1u(A) 求上三角矩阵U、有单位对角线的下三角矩阵L和交换矩阵P, 满足LU=PA。 ■例7.4 如果矩阵A=[123;456;780],输入[L,0]=1u(A),运行结果为： L= 0.1429 1.0000 0 0.5714 0.5000 1.0000 1.0000 0 U= 7.0000 8.0000 0 0 0.8571 3.0000 0 0 4.5000

先找出系数矩阵和右边的b： 解向量是x= ( x1 x2 ) ´，可以用x = A \ b来求得： M AT L A B依据系数矩阵 A的不同而相应地使用不同的方法求解线性系统。如果可能， M AT L A B先分析矩阵的结构。例如，如果 A是对称且正定的，则使用C h o l e s k y分解。 如果没有找到可以替代的方法，则采用高斯消元法和部分主元法。主要是对矩阵进行 L U 因式分解或L U分解。这种方法就是令 A = L U，其中U是一个上三角矩阵， L是一个带有单位对 角线的下三角矩阵。 然而为了保证计算的稳定性可以使用部分主元法。也就是说，L通常是一个改变序列的下三 角矩阵，即有些行进行互换。这样，L就可能显得结构不完整，将这些排列定义为交换矩阵P。 交换矩阵P的大小为n×n，它实际上是一个单位矩阵，按照交换的顺序来交换列向量。交 换矩阵的逆等于它本身的转置。 L U因式分解可以用非交换下三角矩阵 Ll表示出来： 即交换矩阵L由L = P ´ Ll给出。 M AT L A B中用命令l u可以求得U和交换或非交换下三角矩阵L，后一种情况也可给出交换矩阵P。 命令集6 9 L U分解 [ L , U ] = l u ( A ) 求上三角矩阵 U和交换下三角矩阵 L。L是一个带有单位对角线 的下三角矩阵和交换矩阵，即 P的逆矩阵的乘积，见下个命令。 [ L , U , P ] = l u ( A ) 求上三角矩阵 U、有单位对角线的下三角矩阵 L和交换矩阵 P， 满足L U = PA。 ■ 例7 . 4 如果矩阵A = [1 2 3 ; 4 5 6 ;7 8 0]，输入[ L , U ] = l u ( A )，运行结果为： 1 0 0 M ATLAB 5 手册 下载 ■

China-bub.com 第7章线性方程系统 101 下载 这里的交换矩阵的逆为： 0 1 P-1=P 0 1 0 0 下面列出高斯消元法的详细步骤： 1 2 3 0 8 0 A= 45 6 A1= 4 5 6 A2= 0 0.4286 6 8 0 1 2 0 0.8571 3 0 7 8 0 A3= 0 0.8571 3 A4= 0 0.8571 3 0 0.4286 6 0 0 4.5 先从矩阵A开始，第1个主元的位置在(1,1)， 通过第1行和第3行互换，可以得到最大的主 元，也就是7。第1次交换后得到矩阵A1,第1次消元后得到矩阵A2。 第2个主元在(2,2)上，将第2行和第3行交换后得到矩阵A3,消元后得到矩阵A4,和矩阵 U是同一个矩阵。 这个过程进行了两次交换，第1次是第1行和第3行互换，也就是： 00 P1= 010 100 第2次进行第2、3行互换，也就是： /1 0 0 P2=00 010 P1和P2相乘，形成P: 00 1 P1*P2=P= 1 00 0 1 0 P的逆为： 0 1 0 P-1= 0 1 1 00 如果输入[L,0,P]=1u(A),运行就会得到如下的结果： L= 1.0000 0 0 0.1429 1.0000 0 0.5714 0.5000 1.0000 U= 7.0000 8.0000 0 0 0.8571 3.0000 0 0 4.5000

这里的交换矩阵的逆为： 下面列出高斯消元法的详细步骤： 先从矩阵A开始，第1个主元的位置在(1, 1)，通过第1行和第3行互换，可以得到最大的主 元，也就是7。第1次交换后得到矩阵A 1，第1次消元后得到矩阵A 2。 第2个主元在(2, 2)上，将第2行和第3行交换后得到矩阵 A 3，消元后得到矩阵A 4，和矩阵 U是同一个矩阵。 这个过程进行了两次交换，第 1次是第1行和第3行互换，也就是： 第2次进行第2、3行互换，也就是： P 1和P 2相乘，形成P： P的逆为： 如果输入[ L , U , P ] = l u ( A )，运行就会得到如下的结果： 第7章 线性方程系统 1 0 1 下载

102 MATLAB5手册 China-bub.c 下载 P= 1 0 0 ■ 有许多方法可求解Ax=b这样的问题，A是n×n的矩阵，b是一个长度为n的列向量。在本 书中没有对算法进行完整的描述，详细信息可参见MATLAB Help Desk。 命令集70 解方程组的方法 x=bicg(A,b,tol, 用双共轭梯度法解方程组。如果给定tol,则用它来 maxit,M) 指定解的精度，也就是和norm(b一A★x)/norm(b) 比较来判断解是否可以接受。如果给定mnaxit,则用 它作为最大的迭代数。要使用预处理因子，可在矩 阵M中规定它。 b1cg(A,b,···,M1,M2,x0) 操作同上，但是使用矩阵M1和M2作为预处理因子。 矩阵M1、M2和预处理因子M的关系为M-M1·M2。 如果给定x0,则将它作为初始化向量开始迭代。 [x,flag,relres,iter, 求x中的问题解，有关bigc的收敛信息存放在f1ag中。 resvect]=bicg(···) 结果的相对剩余范数保存在iter中。值resvect是每次 迭代的范数，结果向量中的所有变量都可忽略。 bicgstab(···) 用稳定双共轭梯度法解方程组，调用方式和返回的 结果形式和命令bicg一样。 cgs(···) 用复共轭梯度平方法解方程组，调用方式和返回的 结果形式和命令bicg一样。 gmres (A, b,restart,用广义最小残量法解方程组，除了参数restarti可以 ·……) 给出外，调用方式和返回的结果形式和命争icg一样。 pcg(···) 用预处理共轭梯度法解方程组，调用方式和返回的 结果形式和命令bicg一样。 qmr(··…) 用准最小残量法解方程组，调用方式和返回的结果 形式和命令bicg一样。 7.3 行梯形矩阵 LU因式分解的另一种方法就是将系数矩阵A降维成行梯形矩阵的形式。因为它能应用在 长方矩阵上，所以这种方法更常用。这种矩阵能给出许多有关线性系统的信息。 对于每一个m×n的矩阵A都有一个交换矩阵P,一个有单位对角线的下三角矩阵L和一个 m×n的梯形矩阵R,它们满足PA=LR。 如果下列情况成立，则矩阵是梯形矩阵： 1)如果有零行，就放在矩阵的底部： 2)每行中第一非零元素是1，它作为每行的最主要元素：

有许多方法可求解 A x = b这样的问题，A是n×n的矩阵，b是一个长度为n的列向量。在本 书中没有对算法进行完整的描述，详细信息可参见 M ATLAB Help Desk。 命令集7 0 解方程组的方法 x = b i c g ( A , b ,t o l, 用双共轭梯度法解方程组。如果给定 t o l，则用它来 m a x i t,M) 指定解的精度，也就是和n o r m ( b－A*x ) / n o r m ( b ) 比较来判断解是否可以接受。如果给定 m a x i t，则用 它作为最大的迭代数。要使用预处理因子，可在矩 阵M中规定它。 b i c g ( A , b ,. . ., M 1 , M 2 ,x 0) 操作同上，但是使用矩阵M 1和M 2作为预处理因子。 矩阵M 1、M 2和预处理因子M的关系为M = M 1·M 2。 如果给定x 0，则将它作为初始化向量开始迭代。 [ x , f l a g , r e l r e s , i t e r , 求x中的问题解，有关b i g c的收敛信息存放在flag中。 r e s v e c t ] =b i c g (. . .) 结果的相对剩余范数保存在i t e r中。值re s v e c t是每次 迭代的范数，结果向量中的所有变量都可忽略。 b i c g s t a b (. . .) 用稳定双共轭梯度法解方程组，调用方式和返回的 结果形式和命令b i c g一样。 c g s (. . .) 用复共轭梯度平方法解方程组，调用方式和返回的 结果形式和命令b i c g一样。 gmres(A, b, restart, 用广义最小残量法解方程组，除了参数 re s t a rt可以 . . .) 给出外，调用方式和返回的结果形式和命令b i c g一样。 p c g (. . .) 用预处理共轭梯度法解方程组，调用方式和返回的 结果形式和命令b i c g一样。 q m r (. . .) 用准最小残量法解方程组，调用方式和返回的结果 形式和命令b i c g一样。 7.3 行梯形矩阵 L U因式分解的另一种方法就是将系数矩阵 A降维成行梯形矩阵的形式。因为它能应用在 长方矩阵上，所以这种方法更常用。这种矩阵能给出许多有关线性系统的信息。 对于每一个m×n的矩阵A都有一个交换矩阵 P，一个有单位对角线的下三角矩阵 L和一个 m×n的梯形矩阵R，它们满足PA=L R。 如果下列情况成立，则矩阵是梯形矩阵： 1) 如果有零行，就放在矩阵的底部； 2) 每行中第一非零元素是1，它作为每行的最主要元素； 1 0 2 M ATLAB 5 手册 下载 ■

China-pub.com 第7章线性方程系统 103 下载 3)每行的最主要元素放在上一行最主要元素的右边： 4)在有最主要元素1的列中，其他元素都是零。 可以用MATLAB中命令rref来分析系统Ax=b。 命令集71 缩减行阶梯矩阵 rref(A) 用高斯一约当消元法和行主元法求A的缩减行的阶梯矩阵。 rref(A,tol) 和refA)一样，但是使用精度tol,它用来决定什么时候元素可 以忽略不计。 rrefmovie(A) 求缩减行的阶梯矩阵，并给出每一步的求解过程。 通常使用的精度是tol=max(size(A))*eps*nom(A,inf)，在7.6节中定义了命令norm。 ■例7.5 令 1 33 2 69 B= 1 -33 0 -221 (a)如果用命令rref(A)得到的矩阵中有一个或多个零行，则可以知道矩阵A中有一些行线 性相关。同样也可用命令rref(A)来知道矩阵A的列向量的相关性。 输入命令Aref=rref(A),Bref=rref(B),运行结果为： Aref 1.0000 3.0000 0 1.0000 0 0 1.0000 0.3333 0 0 0 0 Bref 1 0 0 0 1 0 0 1 (b)求矩阵A的秩可以在A的缩减行的阶梯矩阵中数非零行的数量。在这个例子中，可以 知道A的秩是2，B的秩是3。可以用命令rankA=rank(A),rankB=rank(B)来确认一下： rankA rankB 3 ■ 命令rref可以用来研究线性方程组，令Ax=b代表一个线性方程组，令矩阵B=(Ab)。用 命令C=rref(B)求得矩阵B的缩减行的阶梯矩阵C,然后有下列结论： ·如果C中有一个或多个全零的行向量，则这个方程组中有冗余信息。这就意味着可以去 掉一个或多个方程 ·如果C中有除了最后一个元素不是零外其余都是零的行，即0,0，…，0,1)，则这个方程组无解

3) 每行的最主要元素放在上一行最主要元素的右边； 4 ) 在有最主要元素1的列中，其他元素都是零。 可以用M AT L A B中命令r r e f来分析系统A x = b。 命令集7 1 缩减行阶梯矩阵 r r e f ( A ) 用高斯—约当消元法和行主元法求 A的缩减行的阶梯矩阵。 r r e f ( A , t o l ) 和r r e f ( A )一样，但是使用精度 t o l，它用来决定什么时候元素可 以忽略不计。 r r e f m o v i e ( A ) 求缩减行的阶梯矩阵，并给出每一步的求解过程。 通常使用的精度是t o l = m a x ( s i z e ( A ) )*e p s*norm(A, inf) ，在7 . 6节中定义了命令n o r m。 ■ 例7 . 5 令 (a) 如果用命令r r e f ( A )得到的矩阵中有一个或多个零行，则可以知道矩阵 A中有一些行线 性相关。同样也可用命令r r e f ( A ´ ) 来知道矩阵A的列向量的相关性。 输入命令A r e f = r r e f ( A )，B r e f = r r e f ( B )，运行结果为： (b) 求矩阵A的秩可以在A的缩减行的阶梯矩阵中数非零行的数量。在这个例子中，可以 知道A的秩是2，B的秩是3。可以用命令r a n k A = r a n k ( A )，r a n k B = r a n k ( B )来确认一下： 命令r r e f可以用来研究线性方程组，令 A x = b代表一个线性方程组，令矩阵 B=(A b)。用 命令C = r r e f ( B )求得矩阵B的缩减行的阶梯矩阵C，然后有下列结论： • 如果C中有一个或多个全零的行向量，则这个方程组中有冗余信息。这就意味着可以去 掉一个或多个方程 • 如果C中有除了最后一个元素不是零外其余都是零的行，即(0, 0, . . .,0, 1)，则这个方程组无解 第7章 线性方程系统 1 0 3 下载 ■

104 MATLAB5手册 China-pub.com 下载 ·如果方程组有唯一解，则在C的最后一个列可以得到该解。 7.4 Cholesky因式分解 如果矩阵A是一个对称正定矩阵，即A=A'且对于每个x≠0都有x'Ax>0,则存在一个上 三角矩阵G,它的对角线元素是正数，且满足G'G=A。 这是LU因式分解的一种特殊情况称为Cholesky因式分解，它只有标准LU因式分解的大约 一半的计算步骤。注意，在有些书中是按照下三角矩阵来定义Cholesky因式分解的。 在用左除l来解对称正定的方程组时，MATLAB会自动用Cholesky因式分解来求解。 命令chol可用来计算正定矩阵A的Cholesky因式分解。 命令集72 Choleskyl因式分解 chol(A) 求矩阵A的Cholesky因子，是一个上三角矩阵。如果A 不是一个正定矩阵，则给出一个错误信息。 [G,err]=chol(A) 求矩阵A的Cholesky因子G。如果A不是一个正定矩阵， 则不给出错误信息，而是将err设为非零值。 R1=cholupdate(R,x) 如果R=chol(A)且x是一个和A的列长度一样的列向量， 则求A+xx'的上三角Cholesky因子R1。 R1=cholupdate(R,x,'-')如果R=chol(A)且x是一个和A的列长度一样的列向量， 则求A一xx'的上三角Cholesky因子R1。 ■例7.6 (a)b=[-1-1-1];A=4*eye(4)+diag(b,-1)+diag(b,1),. G=chol(A) 结果为： A= 4 -1 0 0 -1 -1 0 0 -1 4 -1 0 0 -1 4 G= 2.0000 -0.5000 0 0 0 1.9365 -0.5164 0 0 0 1.9322 -0.5175 0 0 1.9319 用Test=G'*G来检验这个结果，得： Test 4.0000 -1.0000 0 0 -1.0000 4.0000 -1.0000 0 0 -1.0000 4.0000 -1.0000 0 0 -1.0000 4.0000 这样就又得到矩阵A

• 如果方程组有唯一解，则在 C的最后一个列可以得到该解。 7.4 Cholesky因式分解 如果矩阵A是一个对称正定矩阵，即 A=A′且对于每个x≠0都有x′A x> 0，则存在一个上 三角矩阵G，它的对角线元素是正数，且满足 G′G=A。 这是L U因式分解的一种特殊情况称为 C h o l e s k y因式分解，它只有标准 L U因式分解的大约 一半的计算步骤。注意，在有些书中是按照下三角矩阵来定义 C h o l e s k y因式分解的。 在用左除\来解对称正定的方程组时， M AT L A B会自动用C h o l e s k y因式分解来求解。 命令c h o l可用来计算正定矩阵A的C h o l e s k y因式分解。 命令集7 2 C h o l e s k y因式分解 c h o l ( A ) 求矩阵A的C h o l e s k y因子，是一个上三角矩阵。如果 A 不是一个正定矩阵，则给出一个错误信息。 [ G , e r r ] = c h o l ( A ) 求矩阵A的C h o l e s k y因子G。如果A不是一个正定矩阵， 则不给出错误信息，而是将 e r r设为非零值。 R 1 = c h o l u p d a t e ( R , x ) 如果R= c h o l ( A )且x是一个和A的列长度一样的列向量， 则求A+x x′的上三角C h o l e s k y因子R 1。 R 1 = c h o l u p d a t e ( R , x ,’-’) 如果R = c h o l ( A )且x是一个和A的列长度一样的列向量， 则求A－x x′的上三角C h o l e s k y因子R 1。 ■ 例7 . 6 结果为： 用Te s t = G ´*G来检验这个结果，得： 这样就又得到矩阵A。 1 0 4 M ATLAB 5 手册 下载