原文: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之间的相关性。
(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是此欠定方程的一个特解。