线性方程组

算法介绍

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

高斯消元法

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

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

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

矩阵LU分解

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

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

$\boldsymbol{L}=\left[\begin{array}{ccc} 1 & 0 & 0 \\ \vdots & \ddots & 0 \\ \vdots & \cdots & 1 \end{array}\right] \in \mathfrak{R}^{m \times 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=U$ . 其中矩阵L又可以分解成一系列的下三角初等矩阵的乘积, 即 $\underbrace{L_{m} \ldots L_{2} L_{1}}_{L} A=U$ , 换一种思路去理解这个式子, 就可以得到LU分解的步骤. 这个式子理解为对A矩阵进行一系列的初等行变换使其变成一个上三角矩阵 , 经过的行变换的初等矩阵的乘积就是L矩阵, 最后所得到的上三角矩阵就是U矩阵.

例如: $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)的元素变为零, 所以有:

$\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]$

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

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

$\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}$

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

Matlab函数

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

例如方程组 $\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()或者使用\(注意除号方向).

对于方程组 $\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=A^{-1}b$ 来求解方程组.

1
2
3

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

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

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

$\frac{\|\delta \boldsymbol{x}\|}{\|\boldsymbol{x}\|}<\kappa(\boldsymbol{A}) \frac{\|\delta \boldsymbol{A}\|}{\|\boldsymbol{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)

线性系统

线性系统看起来和线性方程组很类似, 但是其研究对象存在着不同. 线性系统研究的是对于向量进行作用的矩阵的一些性质. 在线性方程组 $\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}$ 中我们更关心的是解向量的值, 而在线性系统 $\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向量是很复杂的, 我们需要寻找一些向量和实数, 使其满足 $\boldsymbol{A v}_{i}=\lambda_{i} \boldsymbol{v}_{i}$ , 这样我们可以把b向量进行分解 $\boldsymbol{b}=\sum \alpha_{i} \boldsymbol{v}_{i}$ ,于是上面的运算就转化为 $A \boldsymbol{b}=\sum \alpha_{i} A v_{i}=\sum \alpha_{i} \lambda_{i} v_{i}$