线性方程组
算法介绍
上一篇文章中已经介绍了一种matlab求解线性方程组的方式, 由于线性方程组还可以表示成Ax=b的矩阵形式, 这里介绍如何通过矩阵的形式来求解线性方程组. 在介绍matlab的函数之前, 先了解一下利用矩阵形式求解方程组的相关算法.
高斯消元法
高斯消元发实际上就是我们初中学习的加减消元法, 通过对整个增广矩阵的一系列初等行变换和列变换, 把系数矩阵转化为行最简矩阵, 从而得出方程组的解向量.
例如方程组⎩⎨⎧x+2y+z=22x+6y+z=7x+y+4z=3的求解过程大致如下:
⎣⎡121261114273⎦⎤ ==>⎣⎡122−11−13231⎦⎤ ==>⎣⎡1221−15/2235/2⎦⎤
矩阵LU分解
矩阵的LU分解是把一个矩阵分解成一个上三角矩阵与下三角矩阵乘积的形式. 即 A=L−1U , 这样分解后 线性方程组的矩阵形式就可以表示为L−1Ux=b , 解决这个问题我们分两步进行, 首先令y=Ux, 将原方程组化为L−1y=b的形式, 先解出y, 再求解x. 在这种形式下, y的解是显而易见的, 即 , y=Lb . 下面介绍如何来对矩阵进行LU分解.
LU分解的矩阵除了是上下三角矩阵外, 还有一个特殊要求, L矩阵的对角线上的值必须都为1, 对U矩阵没有这项要求.
L=⎣⎢⎢⎡1⋮⋮0⋱⋯001⎦⎥⎥⎤∈Rm×m
U=⎣⎢⎢⎢⎡⋮00⋯⋱0⋮⋮⋮⎦⎥⎥⎥⎤∈Rm×m
首先把矩阵LU分解换一种表示形式LA=U . 其中矩阵L又可以分解成一系列的下三角初等矩阵的乘积, 即LLm…L2L1A=U , 换一种思路去理解这个式子, 就可以得到LU分解的步骤. 这个式子理解为 对A矩阵进行一系列的初等行变换使其变成一个上三角矩阵 , 经过的行变换的初等矩阵的乘积就是L矩阵, 最后所得到的上三角矩阵就是U矩阵.
例如: A=⎣⎡124136158⎦⎤
为了把A变成上三角矩阵, 首先利用A(1,1)的元素把A(2,1)和A(3,1)的元素变为零, 所以有:
L1A=⎣⎡1−2−4010001⎦⎤⎣⎡124136158⎦⎤=⎣⎡100112134⎦⎤
这里的L1并不是初等矩阵, 而是两个初等矩阵的乘积.
接下来对矩阵L1A继续进行行变换, 把第二行的-2倍加到第三行:
L2(L1A)=⎣⎡10001−2001⎦⎤⎣⎡100112134⎦⎤=⎣⎡10011013−2⎦⎤=U
于是L2L1就是矩阵L, 右侧矩阵则是LU分解中的U矩阵. 分解完成后可以继续求解方程组.
Matlab函数
高斯消元法求解方程组使用rref
函数, 其参数为整个方程组的增广矩阵.
例如方程组⎩⎪⎨⎪⎧x+2y+z=22x+6y+z=7x+y+4z=3⇒⎣⎡121261114273⎦⎤ 的matlab代码为 :
1 2 3
| A = [1 2 1;2 6 1;1 1 4]; b = [2; 7; 3]; R = rref([A b])
|
矩阵的LU分解使用函数lu()
, 例如:
1 2
| A = [1 1 1;2 3 5;4 6 8]; [L, U, P] = lu(A);
|
下面给出一些matlab中关于矩阵分解的函数的超链接:
实际上, 对于矩阵方程, matlab提供了一个强大的函数可以一步求解, mldivide()
或者使用\
(注意除号方向).
对于方程组⎩⎨⎧x+2y+z=2x+6y+z=x+y+4z=273 可以简单的使用如下代码直接求解.
1 2 3
| A = [1 2 1;2 6 1;1 1 4]; b = [2; 7; 3]; x = A\b
|
还可以通过使用inv()
函数, 通过逆矩阵的方式x=A−1b来求解方程组.
1 2 3
| A = [1 2 1;2 6 1;1 1 4]; b = [2; 7; 3]; x = inv(A)*b
|
但是, 这种方式存在着一定的局限性, 因为不是所有的矩阵都是可逆的, 根据逆矩阵的计算公式A−1=det(A)1adj(A) 当矩阵A的行列式等于0时, 矩阵就是不可逆的. 当然此处的不可逆是严格的不可逆, 在工程应用中, 有时候矩阵A不是不可逆的, 即A的行列式不是0 ,但却极其接近于零, 这会导致逆矩阵的计算结果对原矩阵很敏感, 原矩阵中一个小小的改动, 会让逆矩阵的计算结果产生很大的变动, 这样的矩阵我们认为它不是一个好的矩阵. 比如矩阵A=⎣⎡12924.00018367⎦⎤ 和矩阵B=⎣⎡129258367⎦⎤ 相比, 我们更想要的是矩阵B的形式, 矩阵A虽然从数学的角度来看也是一个可以计算的矩阵, 但因为它的行列式过于接近于0, 在计算过程中的误差比较大.
对于线性方程组Ax=b而言, 衡量系数矩阵是不是一个好的矩阵, 使用下面的公式:
∥x∥∥δx∥<κ(A)∥A∥∥δA∥
δA指矩阵A变化的部分, 寻找一个k使得上式成立, k的值越小 说明矩阵A的变化对结果x的变化影响越小, 此时认为矩阵A是一个好的矩阵. matlab使用cond
函数来计算矩阵的条件数, 条件数的大小反映了矩阵求逆运算的敏感程度.
1 2
| A = [ 1 2 3; 2 4.0001 6; 9 8 7]; cond(A) B = [ 1 2 3; 2 5 6; 9 8 7]; cond(B)
|
线性系统
线性系统看起来和线性方程组很类似, 但是其研究对象存在着不同. 线性系统研究的是对于向量进行作用的矩阵的一些性质. 在线性方程组A[31−24]x[xy]=b[511] 中我们更关心的是解向量的值, 而在线性系统A[21−12−5]b[24]=y[xy] 中, 我们更关心的是矩阵A的一些性质. 因为在有些情况下, 直接求解y向量是很复杂的, 我们需要寻找一些向量和实数, 使其满足Avi=λivi , 这样我们可以把b向量进行分解b=∑αivi ,于是上面的运算就转化为 Ab=∑αiAvi=∑αiλivi
matlab中使用eig()
函数来求解矩阵的特征值和特征向量.
A=[21−12−5]
1 2
| [v,d]=eig([2 -12;1 -5])
|
资料链接
参考视频
参考讲义