线性方程组

算法介绍

上一篇文章中已经介绍了一种matlab求解线性方程组的方式, 由于线性方程组还可以表示成Ax=bA x=b的矩阵形式, 这里介绍如何通过矩阵的形式来求解线性方程组. 在介绍matlab的函数之前, 先了解一下利用矩阵形式求解方程组的相关算法.

高斯消元法

高斯消元发实际上就是我们初中学习的加减消元法, 通过对整个增广矩阵的一系列初等行变换和列变换, 把系数矩阵转化为行最简矩阵, 从而得出方程组的解向量.

例如方程组{x+2y+z=22x+6y+z=7x+y+4z=3\left\{\begin{array}{r} x+2 y+z=2 \\ 2 x+6 y+z=7 \\ x+y+4 z=3 \end{array}\right.的求解过程大致如下:

[121226171143]\left[\begin{array}{lll|l} 1 & 2 & 1 & 2 \\ 2 & 6 & 1 & 7 \\ 1 & 1 & 4 & 3 \end{array}\right] ==>[1212213131]\left[\begin{array}{ccc|c} 1 & 2 & 1 & 2 \\ & 2 & -1 & 3 \\ & -1 & 3 & 1 \end{array}\right] ==>[12122135/25/2]\left[\begin{array}{ccc|c} 1 & 2 & 1 & 2 \\ & 2 & -1 & 3 \\ & & 5 / 2 & 5 / 2 \end{array}\right]

矩阵LU分解

矩阵的LU分解是把一个矩阵分解成一个上三角矩阵与下三角矩阵乘积的形式. 即 A=L1U\boldsymbol{A}=\boldsymbol{L}^{-1} \boldsymbol{U} , 这样分解后 线性方程组的矩阵形式就可以表示为L1Ux=bL^{-1} U x=b , 解决这个问题我们分两步进行, 首先令y=Uxy=Ux, 将原方程组化为L1y=bL^{-1}y=b的形式, 先解出yy, 再求解xx. 在这种形式下, yy的解是显而易见的, 即 , y=Lby=Lb . 下面介绍如何来对矩阵进行LU分解.

LU分解的矩阵除了是上下三角矩阵外, 还有一个特殊要求, L矩阵的对角线上的值必须都为1, 对U矩阵没有这项要求.

L=[10001]Rm×m\boldsymbol{L}=\left[\begin{array}{ccc} 1 & 0 & 0 \\ \vdots & \ddots & 0 \\ \vdots & \cdots & 1 \end{array}\right] \in \mathfrak{R}^{m \times m}

U=[000]Rm×m\boldsymbol{U}=\left[\begin{array}{ccc} \vdots & \cdots & \vdots \\ 0 & \ddots & \vdots \\ 0 & 0 & \vdots \end{array}\right] \in \mathfrak{R}^{m \times m}

首先把矩阵LU分解换一种表示形式LA=ULA=U . 其中矩阵L又可以分解成一系列的下三角初等矩阵的乘积, 即LmL2L1LA=U\underbrace{L_{m} \ldots L_{2} L_{1}}_{L} A=U , 换一种思路去理解这个式子, 就可以得到LU分解的步骤. 这个式子理解为 对A矩阵进行一系列的初等行变换使其变成一个上三角矩阵 , 经过的行变换的初等矩阵的乘积就是L矩阵, 最后所得到的上三角矩阵就是U矩阵.

例如: A=[111235468]A=\left[\begin{array}{lll} 1 & 1 & 1 \\ 2 & 3 & 5 \\ 4 & 6 & 8 \end{array}\right]

为了把A变成上三角矩阵, 首先利用A(1,1)的元素把A(2,1)和A(3,1)的元素变为零, 所以有:

L1A=[100210401][111235468]=[111013024]\boldsymbol{L}_{1} \boldsymbol{A}=\left[\begin{array}{ccc} 1 & 0 & 0 \\ -2 & 1 & 0 \\ -4 & 0 & 1 \end{array}\right]\left[\begin{array}{ccc} 1 & 1 & 1 \\ 2 & 3 & 5 \\ 4 & 6 & 8 \end{array}\right]=\left[\begin{array}{ccc} 1 & 1 & 1 \\ 0 & 1 & 3 \\ 0 & 2 & 4 \end{array}\right]

这里的L1L_1并不是初等矩阵, 而是两个初等矩阵的乘积.

接下来对矩阵L1AL_1A继续进行行变换, 把第二行的-2倍加到第三行:

L2(L1A)=[100010021][111013024]=[111013002]=U\boldsymbol{L}_{2}\left(\boldsymbol{L}_{1} \boldsymbol{A}\right)=\left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -2 & 1 \end{array}\right]\left[\begin{array}{ccc} 1 & 1 & 1 \\ 0 & 1 & 3 \\ 0 & 2 & 4 \end{array}\right]=\left[\begin{array}{ccc} 1 & 1 & 1 \\ 0 & 1 & 3 \\ 0 & 0 & -2 \end{array}\right]=\boldsymbol{U}

于是L2L1L_2L_1就是矩阵LL, 右侧矩阵则是LU分解中的U矩阵. 分解完成后可以继续求解方程组.

Matlab函数

高斯消元法求解方程组使用rref函数, 其参数为整个方程组的增广矩阵.

例如方程组{x+2y+z=22x+6y+z=7x+y+4z=3[121226171143]\left\{\begin{aligned} x+2 y+z=2 & \\ 2 x+6 y+z=7 & \\ x+y+4 z=3 \end{aligned} \quad \Rightarrow \quad\left[\begin{array}{lll|l} 1 & 2 & 1 & 2 \\ 2 & 6 & 1 & 7 \\ 1 & 1 & 4 & 3 \end{array}\right]\right. 的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=22x+6y+z=7x+y+4z=3\left\{\begin{array}{rr} x+2 y+z= & 2 \\ 2 x+6 y+z= & 7 \\ x+y+4 z= & 3 \end{array}\right. 可以简单的使用如下代码直接求解.

1
2
3
A = [1 2 1;2 6 1;1 1 4];
b = [2; 7; 3];
x = A\b

还可以通过使用inv()函数, 通过逆矩阵的方式x=A1bx=A^{-1}b来求解方程组.

1
2
3
A = [1 2 1;2 6 1;1 1 4];
b = [2; 7; 3];
x = inv(A)*b

但是, 这种方式存在着一定的局限性, 因为不是所有的矩阵都是可逆的, 根据逆矩阵的计算公式A1=1det(A)adj(A)A^{-1}=\frac{1}{\operatorname{det}(A)} \operatorname{adj}(A) 当矩阵A的行列式等于0时, 矩阵就是不可逆的. 当然此处的不可逆是严格的不可逆, 在工程应用中, 有时候矩阵A不是不可逆的, 即A的行列式不是0 ,但却极其接近于零, 这会导致逆矩阵的计算结果对原矩阵很敏感, 原矩阵中一个小小的改动, 会让逆矩阵的计算结果产生很大的变动, 这样的矩阵我们认为它不是一个好的矩阵. 比如矩阵A=[12324.00016987]A = \left[\begin{array}{ccc} 1 & 2 & 3 \\ 2 & 4.0001 & 6 \\ 9 & 8 & 7 \end{array}\right] 和矩阵B=[123256987]B =\left[\begin{array}{lll} 1 & 2 & 3 \\ 2 & 5 & 6 \\ 9 & 8 & 7 \end{array}\right] 相比, 我们更想要的是矩阵B的形式, 矩阵A虽然从数学的角度来看也是一个可以计算的矩阵, 但因为它的行列式过于接近于0, 在计算过程中的误差比较大.

对于线性方程组Ax=bAx=b而言, 衡量系数矩阵是不是一个好的矩阵, 使用下面的公式:

δxx<κ(A)δAA\frac{\|\delta \boldsymbol{x}\|}{\|\boldsymbol{x}\|}<\kappa(\boldsymbol{A}) \frac{\|\delta \boldsymbol{A}\|}{\|\boldsymbol{A}\|}

δA\delta 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)

线性系统

线性系统看起来和线性方程组很类似, 但是其研究对象存在着不同. 线性系统研究的是对于向量进行作用的矩阵的一些性质. 在线性方程组[3214]A[xy]x=[511]b\underbrace{\left[\begin{array}{cc} 3 & -2 \\ 1 & 4 \end{array}\right]}_{A} \underbrace{\left[\begin{array}{c} x \\ y \end{array}\right]}_{x}=\underbrace{\left[\begin{array}{c} 5 \\ 11 \end{array}\right]}_{b} 中我们更关心的是解向量的值, 而在线性系统[21215]A[24]b=[xy]y\underbrace{\left[\begin{array}{cc} 2 & -12 \\ 1 & -5 \end{array}\right]}_{A} \underbrace{\left[\begin{array}{l} 2 \\ 4 \end{array}\right]}_{b}=\underbrace{\left[\begin{array}{l} x \\ y \end{array}\right]}_{y} 中, 我们更关心的是矩阵A的一些性质. 因为在有些情况下, 直接求解y向量是很复杂的, 我们需要寻找一些向量和实数, 使其满足Avi=λivi\boldsymbol{A v}_{i}=\lambda_{i} \boldsymbol{v}_{i} , 这样我们可以把b向量进行分解b=αivi\boldsymbol{b}=\sum \alpha_{i} \boldsymbol{v}_{i} ,于是上面的运算就转化为 Ab=αiAvi=αiλiviA \boldsymbol{b}=\sum \alpha_{i} A v_{i}=\sum \alpha_{i} \lambda_{i} v_{i}

matlab中使用eig()函数来求解矩阵的特征值和特征向量.

A=[21215]A=\left[\begin{array}{cc} 2 & -12 \\ 1 & -5 \end{array}\right]

1
2
[v,d]=eig([2 -12;1 -5])
% v 特征值 d 特征向量

资料链接

参考视频 参考讲义