高斯消元法求解方程组

引子
- 高斯消元法简介
- 引例
求解过程
编程实现高斯消元法C

1. 引子

1. 高斯消元法简介

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

2. 引例

求解如下方程组：

⎧⎩⎨2x+3y+z=4x−2y+4z=−53x+8y−2z=13

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

2. 求解过程

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

⎧⎩⎨2x+3y+z=11x−2y+3z=63x+8y−2z=13(1)(2)(3)

\begin{cases}2x + 3y + z = 11 & (1)\\x - 2y + 3z = 6 & (2)\\3x + 8y - 2z = 13 & (3)\\\end{cases}

⎛⎝⎜2133−2813−211613⎞⎠⎟(增广矩阵)

\begin{pmatrix}2 & 3 & 1 & 11 \\1 & -2 & 3 & 6 \\ \tag {增广矩阵}3 & 8 & -2 & 13 \\\end{pmatrix}

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

⎧⎩⎨x−2y+3z=67y−5z=−114y−11z=−5(4)(5)(6)

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

⎛⎝⎜100−27143−5−116−1−5⎞⎠⎟

\begin{pmatrix}1 & -2 & 3 & 6 \\0 & 7 & -5 & -1 \\0 & 14 & -11 & -5 \\\end{pmatrix}

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

⎧⎩⎨x−2y+3z=67y−5z=−1−z=−3(7)(8)(9)

\begin{cases}x - 2y + 3z = 6 & (7) \\7y - 5z = -1 & (8) \\-z = -3 & (9) \\\end{cases}

⎛⎝⎜100−2703−5−16−1−3⎞⎠⎟(行阶梯形矩阵)

\begin{pmatrix}1 & -2 & 3 & 6 \\0 & 7 & -5 & -1 \\ \tag{行阶梯形矩阵}0 & 0 & -1 & -3 \\\end{pmatrix}

还可以进一步简化：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪x−2y+3z=6y−57z=−17z=3(7)(8)(9)

\begin{cases}x - 2y + 3z = 6 & (7) \\[2ex]y - \cfrac 57 z = -\cfrac 17 & (8) \\[2ex]z = 3 & (9) \\\end{cases}

⎛⎝⎜⎜⎜⎜⎜⎜100−2103−5716−173⎞⎠⎟⎟⎟⎟⎟⎟(行最简形矩阵)

\begin{pmatrix}1 & -2 & 3 & 6 \\[2ex]0 & 1 & -\frac57 & -\frac17 \\[2ex] \tag{行最简形矩阵}0 & 0 & 1 & 3 \\\end{pmatrix}

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

⎧⎩⎨x=1y=2z=3(7)(8)(9)

\begin{cases}x = 1 & (7) \\y = 2 & (8) \\z = 3 & (9) \\\end{cases}

⎛⎝⎜100010001123⎞⎠⎟(标准形)

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