原文:http://www.cnblogs.com/JCSU/articles/1490363.html
				
				 
				
						一、用矩阵运算法解线性方程组
				
				 
				
						1、矩阵运算规则及MATLAB实例
				
				 
				(1)  矩阵加(减)法:即两矩阵对应元素相加减,要求两矩阵的阶数必须相同。
				      C = A + B
				 
				(2)  矩阵乘法:m×p阶矩阵A与p×n阶矩阵B的乘积C是一个m×n阶的矩阵。元素C(i, j)的值为A矩阵的第i行和B矩阵的第j列对应元素乘积的和。
				      A*B = C
				      矩阵乘法一般不满足交换律,即:A*B ≠ B*A
				例:
				
						
						>> A=[1 2 3]
A =
     1     2     3
>> B=[4; 5; 6]
B =
     4
     5
     6
>> A*B 
						%结果是一个标量
						
								
								
ans =
    32
>> B*A 
						%结果是一个矩阵
						
								
								
ans =
     4     8    12
     5    10    15
     6    12    18
>> 
				
				 
				 
				(3)  矩阵求逆:矩阵的逆阵V定义为满足V*A = I (I为与A同阶的单位矩阵)。只有方阵才可以求逆,而且此方阵的行列式必须不为零,det(A) ≠ 0。如果det(A) = 0,求逆时就会出现无穷大,此时的矩阵称为奇异的。
				例:
				
						
						>> A
A =
     4     8    12
     5    10    15
     6    12    18
>> B
B =
     4     8    12
     6    10    15
     7    12     4
>> det(A) 
						%矩阵A的行列式为零
						
								
								
ans =
     0
>> det(B) 
						%矩阵B的行列式不为零
						
								
								
ans =
   112
>> inv(A) 
						%对矩阵A求逆时出现无穷大(Inf)
						
								
Warning: Matrix is singular to working precision.
ans =
   Inf   Inf   Inf
   Inf   Inf   Inf
   Inf   Inf   Inf
>> inv(B) 
						%对矩阵B求逆正常
						
								
								
ans =
   -1.2500    1.0000         0
    0.7232   -0.6071    0.1071
    0.0179    0.0714   -0.0714
>> 
				
				 
				 
				(4)  矩阵转置:矩阵的行列互换后构成转置矩阵。如果A是m×n阶的,则AT是n×m阶的。
				例:
				
						
						>> A
A =
     4     8    12
     5    10    15
>> B=A'
B =
     4     5
     8    10
    12    15
				
				 
				 
				(5)  矩阵分块:矩阵A可以划分成许多小矩阵。
				
						 
				
				例:
				
						
						>> A
A =
     3    -4     3     2    -1
     0    -6     0    -3    -3
     4    -3     4     2    -2
     1     1     1     0    -1
>> A1=A(1:3, 1:3)
A1 =
     3    -4     3
     0    -6     0
     4    -3     4
>> A2=A(1:3, 4:5)
A2 =
     2    -1
    -3    -3
     2    -2
>> A3=A(4, 1:4)
A3 =
     1     1     1     0
>> A4=A(4, 5)
A4 =
    -1
>> 
				
				 
				 
				(5)  矩阵分解为向量:把矩阵沿行向或列向分解为单列或单行向量,常见的是分解为列向量。
				
						 
				
				其中:
				
						 
				
				 
				 
				(6)  行向量左乘列向量:要求两向量的长度必须一致,结果为一个标量。得出的是向量各分量的平方和。如果这些分量是正交的,则得出的是向量的长度平方。它的平方根就是向量长度,也称为向量的范数或2范数。
				行向量左乘列向量得到的是 ,在工程中具有这样计算形式的问题很多,比如用它计算两个向量x和y之间的相关性。
,在工程中具有这样计算形式的问题很多,比如用它计算两个向量x和y之间的相关性。
				 
				(7)  列向量左乘行向量:如果列向量的长度为n,行向量的长度为m,则相乘会得出一个n×m的矩阵,这种方法常常用来生成和计算一些复杂的大矩阵。
				 
				(8)  矩阵乘法和幂次:A^2 = A*A,  A^n =  A*A*...*A。根据矩阵相乘内阶数必须相等,只有方阵才能乘方和幂次。
				 
				 
				
						2、初等变换乘子矩阵的生成及MATLAB实例
				
				 
				本节通过将原矩阵左乘变换的单位矩阵实现矩阵的初等行变换: 
				(1) 行交换:将矩阵A的第 i , j 两行互换位置。
				
						
						
								
						
						% 行交换:交换矩阵A的第1行和第3行
						
								
								
						
						% 步骤1:得到一个5阶单位阵:E=eye(5)
						
								
						
						% 步骤2:交换E的第1行和第3行得到矩阵
						E1
						
								
						
						% 步骤3:矩阵A左乘矩阵E1得到交换后
						的矩阵
						
								
								
>> A
A =
     3    -4     3     2    -1
     0    -6     0    -3    -3
     4    -3     4     2    -2
     1     1     1     0    -1
    -2     6    -2     1     3
>> 
E1 =
     0     0     1     0     0
     0     1     0     0     0
     1     0     0     0     0
     0     0     0     1     0
     0     0     0     0     1
>> E1*A
ans =
     4    -3     4     2    -2
     0    -6     0    -3    -3
     3    -4     3     2    -1
     1     1     1     0    -1
    -2     6    -2     1     3
>> 
				
				 
				写成函数E1gen(A,i,j):
				
						
						function E=E1gen(A,i,j)
n=size(A); 
						%求矩阵A的行数和列数
						
								
m=min(n); 
						%获取矩阵行数和列数中的最小值
						
								
E=eye(m); 
						%产生单位对角阵
						
								
E(i,i)=0; E(j,j)=0; E(i,j)=1; E(j,i)=1;
				
				 
				 
				(2) 行乘数:将矩阵A的第 i 行乘以 k 。
				 
				
						
								
								
								
						
						
								
								 
						
								% 行乘数:将矩阵A的第4行乘以5
								
										
										
								
								% 步骤1:得到一个5阶单位阵:E=eye(5)
								
										
								
								% 步骤2:将单位阵E中的元素E(4,4)乘以5得到矩阵E2
								
										
								
								% 步骤3:矩阵A左乘矩阵E2得到变换后
								的矩阵
								
										
										
>> A
A =
     3    -4     3     2    -1
     0    -6     0    -3    -3
     4    -3     4     2    -2
     1     1     1     0    -1
    -2     6    -2     1     3
>> E2
E2 =
     1     0     0     0     0
     0     1     0     0     0
     0     0     1     0     0
     0     0     0     5     0
     0     0     0     0     1
>> E2*A
ans =
     3    -4     3     2    -1
     0    -6     0    -3    -3
     4    -3     4     2    -2
     5     5     5     0    -5
    -2     6    -2     1     3
>> 
						
				 
				 
				写成函数E2gen(A,i,k):
				 
				
						
						function E=E2gen(A,i,k)
n=size(A);
m=min(n);%获取矩阵行数和列数中的最小值
E=eye(m);%产生单位对角阵
E(i,i)=k;
				
				 
				 
				(3) 行相加:将矩阵的第 i 行乘以 k,加到第 j 行。
				
						
						
								
						
						% 行乘数加到另一行:将矩阵A的第4行乘以2加到第5行
						
								
								
						
						% 步骤1:得到一个5阶单位阵:E=eye(5)
						
								
						
						% 步骤2:将单位阵E中的元素E(5,4)乘以2得到矩阵E3
						
								
						
						% 步骤3:矩阵A左乘矩阵E3得到变换后的矩阵
						
								
								
>> A
A =
     3    -4     3     2    -1
     0    -6     0    -3    -3
     4    -3     4     2    -2
     1     1     1     0    -1
    -2     6    -2     1     3
>> E3
E3 =
     1     0     0     0     0
     0     1     0     0     0
     0     0     1     0     0
     0     0     0     1     0
     0     0     0     2     1
>> E3*A
ans =
     3    -4     3     2    -1
     0    -6     0    -3    -3
     4    -3     4     2    -2
     1     1     1     0    -1
     0     8     0     1     1
>> 
				
				 
				写成函数E3gen(A,i,j,k):
				 
				
						
						function E=E3gen(A,i,j,k)
n=size(A);
m=min(n);%获取矩阵行数和列数中的最小值
E=eye(m);%产生单位对角阵
E(j,i)=k;
				
				 
				
						实例:
				
				要消去下列矩阵的A(2,1),求乘子矩阵E3
				
						
						A =
     3    -4     3     2    -1
     6    -6     0    -3    -3
     4    -3     4     2    -2
     1     1     1     0    -1
    -2     6    -2     1     3
				
				 
				在第二行加以第一行乘 -A(2,1)/A(1,1) = -2,故令E3 = E3gen(A,1,2,-2)
				
						
						>> E3=E3gen(A,1,2,-2)
E3 =
     1     0     0     0     0
    -2     1     0     0     0
     0     0     1     0     0
     0     0     0     1     0
     0     0     0     0     1
>> E3*A
ans =
     3    -4     3     2    -1
     0     2    -6    -7    -1
     4    -3     4     2    -2
     1     1     1     0    -1
    -2     6    -2     1     3
>> 
				
				 
				
						综上可知:
				
				
						1、
						行阶梯生成等价于矩阵左乘
				
				因此,整个行阶梯形式U的生成过程,可以看作把原矩阵左乘以一系列的初等变换矩阵E1和E3。把这些初等矩阵的连乘积写成EX,则有: U=EX*A,设其逆为L:
				从而有 L*U=A 
				就是说,A可以分解为一个准下三角矩阵L和一个上三角(即行阶梯)矩阵U的乘积。MATLAB提供了三角分解的函数lu,它的调用方法是:[L,U]=lu(A) 
				
						2、
						lu分解是求行阶梯的一个方法
				
用lu函数求出的U实际上就是A的行阶梯形式(不是简化行阶梯形式)。所以,求简化行阶梯形式用rref函数,而求行阶梯形式可以用lu 函数。不过,它和我们用消元运算所得U的数据不一定相同,尽管得出的阶次和阶梯形状相同。但因为行阶梯形式可以有无数种,用不同步骤算出的结果也不同。只有变成简化行阶梯形式,才能进行比较,看它是不是惟一的。 
 
 
3、行列式
(1) 行列式的几何意义:
行列式的几何意义是面积或体积,它的用途很单一,就是判断奇异性,连正负号都不必关心。
(2) 行列式的计算方法:
计算行列式的最好方法还是行阶梯法,可以利用lu分解
                [L,U] = lu(A)
把A分解为一个准下三角矩阵L和一个上三角矩阵U的乘积。因为det(L) = 1,所以U和A的行列式相等。  
               det(A) = det(U)
而三角矩阵U的det(U)很好求。只要把U的主对角线元素连乘就可得到它的行列式。
 
实例:
>> A
A =
    10     8     6     4     1
     2     5     8     9     4
     6     0     9     9     8
     5     8     7     4     0
     9     4     2     9     1
>> [L U]=lu(A)
L =
    1.0000            0           0           0            0
    0.2000   -0.7083    1.0000           0            0
    0.6000    1.0000            0           0            0
    0.5000   -0.8333    0.8000   -0.2953    1.0000
    0.9000    0.6667   -0.6588    1.0000            0
U =
   10.0000    8.0000     6.0000     4.0000    1.0000
            0   -4.8000     5.4000     6.6000    7.4000
            0            0   10.6250    12.8750    9.0417
            0            0            0      9.4824    1.1235
            0            0            0             0   -1.2349
>> du=diag(U)
du =
   10.0000
   -4.8000
   10.6250
    9.4824
   -1.2349
>> result=prod(du)
result =
  5.9720e+003
>> 
 
用det求A的行列式值得到相同的结果:
>> det(A)
ans =
        5972
>> 
 
 
 
4、矩阵的秩和矩阵求逆
 
矩阵的秩:
(1) 按定义,矩阵的秩是矩阵A中行列式不等于零的最高阶子式的阶次。是用以衡量联立方程中有效方程数目的指数。
(2) 按照定义来计算矩阵的秩,可能遇到的问题也是子矩阵的数量很大,每个矩阵的行列式计算又非常麻烦,其计算量也将是不可接受的天文数字。
(3) 计算矩阵的秩的最好方法是行阶梯法,行阶梯化简后非全为零的行数就是该矩阵的秩。用MATLAB函数r=rank(A)可以检验A的秩,rank函数对A是否是方阵没有要求,即可以有m≠n。
 
矩阵求逆:
(1) 对于n×n方阵A,当r=n时,称A是满秩的,若r<n,必有det(A) = 0,称A是欠秩的或奇异的。奇异矩阵不可以求逆。  
(2) 矩阵求逆的最简单方法也是行阶梯化简,其方法是设定一个由A和I组成的增广矩阵C = [A,I],求C的简化行阶梯形式UC = rref([A,I]),得出UC = [I,V]。V就显示出这个逆矩阵的内容。 
 
求逆实例:
>> A=[3 0 3 -6;5 -1 1 -5;-3 1 4 -9;1 -3 4 -4]
A =
     3     0     3    -6
     5    -1     1    -5
    -3     1     4    -9
     1    -3     4    -4
>> C=[A,eye(4)]
C =
     3     0     3    -6     1     0     0     0
     5    -1     1    -5     0     1     0     0
    -3     1     4    -9     0     0     1     0
     1    -3     4    -4     0     0     0     1
>> UC=rref(C) %UC右边四列就是矩阵A的逆
UC =
    1.0000         0         0         0    0.2323   -0.0101   -0.1313   -0.0404
         0    1.0000         0         0    0.5354   -0.3131   -0.0707   -0.2525
         0         0    1.0000         0    0.5859   -0.4747   -0.1717    0.1010
         0         0         0    1.0000    0.2424   -0.2424   -0.1515    0.0303
>> V=UC(:,5:8)
V =
    0.2323   -0.0101   -0.1313   -0.0404
    0.5354   -0.3131   -0.0707   -0.2525
    0.5859   -0.4747   -0.1717    0.1010
    0.2424   -0.2424   -0.1515    0.0303
>> 
 
矩阵求逆命令:V=inv(A)
 
 
5、条件数—衡量奇异程度的量
(1) 在用数值方法计算矩阵的逆时,由于计算中的误差,人们不大可能得到理想的零合理想的全零行,所以矩阵是否奇异,并不是那么绝对的。
(2) 为了评价矩阵接近‘奇异’的程度,采用了‘条件数’(Condition Number)作为常用的衡量指标。它永远大于1。其数值愈接近于1,计算误差愈小。
(3) MATLAB中,条件数用cond(A)计算,它达到10的4次方以上时,求逆的误差就可能相当可观。像现在,条件数达到10的16次方(注:条件数是逆条件数RCOND的倒数),结果是根本不能用的。
 
6、用矩阵‘除法’解线性方程
 
(1) 如果m = n,则线性代数方程Ax = b中的A是方阵,设det(A)≠0,则它的逆阵存在。将上式左右同乘以inv(A) ,由于inv(A)*A = I,得到
X = inv(A)*b (1)
MATLAB创立了矩阵除法的概念,因为 inv(A)相当于将A放到分母上去,所以可以把上式写成
X = A \ b      (2)
‘\’就称为左除,因为inv(A)是乘在b的左方。
 
(2) 左除’\’解线性方程的扩展
左除‘\’的功能远远超过了矩阵求逆函数inv,inv(A)函数要求A必须是方阵,所以(1)式只能用来解‘适定’方程,而(2)式并不要求A为方阵,在A是m×n阶且m < n(欠定)时,它只要求A与b的行数相等且A的秩为m。所以(2)式也可以用来解欠定方程,在下例中可以看出。此外,运算符‘\’还能用来解超定方程。
 
例:左除’\’解欠定方程
>> A=[3,4,3,2,1;0,6,0,3,3;4,3,4,2,2;1,1,1,0,1;2,6,2,1,3]
A =
     3     4     3     2     1
     0     6     0     3     3
     4     3     4     2     2
     1     1     1     0     1
     2     6     2     1     3
>> b=[2;3;2;0;1]
b =
     2
     3
     2
     0
     1
>>  x=A\b
Warning: Matrix is singular to working precision.
x =
       NaN
       NaN
       Inf
    1.0000
    0.0000
>> 
 
得到x=inf,无解。改用行阶梯方法找有效行。左除要求的是系数矩阵的行数与秩相同
>> A
A =
     3     4     3     2     1
     0     6     0     3     3
     4     3     4     2     2
     1     1     1     0     1
     2     6     2     1     3
>> b
b =
     2
     3
     2
     0
     1
>> B=[A,b]
B =
     3     4     3     2     1     2
     0     6     0     3     3     3
     4     3     4     2     2     2
     1     1     1     0     1     0
     2     6     2     1     3     1
>> r=rank(B)
r =
     4
>> [UB,ip]=rref(B)
UB =
     1     0     1     0     0     0
     0     1     0     0     0     0
     0     0     0     1     0     1
     0     0     0     0     1     0
     0     0     0     0     0     0
ip =
     1     2     4     5
>> U0=UB(1:r,1:5)
U0 =
     1     0     1     0     0
     0     1     0     0     0
     0     0     0     1     0
     0     0     0     0     1
>> d=UB(1:4,6)
d =
     0
     0
     1
     0
>> x=U0\d
x =
     0
     0
     0
     1
     0
>> 
 
其中,x是此欠定方程的一个特解。