加减法(中学所学)是我们平常用的解法之一。 例如，现有如下所示的二元一次方程组。
将等式两边同乘以一个实数，上下系数合并，消去其中一元未知数的方法便是熟知的加减法。
之后，把带入式1，解得。
把上式用行列式表示如下，
之后，第2行乘以，上下相减，得到
的形式。随后，从最下面的式子中解出未知数，代入，得到：
前面的操作是指，把系数行列的左下部分(不包括对角元素)全部变成0，求解如下所示的三角行列(upper triangle matrix)。
这样的处理称作前进消去(forward elimination)。 根据前进消去，未知数可通过解出， 把带入行可以解出，之后把 代入行可以解出 。然后反复执行类似的操作直到解出，这时我们便得到了所有解。这样的后半部分的处理称作后退代入(back substitution)。通过前进消去和后退代入求解n元1次方程组的方法便是高斯消去法。
首先，把常数向量添加到的最后一列组成的行列。
其中第i行j列的元素表示为
由于要在C,C++上实现，这里的元素序号从0开始计数。
前进消去可以表示成下式
这里的
根据上式，把用更新得到上三角行列。这一步的C++代码如下。
// 前进消去(forward elimination) // - 将左下角元素变为0 for(int k = 0; k < n-1; ++k){ double akk = A[k][k]; for(int i = k+1; i < n; ++i){ double aik = A[i][k]; for(int j = k+1; j < n+1; ++j){ A[i][j] = A[i][j]-aik*(A[k][j]/akk); } } }
这里的二维数组A[n][n+1]中存储着系数行列和常数向量。
后退代入的公式如下。
C++代码：
// 后代入(back substitution) // - x_n的解为b_n/a_nn，x_n，代入n-1行求出x_(n-1) A[n-1][n] = A[n-1][n]/A[n-1][n-1]; for(int i = n-2; i >= 0; --i){ double ax = 0.0; for(int j = i+1; j < n; ++j){ ax += A[i][j]*A[j][n]; } A[i][n] = (A[i][n]-ax)/A[i][i]; }
计算结果存储于A[i][n]中。
高斯消去法的完整代码如下。
/*! * 高斯消去法(无主元交换) * @param[inout] A 为n×n的系数项与n×1的常数项(b)合并后n×(n+1)的行列 * @param[in] n n元一次方程组 * @return 1:成功 */ int GaussElimination(vector< vector<double> > &A, int n) { // 前进消去(forward elimination) // - 左下角元素为0 for(int k = 0; k < n-1; ++k){ double akk = A[k][k]; for(int i = k+1; i < n; ++i){ double aik = A[i][k]; for(int j = k+1; j < n+1; ++j){ A[i][j] = A[i][j]-aik*(A[k][j]/akk); } } } // 后退代入(back substitution) A[n-1][n] = A[n-1][n]/A[n-1][n-1]; for(int i = n-2; i >= 0; --i){ double ax = 0.0; for(int j = i+1; j < n; ++j){ ax += A[i][j]*A[j][n]; } A[i][n] = (A[i][n]-ax)/A[i][i]; } return 1; }
上面使用STL的vector代替2维数组。
我们先重新看一遍高斯消去法中的前进消去公式。
会发现右侧的第二项中包含除以的计算。当前进消去过程中对角元素中有0时，则会停止计算。即使不为0，当有绝对值很小的数出现时也会有误差出现。
针对上述问题，这里我们对主元消元法(pivotal elimination method)进行说明。首先，我们知道联立方程的公式的顺序是可以交换的。也就是说，处理第行时，把其对角元素与同列还未处理的元素 的值进行比较，把绝对值最大的那一行与行互换后，应用到前进消去公式中。上述操作便是主元消元。 像这样的只交换行的操作称为部分主元消元(partial pivoting)， 列方向也交换的操作称为完全主元消元(complete pivoting)。值得注意的是，部分主元消元中未知数的顺序不发生变化，完全主元消元中，未知数的顺序会发生变化。
包含主元消元的高斯消去法代码如下
/*! * 主元选择(Pivoting) * - 只交换行的部分主元消去法 * @param[inout] A * @param[in] n n元一次方程组 * @param[in] k 对象行 * @return 1:成功 */ inline void Pivoting(vector< vector<double> > &A, int n, int k) { // 检索k行后k列的绝对值最大元素所在行 int p = k; // 绝对值最大的行 double am = fabs(A[k][k]);// 最大値 for(int i = k+1; i < n; ++i){ if(fabs(A[i][k]) > am){ p = i; am = fabs(A[i][k]); } } // 若k != p，行交换(主元选择) if(k != p) swap(A[k], A[p]); } /*! * 高斯消去法(pivoting) * @param[inout] A * @param[in] n n元一次方程组 * @return 1:成功 */ int GaussEliminationWithPivoting(vector< vector<double> > &A, int n) { // 前进消去(forward elimination) for(int k = 0; k < n-1; ++k){ // 主元消去 Pivoting(A, n, k); double akk = A[k][k]; for(int i = k+1; i < n; ++i){ double aik = A[i][k]; for(int j = k; j < n+1; ++j){ A[i][j] = A[i][j]-aik*(A[k][j]/akk); } } } // 后退代入(back substitution) A[n-1][n] = A[n-1][n]/A[n-1][n-1]; for(int i = n-2; i >= 0; --i){ double ax = 0.0; for(int j = i+1; j < n; ++j){ ax += A[i][j]*A[j][n]; } A[i][n] = (A[i][n]-ax)/A[i][i]; } return 1; }
对于n元1次方程组，用替换，用替换。增广行列可以表示为
由于，两边同乘，得到
也就是说，增广行列的左侧一半变成了单位行列，右侧一半可以得到逆行列。
此时，逆行列为，
像这样的，扩增前进消去，执行上述处理求得逆行列的方法，便是高斯・约旦法。
高斯・约当法公式如下。
求得逆行列代码如下
/*! * 高斯・约当法(无Pivoting) * @param[inout] A n×2n的增广行列 * @param[in] n n元一次方程 * @return 1:成功 */ int GaussJordan(vector< vector<double> > &A, int n) { // 将增广行列的右半部分变为单位行列 for(int i = 0; i < n; ++i){ for(int j = 0; j < n; ++j){ A[i][j+n] = (i == j ? 1.0 : 0.0); } } // 根据Gauss-Jordan method计算逆行列 for(int k = 0; k < n; ++k){ double akk = A[k][k]; // 为使得对角元素为1，k行所有元素除以a_kk for(int j = 0; j < 2*n; ++j){ A[k][j] /= akk; } // k列的非对角元素变为0 for(int i = 0; i < n; ++i){ if(i == k) continue; double aik = A[i][k]; for(int j = 0; j < 2*n; ++j){ A[i][j] -= A[k][j]*aik; } } } return 1; }
