高斯消元法求解方程组

---

1. 引子

1. 高斯消元法简介

数学上，高斯消元法（Gaussian Elimination），是线性代数规划中的一个算法，可用来为线性方程组求解。

2. 引例

求解如下方程组：

$$\begin{cases} 2x + 3y + z = 4 \\ x - 2y + 4z = -5 \\ 3x + 8y - 2z = 13 \\ \end{cases}$$

2. 求解过程

可以将方程组和矩阵联系起来。如下：

$$\begin{cases} 2x + 3y + z = 11 & (1)\\ x - 2y + 3z = 6 & (2)\\ 3x + 8y - 2z = 13 & (3)\\ \end{cases}$$
$$\begin{pmatrix} 2 & 3 & 1 & 11 \\ 1 & -2 & 3 & 6 \\ \tag {增广矩阵} 3 & 8 & -2 & 13 \\ \end{pmatrix}$$

令 $(1) - 2(2)$、$3 - 3(2)$ 得：

$$\begin{cases} x - 2y + 3z = 6 & (4)\\ 7y - 5z = -1 & (5)\\ 14y - 11z = -5 & (6)\\ \end{cases}$$
$$\begin{pmatrix} 1 & -2 & 3 & 6 \\ 0 & 7 & -5 & -1 \\ 0 & 14 & -11 & -5 \\ \end{pmatrix}$$

令 $(6)-2(7)$ 得：

$$\begin{cases} x - 2y + 3z = 6 & (7) \\ 7y - 5z = -1 & (8) \\ -z = -3 & (9) \\ \end{cases}$$
$$\begin{pmatrix} 1 & -2 & 3 & 6 \\ 0 & 7 & -5 & -1 \\ \tag{行阶梯形矩阵} 0 & 0 & -1 & -3 \\ \end{pmatrix}$$

还可以进一步简化：

$$\begin{cases} x - 2y + 3z = 6 & (7) \\[2ex] y - \cfrac 57 z = -\cfrac 17 & (8) \\[2ex] z = 3 & (9) \\ \end{cases}$$
$$\begin{pmatrix} 1 & -2 & 3 & 6 \\[2ex] 0 & 1 & -\frac57 & -\frac17 \\[2ex] \tag{行最简形矩阵} 0 & 0 & 1 & 3 \\ \end{pmatrix}$$

$(7)+2(8)$、$(8)+\cfrac57(9)$ 继续进行化简：

$$\begin{cases} x = 1 & (7) \\ y = 2 & (8) \\ z = 3 & (9) \\ \end{cases}$$
$$\begin{pmatrix} 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 2 \\ \tag{标准形} 0 & 0 & 1 & 3 \\ \end{pmatrix}$$

可见，当增广矩阵简化为标准形后，解即为矩阵最后一列的值。因此，我们可以脱离原方程组，对矩阵进行化简，得到方程组的解。

3. 编程实现高斯消元法（C#)

    using System.Collections.Generic;
    using System.Linq;

    namespace GaussianEliminationLib
    {
        public static class GaussianElimination
        {
            public static List<double> Resolve(IEnumerable<IEnumerable<double>> equationSet)
            {
                // 存储方程组矩阵
                List<List<double>> equationSetList = new List<List<double>>();
                foreach(IEnumerable<double> equation in equationSet)
                {
                    equationSetList.Add(equation.ToList());
                }

                for (int i = 0; i < equationSetList.Count; ++i)
                {
                    // 化简成行最简形
                    Simplitify(equationSetList, i, i);
                }
                // 化简成标准形
                return SimplitifyAgain(equationSetList);
            }

            private static void Simplitify(List<List<double>> equationSet, int rowStart, int colStart)
            {
                // 选择当前列最大行
                int pivotRow = FindPivotRow(equationSet, rowStart, colStart);
                SwapRows(equationSet, rowStart, pivotRow);

                int cols = equationSet[rowStart].Count;
                double maxValue = equationSet[rowStart][colStart];
                for (int i = rowStart; i < equationSet.Count; ++i)
                {
                    for (int j = colStart; j < cols; ++j)
                    {
                        equationSet[i][j] /= maxValue;
                    }
                }

                for (int i = rowStart + 1; i < equationSet.Count; ++i)
                {
                    double primer = equationSet[i][colStart];
                    for (int j = colStart; j < cols; ++j)
                    {
                        equationSet[i][j] -= primer * equationSet[rowStart][j];
                    }
                }
            }

            private static List<double> SimplitifyAgain(List<List<double>> equationSet)
            {
                int solutionsCount = equationSet.Count;
                List<double> solutions = new List<double>(solutionsCount);
                for (int i = 0; i < equationSet.Count; ++i)
                {
                    solutions.Add(0.0);
                }

                for (int i = equationSet.Count - 1; i >= 0; --i)
                {
                    double value = 0.0;
                    for (int j = 0; j < solutionsCount; ++j)
                    {
                        value += solutions[j] * equationSet[i][j];
                    }
                    solutions[i] = equationSet[i][equationSet.Count] - value;
                }

                return solutions;
            }

            private static int FindPivotRow(List<List<double>> equationSet, int rowStart, int col)
            {
                int row = rowStart;
                double value = equationSet[rowStart][col];
                for (int i = rowStart; i < equationSet.Count(); ++i)
                {
                    if (value < equationSet[i][col])
                    {
                        row = i;
                        value = equationSet[i][col];
                    }
                }
                return row;
            }   

            private static void SwapRows(List<List<double>> equationSet, int currentRow, int pivotRow)
            {
                List<double> row = equationSet[currentRow];
                equationSet[currentRow] = equationSet[pivotRow];
                equationSet[pivotRow] = row;
            }
        }
    }