『线性方程组 高斯消元算法』

<更新提示>

<第一次更新>


<正文>

线性方程组

形如

\[\begin{cases} a_{11}x_{1}+a_{12}x_{2}+...+a_{1n}x_{n}=b_1 \\ a_{21}x_{1}+a_{22}x_{2}+...+a_{2n}x_{n}=b_2 \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ... \\ a_{m1}x_{1}+a_{m2}x_{2}+...+a_{mn}x_{n}=b_m \end{cases} \]

这样的由\(m\)\(n\)元一次方程组成的方程组我们称为线性方程组。

对于这样的方程,我们通常使用高斯消元算法来求解。

高斯消元

系数矩阵和增广矩阵

首先,我们需要对该方程进行模型转换。

我们可以将每一个方程的系数\(a\)写作一个大小为\(m*n\)的矩阵。

\[\left[ \begin{matrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} \\ \end{matrix} \right] \]

如果把常数\(b\)也写进来,就可以得到该方程组的增广矩阵。

\[\left[ \begin{array}{cccc|c} a_{11} & a_{12} & \cdots & a_{1n} & b_1 \\ a_{21} & a_{22} & \cdots & a_{2n} & b_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} & b_m \\ \end{array} \right] \]

这是高斯消元算法的预备知识,得到了增广矩阵,就方便了我们对方程组进行操作了。

初等行变换

我们知道,解线性方程组一般来说可以用两种方法:代入消元法加减消元法

事实上,在程序中加减消元法会显得更好操作一些,尤其是在我们得到的增广矩阵中。所以,我们定义如下两种矩阵操作:

\(1.\)交换两行的位置
\(2.\)将某一行乘上若干倍加到另一行上

利用这两个操作,我们就能求解该方程,我们称这两种操作为初等行变换,利用增广矩阵和初等行变换求解线性方程组的方法就称为高斯消元算法

高斯消元

解决完预备知识之后,我们考虑如下一种可行性算法:

\(1.\) 枚举\(i\),找到未知数系数\(x_i\)不为零的一个方程
\(2.\) 利用初等行变换\(2.\),用该方程消去其他所有方程中\(x_i\)的系数
\(3.\) 利用初等行变换\(1.\),将该方程放在增广矩阵的第\(i\)

不难发现,当\(i\)被枚举到后,除了第\(i\)个方程以外,其他所有方程\(x_i\)的系数就都是\(0\)了。那么最后我们就能得到这样的一个增广矩阵:

\[\left[ \begin{array}{cccc|c} a'_{1} & 0 & \cdots & 0 & b'_1 \\ 0 & a'_{2} & \cdots & 0 & b'_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & a'_{m} & b'_m \\ \end{array} \right] \]

那么这个方程就已经给出了原方程组的解,\(\forall \ i \ ,x_i=\frac{b'_i}{a'_i}\)

其他问题

我们已经介绍了基础的高斯消元算法,但是在实际运用中可能会出现很多的问题,接下来我们一一解决。

无穷解

在高斯消元的最后,如果找到了全\(0\)方程,就说明该方程与某个\(x\)无关,原方程组有无穷多个解。

利用基础数学知识可以知道,当存在\(p\)个全\(0\)方程时,原方程组有\(n-p\)个主元,\(p\)个自由元。换言之,\(n-p\)个未知数任意取值,都可以计算出与之对应的剩下\(p\)个未知数,满足原方程组。

无解

在高斯消元的第一步中,如果找到了形如\(0=b(b \not = 0)\)的方程,则说明原方程组存在矛盾,方程组无解。

精度

由上可知,高斯消元算法的核心是利用加减消元来消去其他方程的系数,而频繁地使用乘除法难免出现精度误差。

事实上,可以证明,当每一次选取\(x_i\)系数最大的方程进行消元操作时,精度的损失将会降到最小。

所以通常来说,我们每一次选取\(x_i\)系数最大的方程进行消元。

\(Code:\)

inline bool Gauss(void)
{
    for (int i=1;i<=n;i++)
    {
        int row=i;
        for (int j=i+1;j<=n;j++)
            if ( fabs(a[j][i]) > fabs(a[row][i]) )
                row=j;
        //当找不到xi系数非0方程时,必然是无解或无穷解情况之一,直接返回即可
        if ( fabs(a[row][i]) < eps )
            return false;
        if ( row ^ i )swap(a[i],a[row]) , swap(b[i],b[row]);
        for (int j=1;j<=n;j++)
        {
            if (i==j)continue;
            double rate = a[j][i] / a[i][i];
            for (int k=i;k<=n;k++)
                a[j][k] -= a[i][k] * rate;
            b[j] -= b[i] * rate;
        }
    }
    return true;
}

<后记>

posted @ 2019-04-15 20:15  Parsnip  阅读(769)  评论(0编辑  收藏  举报